Как учесть повторные измерения и корреляцию между факторами в древовидной форме ANOVA? - PullRequest
Новая
       27

Как учесть повторные измерения и корреляцию между факторами в древовидной форме ANOVA?

0 голосов
/ 25 июня 2019

Я новичок в статистике и не понимаю различий между различными функциями, используемыми для линейных моделей. Мне нужна помощь, чтобы выбрать, какую функцию использовать и как ее написать, чтобы учесть все параметры.

Мой экспериментальный дизайн состоит из трех факторов: Фактор 1. происхождение генотипа (3 уровня); Фактор 2. генотип дерева (10 уровней); Фактор 3. Тип субстрата (3 уровня).

Дополнительно я повторил меры: Я измерял рост деревьев 3 раза в течение сезона (4 дня, 57 дней и 91 день после начала эксперимента).

Я пытался использовать gls (nlme):

Сначала я попробовал просто включить каждый фактор 

mod1 <- gls(data = my_data, model = Height ~ Origin*Genotype*Substrate*Date)

Error in glsEstimate(glsSt, control = glsEstControl) : 
  l'ajustement "gls" calculé est singulier, de rang 61

Но похоже, что он не работает из-за "Происхождения", поэтому я удалил его. 

mod1 <- gls(data = my_data, model = Height ~ Genotype*Substrate*Date)

Я попытался указать корреляцию между повторными измерениями. 

mod2 <- gls(data = my_data, model = Height ~ Genotype*Substrate*Date,
            correlation = corCAR1(form = ~ Date | Sample_ID))

Это работает, но я не знаю, как включить корреляцию между генотипом дерева и его происхождением.

Кроме того, добавление корреляции не улучшает модель: 

> AICc(mod1, mod2)
     df     AICc
mod1 61 2385.981
mod2 62 2389.311

Я делаю это неправильно?

Я читал, вы также можете использовать lmer (lme4) для такого анализа, было бы лучше в моем случае?

Вот мои данные 

