1

私は次の論文で 3 因子のネストされた ANOVA 分析を再現しようとしています: Underwood, AJ (1993) The Mechanics of spatially replicated sampling programs to detect environment impacts in a variable world.

この例のデータ (表 3、Underwood 1993 から) は、次の方法で作成できます。

dat <-
structure(list(B = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
), .Label = c("A", "B"), class = "factor"), C = structure(c(2L,
2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("C", "I"), class = "factor"),
    Times = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L,
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
    4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor"),
    Locations = c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L,
    1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L,
    3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L,
    2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L,
    1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L,
    3L), X = c(59L, 51L, 45L, 46L, 40L, 32L, 39L, 32L, 25L, 51L,
    44L, 37L, 55L, 47L, 41L, 31L, 38L, 45L, 41L, 47L, 55L, 43L,
    36L, 29L, 23L, 30L, 37L, 57L, 50L, 43L, 36L, 44L, 51L, 39L,
    29L, 23L, 38L, 44L, 52L, 31L, 38L, 45L, 42L, 35L, 28L, 52L,
    44L, 37L, 51L, 43L, 37L, 38L, 31L, 24L, 60L, 52L, 46L, 30L,
    37L, 44L, 41L, 34L, 27L, 53L, 46L, 39L, 40L, 34L, 26L, 21L,
    27L, 35L), Times.unique = structure(c(5L, 5L, 5L, 5L, 5L,
    5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L,
    7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
    8L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L,
    4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("A_1", "A_2", "A_3",
    "A_4", "B_1", "B_2", "B_3", "B_4"), class = "factor")), .Names = c("B",
"C", "Times", "Locations", "Y", "Times.unique"), row.names = c(NA,
-72L), class = "data.frame")

dat

データ フレーム データには 4 つの要素があります。

B - 「A」と「B」の 2 つのレベルがあります (前 v 後)

時間 - 8 つのレベル、「B」の前に 4 つ、「A」の後に 4 つ、それぞれに 1:4 としてコード化されています。変数 Times.unique は同じものですが、各時間 (前後) に固有のコードがあることに注意してください。

場所 - 3 つのレベルがあり、前後に毎回測定されます。

C - 2 つのレベル コントロール (C) と (I) があります。注: 2 か所はコントロールで、1 か所はインパクトです

混合モデル (lmer) を使用してそのような設計を分析する方法については明確ですが、いくつかのシミュレーションを実行して彼の方法を比較できるように、彼の例を正確に再現したいと思います。

特に、表 4 の列「a」に示されている SS 値を再現しようとしています。彼は、次の項について SS 値と df 値をもつ計画を当てはめます。

B -> SS = 66.13、df = 1

回 (B) -> SS = 280.64、自由度 = 6

ロケーション -> SS = 283.86、df = 2

B x ロケーション -> SS = 29.26、df = 2

時間(B) x ロケーション -> SS = 575.45、df = 12

残差 -> SS = 2420.00、自由度 = 48

トータル -> SS = 6208.34、df = 71

Times(B) という用語は、Before/After 処理「B」内にネストされた Times を表していると思います。この例では、場所がコントロールと影響の処理によるものであることを無視し、因子 C を完全に除外します。

このネストされた anova を再現するために考えられるすべての組み合わせを試しました。一意の Times コーディングと、B 内で 1:4 としてコーディングされた Times (前と後) の両方を使用します。%in%、/、Error() 引数、および Anova from car を使用して、計算される SS のタイプを変更しようとしました。%in% および / ネストされた適合の例は次のとおりです。

aov(Y~B+Locations+Times%in%B+B:Locations+Times%in%B:Locations, data=dat)
aov(Y~B+Locations+B/Times+B:Locations+B/Times:Locations, data=dat)

特に 2 つの交互作用項については、アンダーウッドの SS 値を正確に再現できないようです。友人が、SS 値を正確に再現できる statistix にモデルを適合させてくれたので、このモデルの上記の SS 値を取得することが可能です。

このモデルをRに適合させるのを手伝ってくれる人はいますか? それをより大きなシミュレーションに組み込みたいのですが、Underwood 1993 SS の値が正確に再現されるように、モデルを R で実行できるようにする必要がありますか?

4

1 に答える 1

1

あなたの問題は、dat$Locationsそれが要因(3つの一意の場所)である必要がある場合、整数であることです。ヒントの 1 つは、あなたの ANOVA ラインは Locations が 1 df しか占めていないと考えているのに対し、Underwood は 2 を与えているということです。

次の行を追加するだけです。

dat$Locations = factor(dat$Locations)

そして、あなたのコード行はアンダーウッドの結果を完全に再現します:

aov(Y~B+Locations+B/Times+B:Locations+B/Times:Locations, data=dat)
#Call:
#   aov(formula = Y ~ B + Locations + B/Times + B:Locations + B/Times:Locations, 
#    data = dat)
#
#Terms:
#                        B Locations   B:Times B:Locations B:Locations:Times
#Sum of Squares    66.1250 2836.8611  280.6389     29.2500          575.4444
#Deg. of Freedom         1         2         6           2                12
#                Residuals
#Sum of Squares  2420.0000
#Deg. of Freedom        48
于 2013-10-02T05:21:19.793 に答える