myscratchpad

# 2件法の01データの場合はKuder-Richardsonの公式というのを使うらしい。全く知らなかった。知らないことが多いものである。日本語でぐぐると99件しかヒットしない。皆知らないのだろう。
 # 最近4件法や5件法は順序尺度として扱うべしとの意見があるが、信頼性を出すときはどうするんだろう。クロンバックでいいのだろうか
 # このサイトに少し解説。以下適当な訳と覚書
  # 1. 多くの研究でそうしているように、クロンバックのアルファを使う。その仮定が満たされている限り
  # 2. ポリコリック相関行列からクロンバックのアルファ的なものを計算する (B D Zumbo, A M Gadermann, and C Zeisser. Ordinal versions of coefficients alpha and theta for likert rating scales. Journal of Modern Applied Statistical Methods, 6: 21-29, 2007.)
 # 3. 01データならKuder-Richardsonの公式を使う
 # 4. 他にもいろんな信頼性指標がある (R E Zinbarg, W Revelle, I Yovel, and W Li. Cronbach’s α, Revelle’s β, and McDonald’s ωh: Their relations with each other and two alternative conceptualizations of reliability. Psychometrika, 70(1): 123-133, 2005.)
 # A good review is Bruce Thompson. Score Reliability. Contemporary Thinking on Reliability issues. Sage Publications, 2003. だそうだ。まあとりあえずはクロンバックを使おう。気が向いたらRevelle先生のやつを読もう
# 関連記事：Rで級内相関

# Kuder-Richardsonの公式
 # 以下参考URL。感謝
 # http://www.eflclub.com/elvin/publications/2003/itemanalysis.html
 # http://www.gifted.uconn.edu/siegle/research/Instrument%20Reliability%20and%20Validity/Reliability.htm
 # http://www.jstage.jst.go.jp/article/rika/17/3/17_181/_article/-char/ja

 dat <- data.frame(
 q001 = c(1,1,0,0,1,0,0,0,0),
 q002 = c(1,1,0,1,1,1,0,0,0),
 q003 = c(1,1,1,0,1,0,0,0,0),
 q004 = c(1,1,1,1,1,0,0,0,0),
 q005 = c(1,1,1,1,1,1,0,1,0),
 q006 = c(1,0,1,1,1,0,0,1,0),
 q007 = c(1,1,1,0,0,0,1,0,0),
 q008 = c(1,1,1,1,0,1,1,0,0),
 q009 = c(1,1,1,1,0,0,0,0,0),
 q010 = c(1,1,0,0,0,0,0,0,0)
 )
 dat

 (ni <- ncol(dat)) # 項目数
 (ns <- nrow(dat)) # 参加者数
 (mt <- mean(rowSums(dat))) # 合計得点の平均点
 (sdt <- sd(rowSums(dat))) # 合計得点の標準偏差
 (sds <- sqrt(var(rowSums(dat))*((ns-1)/ns))) # 合計得点の標本標準偏差

 pc <- colSums(dat)/ns # 各項目ごとの正答率 (1の個数/全参加者数)
 pw <- 1-pc # 各項目ごとの誤答率 (0の個数/全参加者数)
 pcpw <- pc*pw
 spcpw <- sum(pcpw)


# Kuder-Richardson 20
 (ni/(ni-1)) * (((sds*sds)-spcpw)/(sdt*sdt)) # 不偏標準偏差を使ったとき
 (ni/(ni-1)) * (((sds*sds)-spcpw)/(sds*sds)) # 標本標準偏差を使ったとき

# Kuder-Richardson 21
# おそらく、ローデータがないときに使う簡易版
ni/(ni-1) * (1-(mt-mt*mt/ni)/(sdt*sdt)) # 不偏標準偏差を使ったとき
ni/(ni-1) * (1-(mt-mt*mt/ni)/(sds*sds)) # 標本標準偏差を使ったとき

# よくわからんけど標本標準偏差を使うのかな



# おまけ。カッパ係数
 # 参考:http://www.mizumot.com/stats/kappa.htm
