стандартное отклонение для коэффициентов случайных эффектов - PullRequest
Новая
       58

стандартное отклонение для коэффициентов случайных эффектов

0 голосов
/ 09 октября 2018

У меня есть вопрос о том, как получить значение неопределенности для коэффициентов случайного эффекта в nlme.

Я составил набор данных для отклика урожайности зерна на внесение удобрений в разных местах (агрономические эксперименты).Модель представляет собой модель квадратичного плато.Параметры модели: максимальная урожайность (Ymax), удобрение для максимальной урожайности (M) и начальный уклон (k).Для данного местоположения M составляет 202,4, а Ymax составляет 14306,79.Я хочу сообщить об этих коэффициентах с некоторой степенью неопределенности.Как я могу получить значение SD или SE для этих коэффициентов?

Вот воспроизводимый пример. 

library(nlme)
#dataframe
loc <- c("Loc1", "Loc2", "Loc3")
block <- c("Block_1", "Block_2", "Block_3", "Block_4")
treat <- as.numeric(c("0","40","80","120","160","200","240","280"))
empty <-expand.grid(treat,  block, loc)
response <- as.numeric(c(7064,  9250,   12306,  13549,  13300,  13973,    
14749,  14086, 7680, 11426, 12874,  12556,  14274,  14289,  15295,  14587,
                     8445,  11588,  13223,  13322,  13508,  13616,  13747,  13352,
                     9454,  11104,  12462,  13373,  14060,  14576,  14133,  14427,
                     5463,  8689,   10194,  11996,  13475,  12544,  12856,  11557,
                     5251,  7537,   12074,  12438,  12120,  11312,  9908,   12841,
                     4643,  7513,   10499,  12423,  12177,  12545,  12876,  13047,
                     4992,  9458,   11071,  12104,  13552,  12602,  13210,  14428,
                     4061,  3959,   5871,   8016,   9472,   11554,  12525,  12636,
                     4598,  7717,   7274,   8476,   9433,   10768,  10275,  8200,
                     4862,  5727,   6468,   8532,   10662,  12054,  12227,  12672,
                     5218,  7878,   8238,   10303,  10331,  13337,  12877,  11661))

resp.data <- cbind(empty, response)
resp.data <- resp.data[c("Var3", "Var2", "Var1", "response")]
names(resp.data) <- c("loc", "block", "treat", "response")

f.quad.loc <- nlme(response ~ ymax*(treat > M) +
                  (ymax - (k/2) * (M-treat)^2) * (treat < M),
               data = resp.data,
               fixed = list(ymax ~ 1,
                            M ~ 1,
                            k ~ 1),
               random = list(loc = pdDiag(ymax + M + k ~ 1)),
               start = c(12000,
                         100,
                         0.2),
               na.action = na.omit,
               verbose = F)


r.effects <- random.effects(f.quad.loc)
f.effects <- fixed.effects(f.quad.loc)
r.effects$ymax <- r.effects$ymax + f.effects[1]
r.effects$M <- r.effects$M + f.effects[2]
r.effects$k <- r.effects$k+ f.effects[3]

r.effects

Заранее большое спасибо.

...