手計算でもaov
でも一元配置分散分析は出来るので、Rで実装してみようと思った次第です。F分布の累積分布を除いてほぼほぼベクトル演算?を使っています。データは水準の繰り返し数のトータルから適当に生成しています。統計学的観点はほぼ0です。
実装
# データ設定 #### ## 水準 catv <- c( rep("A",5),rep("B",5),rep("C",10),rep("D",3)) ## データの範囲 data.range <- 10:30 ## データ datav <- sample(data.range,length(catv),T) # 計算 #### ## 総平均 mean.all <- mean(datav) ## 水準ごとの平均 mean.cat <- tapply(datav, catv, mean) ## 水準平方和 catv.sumsq <- sum( table(catv) * (mean.cat - mean.all)^2 ) ## 自由度: 水準の種類数 - 1 catv.df <- length(table(catv))-1 ## 全平方和 total.sumsq <- sum((datav - mean.all)^2) ## 残差平方和 ### 残差平方和 := 全平方和 - 水準平方和 zansa.sumsq <- total.sumsq - catv.sumsq ## 自由度: 水準ごとの繰り返し数-1 をそれぞれ足す zansa.df <- sum(table(catv)-1) ### 別解: 残差平方和 := 要素 - 水準平均 の平方和 #abs( zansa.sumsq - sum( (datav - rep(mean.cat,table(catv)))^2 ) ) < 0.01 # 答え合わせ #### # aov summary( aov(datav ~ catv) ) # 自力で計算 cat("\t ","Df","Sum Sq"," Mean Sq","F value","Pr(>F)","\n" ," catv"," ",catv.df,"",round(catv.sumsq,2), round(catv.sumsq/catv.df,2) ," ",round((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),3)," ",round(1-pf((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),catv.df,zansa.df ),3),"\n" ," Residuals",zansa.df,"",round(zansa.sumsq,2), "",round(zansa.sumsq/zansa.df,2) )
実行例
> # 答え合わせ #### > # aov > summary( aov(datav ~ catv) ) Df Sum Sq Mean Sq F value Pr(>F) catv 3 170.9 56.96 2.349 0.105 Residuals 19 460.8 24.25 > > # 自力で計算 > cat("\t ","Df","Sum Sq"," Mean Sq","F value","Pr(>F)","\n" + ," catv"," ",catv.df,"",round(catv.sumsq,2), round(catv.sumsq/catv.df,2) + ," ",round((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),3)," ",round(1-pf((catv.sumsq/catv.df)/(zansa.sumsq/zansa.df),catv.df,zansa.df ),3),"\n" + ," Residuals",zansa.df,"",round(zansa.sumsq,2), "",round(zansa.sumsq/zansa.df,2) + ) Df Sum Sq Mean Sq F value Pr(>F) catv 3 170.89 56.96 2.349 0.105 Residuals 19 460.77 24.25
学び
table関数の出力結果はベクトルのようである。
確信はないのだけれどtable関数の出力はベクトルっぽく扱えます。
例
> table( c(rep("A",2),rep("B",4)) ) A B 2 4 > table( c(rep("A",2),rep("B",4)) ) +3 A B 5 7
rep関数の繰り返し回数はベクトルを渡すことができる
水準ごとの平均を水準内の繰り返し数分だけ繰り返す、といった処理をしたいときにrepとtableを組み合わせて対処することができます。
例
10を1回、20を2回、30を3回、40を4回繰り返す。
> rep(c(10,20,30,40),times = c(1,2,3,4)) [1] 10 20 20 30 30 30 40 40 40 40