a <- c(2,1,3,2,2,3,2,1,3,3,3,3,2,1,2,3,1,3,1,2)
b <- c(2,1,3,2,3,3,3,1,3,3,3,3,2,1,2,3,2,3,1,2)
x <- table(a,b) # データフレームではなくtableで集計表にする
x

library(psych)
print.digits(cohen.kappa(x), digits=5)
 # 普通はunweighted kappaのestimateを見る

# psychパッケージ以外。まあpsychパッケージが一番いいけど
library(vcd)
Kappa(x)

cprp <- function(x, p=1) {
xn <- addmargins(x)
cp <- t(t(xn)/xn[nrow(xn),]*p)[-nrow(xn), -ncol(xn)]
rp <- (xn/xn[, ncol(xn)]*p)[-nrow(xn), -ncol(xn)]
pt <- prop.table(x)
return(list(adm=xn, cp=cp, rp=rp, pt=pt))
}

x <- matrix(c(10,10,20, 15, 5,20, 25,10,5, 5,5,10), nr=3,
dimnames=list(paste("r", 1:3, sep=""), paste("c", 1:4, sep="")))
x
cprp(x)

x <- rnorm(1000)

## 正規性の検定
hist(x)
qqnorm(x)
qqline(x)
shapiro.test(x) # p > .05なら正規分布

ks.test(x=x, y="pnorm", mean=mean(x), sd=sd(x))

## ks.testは「2群の分布が等しい」という帰無仮説を検定するもの。正規性の検定の場合、yに正規分布を指定している
x2 <- rexp(1000) # 正規分布じゃないやつ
ks.test(x=x, y=x2)


## 多変量正規性
library(mvnormtest)
data(EuStockMarkets)
C <- t(EuStockMarkets[15:29,1:4])
mshapiro.test(C)

下のspssのページをもとにした。
http://www.nyu.edu/its/statistics/Docs/intracls.html
下のページも参照。感謝
http://www.hs.hirosaki-u.ac.jp/~pteiki/research/stat/S/icc/

dat <- data.frame(
rater1=c(9,6,8,7,10,6),
rater2=c(2,1,4,1,5,2),
rater3=c(5,3,6,2,6,4),
rater4=c(8,2,8,6,9,7)
)

# psycパッケージのICC関数を使う。psychパッケージはマジ便利
library(psych)
ICC(dat)
## 勉強のため、ICCのコードの一部を貼り付けておく
x <- dat
n.obs <- dim(x)[1]
nj <- dim(x)[2]
x.s <- stack(x)
x.df <- data.frame(x.s, subs = rep(paste("S", 1:n.obs, sep = ""),
nj))
aov.x <- aov(values ~ subs + ind, data = x.df)
s.aov <- summary(aov.x)
stats <- matrix(unlist(s.aov), ncol = 3, byrow = TRUE)
MSB <- stats[3, 1]
MSW <- (stats[2, 2] + stats[2, 3])/(stats[1, 2] + stats[1,
3])
MSJ <- stats[3, 2]
MSE <- stats[3, 3]
(ICC1 <- (MSB - MSW)/(MSB + (nj - 1) * MSW)) # ICC1. Single raters absolute
(ICC2 <- (MSB - MSE)/(MSB + (nj - 1) * MSE + nj * (MSJ - MSE)/n.obs)) # ICC2. Single random raters
(ICC3 <- (MSB - MSE)/(MSB + (nj - 1) * MSE)) # ICC3. Single fixed raters
(ICC12 <- (MSB - MSW)/(MSB)) # ICC1k. Average random raters
(ICC22 <- (MSB - MSE)/(MSB + (MSJ - MSE)/n.obs)) # ICC2k. Average random raters
(ICC32 <- (MSB - MSE)/MSB) # ICC3k. Average fixed raters


head(x.df)
aov.x
s.aov
round(stats, 3)
alpha(dat) # 質問紙の信頼性係数を求めるときと同じだが、列には項目ではなく評定者がならぶ

gplotパッケージのplotmeans関数を使ってみよう
こちらにも解説がある

もしくはpsychパッケージのerror.bars.by関数か

凡例表示を覚えないとどうしようもないが