#################################################
#豊田秀樹(1994) 違いを見ぬく統計学　実験計画と分散分析入門　講談社
#第5章 p.148-152
#################################################

##設定
source("http://blue.zero.jp/yokumura/R/testtheory/ctt.txt")
library(lme4)
library(plotrix)


##データの読み込み 
dat <- read.table("http://blue.zero.jp/yokumura/R/testtheory/skate.dat", header=T)
dat


##3因子変量モデルの分散成分
mod1 <- lmer(value ~ 1 + (1|judge) + (1|subj) + (1|item)+
             (1|subj:judge) + (1|subj:item) + (1|judge:item), data=dat)
mod1


##分散成分の抽出
varcomp <- VarCorr(mod1)

output <- rep(NA, length(varcomp)+1)
for(i in 1:length(varcomp)){
    output[i] <- varcomp[[i]][,1]
}
names(output) <- c(names(varcomp), "Residual")
output <- sort(output, na.last=T, decreasing=T)
output["Residual"] <- attr(varcomp, "sc")[[1]]^2


##信頼性の推定
var.a  <- output[["subj"]]
var.ab <- output[["subj:judge"]]
var.ac <- output[["subj:item"]]
var.e  <- output[["Residual"]]
  
reliability1 <- function(b, c){
  denom <- var.a + var.ab / b + var.ac / c + var.e / (b * c)
  alpha1  <- var.a / denom
  return(alpha1)
}
reliability1(b=9, c=1)    #「えこひいき」「得手不得手」も誤差の内 
                          #実力を出す機会1回，審査委員9人 

reliability2 <- function(b, c){
  denom <- var.a + var.ab / b + var.ac / c + var.e / (b * c)
  alpha2  <- (var.a + var.ab / b)/ denom
  return(alpha2)
}
reliability2(b=9, c=1)    #「えこひいき」も実力の内 
                          #実力を出す機会1回，審査委員9人 

reliability3 <- function(b, c){
  denom <- var.a + var.ab / b + var.ac / c + var.e / (b * c)
  alpha3  <- (var.a + var.ac / c)/ denom
  return(alpha3)
}
reliability3(b=9, c=1)    #「得手不得手」も実力の内 
                          #実力を出す機会1回，審査委員9人 

reliability3(b=9, c=2)    #「得手不得手」も実力の内 
                          #実力を出す機会2回，審査委員9人 


##オリジナル関数を使用
# 3因子変量モデルを実行して，信頼性係数を算出する 
#------------------------------------------------------------------
#DG.2facets(data)
# dataには，"value" "subj" "judge" "item" の4つの変数名が必要
# subjが対象，itemとjudgeが相，valueが従属変数 
#------------------------------------------------------------------
DG.2facets(data=dat)
graphics.off()

DG.2facets(data=dat, 
           sim.item = seq(1, 10, by=1), 
           sim.judge = seq(1, 12, by=1))  #観点が1-10, 評定者が1-12に変化させたときの信頼性を推定 
graphics.off()

DG.2facets(data=dat, 
           sim.item = seq(1, 20, by=1), 
           sim.judge = seq(1, 5, by=1))  #観点が1-10, 評定者が1-5に変化させたときの信頼性を推定 
graphics.off()