課題内容

ある農学者は、植物の成長に対する肥料と日照時間の影響を調査しています。彼は異なる種類の肥料（肥料A、肥料B、肥料C）と異なる日照時間（8時間、10時間、12時間）を使って実験を行いました。それぞれの組み合わせで植物を育て、その成長（高さ）を測定しました。以下は実験から得られた植物の成長データです（単位：cm）。

肥料/日照時間 8時間 10時間 12時間肥料A 12, 14, 13 15, 16, 17 18, 19, 20 肥料B 14, 15, 16 17, 18, 19 20, 21, 22 肥料C 16, 17, 18 19, 20, 21 22, 23, 24

データの整理：5点

上記のデータを使って、二元配置の分散分析を行うためのデータフレームを作成してください。各観測値に対応する肥料の種類と日照時間を列に追加してください。

#データの作成1(日照時間をファクターとして定義した場合)
肥料 <- rep(c('A', 'B', 'C'), each = 9)
日照時間 <- factor(rep(rep(c(8, 10, 12), each = 3), 3))
成長 <- c(12, 14, 13, 15, 16, 17, 18, 19, 20, 
          14, 15, 16, 17, 18, 19, 20, 21, 22, 
          16, 17, 18, 19, 20, 21, 22, 23, 24)

df1 <- data.frame(肥料, 日照時間, 成長)
print(df1)

##    肥料 日照時間 成長
## 1     A        8   12
## 2     A        8   14
## 3     A        8   13
## 4     A       10   15
## 5     A       10   16
## 6     A       10   17
## 7     A       12   18
## 8     A       12   19
## 9     A       12   20
## 10    B        8   14
## 11    B        8   15
## 12    B        8   16
## 13    B       10   17
## 14    B       10   18
## 15    B       10   19
## 16    B       12   20
## 17    B       12   21
## 18    B       12   22
## 19    C        8   16
## 20    C        8   17
## 21    C        8   18
## 22    C       10   19
## 23    C       10   20
## 24    C       10   21
## 25    C       12   22
## 26    C       12   23
## 27    C       12   24

# データの作成2(日照時間を数値として定義した場合)
肥料 <- rep(c('A', 'B', 'C'), each=3*3)
日照時間 <- rep(c(8, 10, 12), each=3, times=3)
成長 <- c(12, 14, 13, 15, 16, 17, 18, 19, 20, 
          14, 15, 16, 17, 18, 19, 20, 21, 22, 
          16, 17, 18, 19, 20, 21, 22, 23, 24)

df2 <- data.frame(肥料, 日照時間, 成長)
print(df2)

##    肥料 日照時間 成長
## 1     A        8   12
## 2     A        8   14
## 3     A        8   13
## 4     A       10   15
## 5     A       10   16
## 6     A       10   17
## 7     A       12   18
## 8     A       12   19
## 9     A       12   20
## 10    B        8   14
## 11    B        8   15
## 12    B        8   16
## 13    B       10   17
## 14    B       10   18
## 15    B       10   19
## 16    B       12   20
## 17    B       12   21
## 18    B       12   22
## 19    C        8   16
## 20    C        8   17
## 21    C        8   18
## 22    C       10   19
## 23    C       10   20
## 24    C       10   21
## 25    C       12   22
## 26    C       12   23
## 27    C       12   24

ファクターと数値の違い: 日照時間がファクターとして扱われるか、数値として扱われるかの違い。モデルの解釈: ファクターとして扱うと各レベル間の比較が行われるが、数値として扱うと連続変数としての影響が評価される。

仮説の設定：

主効果：肥料の種類（肥料A、肥料B、肥料C）による植物の成長の違いを検証します。主効果：日照時間（8時間、10時間、12時間）による植物の成長の違いを検証します。

交互作用効果：肥料の種類と日照時間の交互作用が植物の成長に及ぼす影響を検証します。

二元配置の分散分析の実施：10点

二元配置の分散分析を実行し、各効果の有意性を検定してください。分散分析表を作成し、F値とp値を報告してください。

# CRANミラーの設定
options(repos = c(CRAN = "https://cloud.r-project.org"))
# 必要なライブラリの読み込み
install.packages("car")

## 
## The downloaded binary packages are in
##  /var/folders/nm/m79f0x2d79768ww18g9h71340000gp/T//RtmpjuqcRt/downloaded_packages

library(car)

## Loading required package: carData

