# 組み込みデータを使って二元配置分散分析 # 二元配置の平均と標準誤差または標準偏差の折れ線グラフを描く関数 twowayplot <- function(VAL, F1, F2, DATASET, TYPE=1, MAIN="") { if (nchar(MAIN)<1) { MAINA <- sprintf("Two way plot of means with %s", ifelse(TYPE==1, "se", "sd"))} else { MAINA <- MAIN} MM <- by(DATASET[, deparse(substitute(VAL))], list(DATASET[, deparse(substitute(F1))], DATASET[, deparse(substitute(F2))]), mean, na.rm=TRUE) SD <- by(DATASET[, deparse(substitute(VAL))], list(DATASET[, deparse(substitute(F1))], DATASET[, deparse(substitute(F2))]), sd, na.rm=TRUE) NN <- by(DATASET[, deparse(substitute(VAL))], list(DATASET[, deparse(substitute(F1))], DATASET[, deparse(substitute(F2))]), length) SE <- SD/sqrt(NN) XLAB <- levels(DATASET[, deparse(substitute(F1))]) CLAB <- levels(DATASET[, deparse(substitute(F2))]) matplot(MM, pch=16, type="b", col=1:length(CLAB), xlab=deparse(substitute(F1)), ylab=deparse(substitute(VAL)), ylim=c(min(MM-SD, na.rm=TRUE), max(MM+SD, na.rm=TRUE)), lty=1:length(CLAB), lwd=1, axes=FALSE, main=MAINA) if (TYPE==1) { arrows(1:length(XLAB), MM-SE, 1:length(XLAB), MM+SE, angle=90, code=3, length=0.1) } else { arrows(1:length(XLAB), MM-SD, 1:length(XLAB), MM+SD, angle=90, code=3, length=0.1) } axis(1, 1:length(XLAB), XLAB) axis(2) legend("topright", pch=16, lty=1:length(CLAB), lwd=1, col=1:length(CLAB), title=deparse(substitute(F2)), legend=CLAB) return(list(MEAN=MM, SD=SD, SE=SE, N=NN)) } # データはMASSパッケージから。 library(MASS) # タイプIIIの平方和は,carパッケージのAnova()でtype=3を付ける library(car) # タイプIIIの平方和の結果をSASなどと一致させるために必要な設定 options(contrasts = c("contr.sum", "contr.sum")) # グラフはフォント入りでpdfへ出力 cairo_pdf("twowayanova-samples.pdf", width=10, height=10, onefile=TRUE, family="Meiryo") # MASSパッケージのcabbages # ソース:http://web.nchu.edu.tw/~numerical/course1012/ra/Applied_Regression_Analysis_A_Research_Tool.pdf # アウトカム:収穫したキャベツの重さ(HeadWt, kg) # 要因1:キャベツの品種(Cult, c39/c52) # 要因2:植え付けた日(Date, d16/d20/d21) twowayplot(HeadWt, Cult, Date, cabbages) Anova(lm(HeadWt ~ Cult*Date, data=cabbages), type=3) # 主効果,交互作用効果とも有意な例 anova(lm(HeadWt ~ Cult*Date, data=cabbages)) # サイズ不斉でタイプ3でないので結果が異なり,これは不適切 # MASSパッケージのcrabs # 西オーストラリアの蟹の形態 # ソース:http://duch.mimuw.edu.pl/~pokar/StatystykaII/DANE/CampbellMahon.pdf # アウトカム:甲羅の前面突起部の幅(FL, mm) # 要因1:種類(sp, B/O)Bは青色,Oはオレンジ色 # 要因2:性別(sex, F/M) twowayplot(FL, sp, sex, crabs) Anova(lm(FL ~ sp*sex, data=crabs), type=3) # spの主効果と交互作用効果が有意な例 Anova(lm(FL ~ sp*sex, data=crabs), type=2) # 各群50個体ずつなのでタイプは2でも3でも同じ anova(lm(FL ~ sp*sex, data=crabs)) # baseのanova()関数でも同じ # MASSパッケージのPima.tr # アリゾナ州フェニックス近くに居住する21歳以上のピマインディアンの女性の糖尿病関連データ # ソース:http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2245318/pdf/procascamc00018-0276.pdf # データ元:http://archive.ics.uci.edu/ml/datasets/Pima+Indians+Diabetes # アウトカム:Body Mass Index(bmi, kg/m^2) # 要因1:WHO基準による糖尿病判定(type, No/Yes) # 要因2:これまでの妊娠回数(npreg,回)を無産(NP),一子(PP),多子(MP)にカテゴリ化(parrous, NP/PP/MP) Pima <- Pima.tr parrous <- ifelse((Pima$npreg==0), 1, ifelse((Pima$npreg<2), 2, 3)) Pima$parrous <- as.factor(parrous) levels(Pima$parrous) <- c("NP","PP","MP") twowayplot(bmi, type, parrous, Pima) Anova(lm(bmi ~ type*parrous, data=Pima.tr), type=3) # 主効果も交互作用効果も有意でない例(1) # MASSパッケージのsurvey # アデレード大学学生の質問紙調査データ # ソース:http://www.planta.cn/forum/files_planta/modern_applied_statistics_with_s_192.pdf # データ元?:http://www.stats.ox.ac.uk/pub/MASS4/ # アウトカム:1分間の心拍数(Pulse) # 要因1:性別(Sex, ) # 要因2:運動頻度(Exer, ) twowayplot(Pulse, Sex, Exer, survey) Anova(lm(Pulse ~ Sex*Exer, data=survey), type=3) # 主効果も交互作用効果も有意でない例(2) # グラフ出力終了 dev.off()