R version 3.2.4 (2016-03-10) -- "Very Secure Dishes" Copyright (C) 2016 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 (64-bit) R は、自由なソフトウェアであり、「完全に無保証」です。 一定の条件に従えば、自由にこれを再配布することができます。 配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 R は多くの貢献者による共同プロジェクトです。 詳しくは 'contributors()' と入力してください。 また、R や R のパッケージを出版物で引用する際の形式については 'citation()' と入力してください。 'demo()' と入力すればデモをみることができます。 'help()' とすればオンラインヘルプが出ます。 'help.start()' で HTML ブラウザによるヘルプがみられます。 'q()' と入力すれば R を終了します。 > # 組み込みデータを使って二元配置分散分析 > # 二元配置の平均と標準誤差または標準偏差の折れ線グラフを描く関数 > 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） > # 要因１：キャベツの品種（Cult, c39/c52） > # 要因２：植え付けた日（Date, d16/d20/d21） > twowayplot(HeadWt, Cult, Date, cabbages) $MEAN : c39 : d16 [1] 3.18 ------------------------------------------------------------ : c52 : d16 [1] 2.26 ------------------------------------------------------------ : c39 : d20 [1] 2.8 ------------------------------------------------------------ : c52 : d20 [1] 3.11 ------------------------------------------------------------ : c39 : d21 [1] 2.74 ------------------------------------------------------------ : c52 : d21 [1] 1.47 $SD : c39 : d16 [1] 0.9566144 ------------------------------------------------------------ : c52 : d16 [1] 0.4452215 ------------------------------------------------------------ : c39 : d20 [1] 0.2788867 ------------------------------------------------------------ : c52 : d20 [1] 0.7908505 ------------------------------------------------------------ : c39 : d21 [1] 0.9834181 ------------------------------------------------------------ : c52 : d21 [1] 0.2110819 $SE : c39 : d16 [1] 0.302508 ------------------------------------------------------------ : c52 : d16 [1] 0.1407914 ------------------------------------------------------------ : c39 : d20 [1] 0.08819171 ------------------------------------------------------------ : c52 : d20 [1] 0.2500889 ------------------------------------------------------------ : c39 : d21 [1] 0.3109841 ------------------------------------------------------------ : c52 : d21 [1] 0.06674995 $N : c39 : d16 [1] 10 ------------------------------------------------------------ : c52 : d16 [1] 10 ------------------------------------------------------------ : c39 : d20 [1] 10 ------------------------------------------------------------ : c52 : d20 [1] 10 ------------------------------------------------------------ : c39 : d21 [1] 10 ------------------------------------------------------------ : c52 : d21 [1] 10 > Anova(lm(HeadWt ~ Cult*Date, data=cabbages), type=3) # 主効果，交互作用効果とも有意な例 Anova Table (Type III tests) Response: HeadWt Sum Sq Df F value Pr(>F) (Intercept) 403.52 1 856.0629 < 2.2e-16 *** Cult 5.89 1 12.4969 0.0008451 *** Date 7.71 2 8.1744 0.0007920 *** Cult:Date 6.89 2 7.3046 0.0015571 ** Residuals 25.45 54 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(lm(HeadWt ~ Cult*Date, data=cabbages)) # サイズ不斉でタイプ3でないので結果が異なり，これは不適切 Analysis of Variance Table Response: HeadWt Df Sum Sq Mean Sq F value Pr(>F) Cult 1 5.8907 5.8907 12.4969 0.0008451 *** Date 2 7.7063 3.8532 8.1744 0.0007920 *** Cult:Date 2 6.8863 3.4432 7.3046 0.0015571 ** Residuals 54 25.4540 0.4714 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > # MASSパッケージのcrabs > # 西オーストラリアの蟹の形態 > # ソース：http://duch.mimuw.edu.pl/~pokar/StatystykaII/DANE/CampbellMahon.pdf > # アウトカム：甲羅の前面突起部の幅（FL, mm） > # 要因１：種類（sp, B/O）Bは青色，Oはオレンジ色 > # 要因２：性別（sex, F/M) > twowayplot(FL, sp, sex, crabs) $MEAN : B : F [1] 13.27 ------------------------------------------------------------ : O : F [1] 17.594 ------------------------------------------------------------ : B : M [1] 14.842 ------------------------------------------------------------ : O : M [1] 16.626 $SD : B : F [1] 2.627814 ------------------------------------------------------------ : O : F [1] 2.973997 ------------------------------------------------------------ : B : M [1] 3.202492 ------------------------------------------------------------ : O : M [1] 3.514971 $SE : B : F [1] 0.3716291 ------------------------------------------------------------ : O : F [1] 0.4205867 ------------------------------------------------------------ : B : M [1] 0.4529008 ------------------------------------------------------------ : O : M [1] 0.497092 $N : B : F [1] 50 ------------------------------------------------------------ : O : F [1] 50 ------------------------------------------------------------ : B : M [1] 50 ------------------------------------------------------------ : O : M [1] 50 > Anova(lm(FL ~ sp*sex, data=crabs), type=3) # spの主効果と交互作用効果が有意な例 Anova Table (Type III tests) Response: FL Sum Sq Df F value Pr(>F) (Intercept) 48566 1 5064.0933 < 2.2e-16 *** sp 466 1 48.6270 4.627e-11 *** sex 5 1 0.4755 0.49128 sp:sex 81 1 8.4091 0.00416 ** Residuals 1880 196 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > Anova(lm(FL ~ sp*sex, data=crabs), type=2) # 各群50個体ずつなのでタイプは2でも3でも同じ Anova Table (Type II tests) Response: FL Sum Sq Df F value Pr(>F) sp 466.35 1 48.6270 4.627e-11 *** sex 4.56 1 0.4755 0.49128 sp:sex 80.64 1 8.4091 0.00416 ** Residuals 1879.69 196 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(lm(FL ~ sp*sex, data=crabs)) # baseのanova()関数でも同じ Analysis of Variance Table Response: FL Df Sum Sq Mean Sq F value Pr(>F) sp 1 466.35 466.35 48.6270 4.627e-11 *** sex 1 4.56 4.56 0.4755 0.49128 sp:sex 1 80.64 80.64 8.4091 0.00416 ** Residuals 196 1879.69 9.59 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > # 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） > # 要因１：WHO基準による糖尿病判定（type, No/Yes） > # 要因２：これまでの妊娠回数（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) $MEAN : No : NP [1] 33.1375 ------------------------------------------------------------ : Yes : NP [1] 36.68333 ------------------------------------------------------------ : No : PP [1] 29.84722 ------------------------------------------------------------ : Yes : PP [1] 34.96667 ------------------------------------------------------------ : No : MP [1] 31.21375 ------------------------------------------------------------ : Yes : MP [1] 34.15532 $SD : No : NP [1] 8.245716 ------------------------------------------------------------ : Yes : NP [1] 4.819814 ------------------------------------------------------------ : No : PP [1] 6.981956 ------------------------------------------------------------ : Yes : PP [1] 7.175653 ------------------------------------------------------------ : No : MP [1] 5.608306 ------------------------------------------------------------ : Yes : MP [1] 4.225945 $SE : No : NP [1] 2.061429 ------------------------------------------------------------ : Yes : NP [1] 1.39136 ------------------------------------------------------------ : No : PP [1] 1.163659 ------------------------------------------------------------ : Yes : PP [1] 2.391884 ------------------------------------------------------------ : No : MP [1] 0.6270277 ------------------------------------------------------------ : Yes : MP [1] 0.6164175 $N : No : NP [1] 16 ------------------------------------------------------------ : Yes : NP [1] 12 ------------------------------------------------------------ : No : PP [1] 36 ------------------------------------------------------------ : Yes : PP [1] 9 ------------------------------------------------------------ : No : MP [1] 80 ------------------------------------------------------------ : Yes : MP [1] 47 > Anova(lm(bmi ~ type*parrous, data=Pima.tr), type=3) # 主効果も交互作用効果も有意でない例(1) Anova Table (Type III tests) Response: bmi Sum Sq Df F value Pr(>F) (Intercept) 17870.2 1 513.7363 < 2e-16 *** type 125.9 1 3.6207 0.05853 . parrous 60.3 1 1.7339 0.18945 type:parrous 19.3 1 0.5551 0.45715 Residuals 6817.8 196 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 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） > # 要因１：性別（Sex, ） > # 要因２：運動頻度（Exer, ） > twowayplot(Pulse, Sex, Exer, survey) $MEAN : Female : Freq [1] 73.60976 ------------------------------------------------------------ : Male : Freq [1] 70.67925 ------------------------------------------------------------ : Female : None [1] 71.42857 ------------------------------------------------------------ : Male : None [1] 80.5 ------------------------------------------------------------ : Female : Some [1] 77 ------------------------------------------------------------ : Male : Some [1] 75.0303 $SD : Female : Freq [1] 12.49175 ------------------------------------------------------------ : Male : Freq [1] 9.593214 ------------------------------------------------------------ : Female : None [1] 11.41428 ------------------------------------------------------------ : Male : None [1] 15.20417 ------------------------------------------------------------ : Female : Some [1] 10.27026 ------------------------------------------------------------ : Male : Some [1] 13.50112 $SE : Female : Freq [1] 1.784536 ------------------------------------------------------------ : Male : Freq [1] 1.189892 ------------------------------------------------------------ : Female : None [1] 3.441534 ------------------------------------------------------------ : Male : None [1] 4.216877 ------------------------------------------------------------ : Female : Some [1] 1.348551 ------------------------------------------------------------ : Male : Some [1] 2.134715 $N : Female : Freq [1] 49 ------------------------------------------------------------ : Male : Freq [1] 65 ------------------------------------------------------------ : Female : None [1] 11 ------------------------------------------------------------ : Male : None [1] 13 ------------------------------------------------------------ : Female : Some [1] 58 ------------------------------------------------------------ : Male : Some [1] 40 > Anova(lm(Pulse ~ Sex*Exer, data=survey), type=3) # 主効果も交互作用効果も有意でない例(2) Anova Table (Type III tests) Response: Pulse Sum Sq Df F value Pr(>F) (Intercept) 594993 1 4479.1870 < 2e-16 *** Sex 52 1 0.3879 0.53419 Exer 695 2 2.6170 0.07572 . Sex:Exer 512 2 1.9261 0.14862 Residuals 24574 185 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > # グラフ出力終了 > dev.off() null device 1 >