# 二元配置の分散分析モデルの作成)(df1について)
model <- aov(成長 ~ 肥料 * 日照時間, data=df1)

# 分散分析表の作成
anova_table <- Anova(model, type=2) 
print(anova_table)

## Anova Table (Type II tests)
## 
## Response: 成長
##               Sum Sq Df F value   Pr(>F)    
## 肥料              72  2      36 5.12e-07 ***
## 日照時間         162  2      81 1.00e-09 ***
## 肥料:日照時間      0  4       0        1    
## Residuals         18 18                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#二元配置の分散分析モデルの作成(df2について)
model2<-aov(成長 ~ 肥料 * 日照時間, data=df2)
summary(model2)

##               Df Sum Sq Mean Sq F value   Pr(>F)    
## 肥料           2     72   36.00      42 4.58e-08 ***
## 日照時間       1    162  162.00     189 5.71e-12 ***
## 肥料:日照時間  2      0    0.00       0        1    
## Residuals     21     18    0.86                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 効果量の計算
# 必要なパッケージのインストール
install.packages("effectsize")

## 
##   There is a binary version available but the source version is later:
##            binary source needs_compilation
## effectsize  0.8.8  0.8.9             FALSE

## installing the source package 'effectsize'

# パッケージの読み込み
library(car)
library(effectsize)

# 効果量の計算 
effect_sizes <- eta_squared(model)
print(effect_sizes)

## # Effect Size for ANOVA (Type I)
## 
## Parameter     | Eta2 (partial) |       95% CI
## ---------------------------------------------
## 肥料          |           0.80 | [0.63, 1.00]
## 日照時間      |           0.90 | [0.81, 1.00]
## 肥料:日照時間 |       6.66e-31 | [0.00, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Pythonで行う場合

anova_results = anova_lm(model, typ=3) のようにypt=3として、タイプⅢ平方和を使用する。

平方和の違い

Type I 平方和（Sequential Sum of Squares）

　概要: 　　　　Type I 平方和は、変数をモデルに追加する順序に依存します。　　　最初に追加された変数の平方和を計算し、その後に追加される変数の平方和を計算します。　　　モデルの順序が重要であり、各変数が追加された時点での残差に基づいて計算されます。

　特徴: 　　　順序に依存。　　　モデルに追加される順序によって結果が変わる可能性がある。

Type II 平方和（Hierarchical or Partially Sequential Sum of Squares）　　概要: 　　　Type II 平方和は、他の要因をすべて含んだ後に各要因の主効果を評価しますが、交互作用の影響は無視します。　　　各要因の効果は、他の要因をすべて含んだ後に計算されます。

　特徴: 　　　順序に依存しない。　　　　交互作用の影響を無視。

Type III 平方和（Marginal Sum of Squares）　　概要: 　　　Type III 平方和は、すべての要因と交互作用を含んだ状態で各要因の効果を評価します。　　　　各要因の効果は、他の要因および交互作用の影響をすべて含んだ後に計算されます。

　特徴: 　　　　順序に依存しない。　　　　交互作用を含む完全なモデルに基づく。

※例題の場合、本来であれば交互作用を考える必要があるので、TypeⅢ平方和を選択する必要があるが、Rで計算する場合、エラーが出てしまう。このため、対象を設定するのが望ましい。

# データの作成
肥料 <- rep(c('A', 'B', 'C'), each = 9)
日照時間 <- factor(rep(rep(c(8, 10, 12), each = 3), 3))
成長 <- c(12, 14, 13, 15, 16, 17, 18, 19, 20, 
          14, 15, 16, 17, 18, 19, 20, 21, 22, 
          16, 17, 18, 19, 20, 21, 22, 23, 24)

df3 <- data.frame(肥料, 日照時間, 成長)

# モデルの作成
model3 <- aov(成長 ~ 肥料 * 日照時間, data = df3)

# 対照の設定
options(contrasts = c("contr.sum", "contr.poly"))

# Type III 平方和の計算
anova_results_III <- Anova(model3, type="III")
print(anova_results_III)

## Anova Table (Type III tests)
## 
## Response: 成長
##               Sum Sq Df F value    Pr(>F)    
## (Intercept)      507  1     507 1.233e-14 ***
## 肥料              24  2      12 0.0004878 ***
## 日照時間          54  2      27 3.815e-06 ***
## 肥料:日照時間      0  4       0 1.0000000    
## Residuals         18 18                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

結果としては、タイプⅡ平方和を使った場合と大差ない結果が得られる。