my_data <- structure(list(Sample_ID = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 
20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 
33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 
46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 
59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 
72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 
85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 1L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 
20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 
33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 
46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 
59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 
72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 
85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L, 1L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 
20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 
33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 
46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 
59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 
72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 
85L, 86L, 87L, 88L, 89L, 90L, 91L, 92L), .Label = c("08C1", "08C2", 
"08C3", "08M1", "08M2", "08M3", "08W1", "08W2", "08W3", "08W4", 
"09C1", "09C2", "09C3", "09M1", "09M2", "09M3", "09W1", "09W2", 
"09W3", "09W4", "10C1", "10C2", "10C3", "10M1", "10M2", "10M3", 
"10W1", "10W2", "10W3", "13C1", "13C2", "13C3", "13M1", "13M2", 
"13M3", "13W1", "13W2", "13W3", "16C1", "16C2", "16M1", "16M2", 
"16W1", "16W2", "16W3", "21C1", "21C2", "21C3", "21M1", "21M2", 
"21M3", "21W1", "21W2", "21W3", "23C1", "23C2", "23C3", "23M1", 
"23M2", "23M3", "23M4", "23W1", "23W2", "23W3", "25C1", "25C2", 
"25C3", "25M1", "25M2", "25M3", "25W1", "25W2", "25W3", "29C1", 
"29C2", "29C3", "29M1", "29M2", "29M3", "29W1", "29W2", "29W3", 
"33C1", "33C2", "33C3", "33M1", "33M2", "33M3", "33W1", "33W2", 
"33W3", "33W4"), class = "factor"), Origin = structure(c(3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 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, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 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, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 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, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L), .Label = c("Moly", "Road", "West"), class = "factor"), 
    Genotype = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 2L, 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, 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, 7L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 
    10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 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, 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, 
    7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 
    9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
    10L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 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, 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, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 
    10L, 10L, 10L, 10L, 10L), .Label = c("PB08", "PB09", "PB10", 
    "PB13", "PB16", "PB21", "PB23", "PB25", "PB29", "PB33"), class = "factor"), 
    Substrate = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 
    3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 
    2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 
    1L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 
    1L, 1L, 1L, 2L, 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, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 
    3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 
    1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 
    3L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 
    3L, 3L, 1L, 1L, 1L, 2L, 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, 3L, 1L, 1L, 1L, 2L, 2L, 
    2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
    1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 
    3L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 
    2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 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, 3L), .Label = c("Control", 
    "Molybdenum", "Westwood"), class = "factor"), Date = structure(c(17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 17280, 
    17280, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 17333, 
    17333, 17333, 17333, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 17364, 
    17364, 17364, 17364, 17364, 17364), class = "Date"), Height = c(129.5, 
    140, 122, 112, 109, 127, 109, 134.5, 96.5, 109, 114.5, 104, 
    107, 117, 79, 119.5, 101.5, 63.5, 122, 74, 109.5, 112, 91.5, 
    104, 109, 104, 117, 114.5, 117, 114, 124.5, 111.5, 122, 99, 
    134.5, 84, 142, 117, 112, 71, 74, 76, 117, 48.5, 114.5, 81.5, 
    89, 109, 114.5, 78.5, 106.5, 89, 119.5, 119.5, 96.5, 96.5, 
    99, 91.5, 106.5, 33, 101.5, 89, 112, 117, 104, 122, 78.5, 
    89, 66, 76, 106.5, 104, 104, 94, 127, 132, 122, 122, 109, 
    129.5, 139.5, 78.5, 104, 139.5, 122, 122, 81.5, 134.5, 91.5, 
    139.5, 122, 114.5, 162.5, 193, 167.5, 137, 160, 172.5, 117, 
    167.5, 114.5, 172.5, 127, 109, 134.5, 144.5, 101.5, 137, 
    117, 111.5, 147.5, 134.5, 139.5, 165, 129.5, 152.5, 155, 
    137, 147.5, 142, 137, 147.5, 157.5, 160, 124.5, 142, 160, 
    144.5, 157.5, 127, 180.5, 147.5, 147.5, 142, 172.5, 122, 
    157.5, 157.5, 137, 172.5, 160, 111.5, 157.5, 144.5, 175.5, 
    157.5, 134.5, 150, 134.5, 137, 157.5, 111.5, 144.5, 122, 
    139.5, 132, 147.5, 167.5, 152.5, 142, 134.5, 152.5, 160, 
    160, 142, 178, 134.5, 185.5, 178, 160, 160, 160, 162.5, 139.5, 
    160, 185.5, 180.5, 170, 142, 183, 129.5, 165, 172.5, 172.5, 
    175, 193, 170, 137, 160, 172.5, 122, 167.5, 114.5, 172.5, 
    127, 124.5, 134.5, 144.5, 101.5, 139.5, 122, 114, 147.5, 
    134.5, 139.5, 165, 132, 155, 155, 142, 147.5, 142, 137, 147.5, 
    157.5, 160, 124.5, 142, 160, 144.5, 160, 127, 180.5, 150, 
    147.5, 142, 172.5, 124.5, 157.5, 157.5, 137, 175, 160, 114, 
    157.5, 144.5, 175.5, 157.5, 134.5, 152, 134.5, 144.5, 157.5, 
    119, 147, 124.5, 155, 132, 147.5, 167.5, 152.5, 147, 137, 
    155, 160, 160, 142, 178, 134.5, 185.5, 178, 160, 160, 160, 
    162.5, 139.5, 162.5, 188, 183, 172.5, 144.5, 183, 129.5, 
    165, 172.5, 172.5)), row.names = c(NA, -276L), class = "data.frame")

Итак, вот вопросы, с которыми мне нужна помощь (в таком порядке):

  1. Какая модель подойдет для моего случая?

  2. Есть ли эффект происхождения генотипа: как я могу включить в модель корреляцию между генотипом (фактор 2) и его происхождением (фактор 1)?

  3. Есть ли время: как я могу включить повторные измерения в модель?

...