Прогнозы NLME при разных уровнях воздействия в R - PullRequest
Новая
       37

Прогнозы NLME при разных уровнях воздействия в R

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

У меня есть простой пример эксперимента на полевом участке с 3 точками, 4 блоками в каждой точке и 8 нормами внесения удобрений на этих участках.ответ был измерен на этих участках.

Это данные: 

#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, 1071,    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")

Я настроил функцию экспоненциального плато с местоположениями и блоками в качестве случайных эффектов.

#exponential plateau function
expfunc <- function(beta0, beta1, beta2, x){
             y <- beta0 * (1 - exp(-exp(beta1) / 1000 * x + beta2))
             return(y)}

# model fit with blocks and locations as random effects
exp.loc_block <- nlme(response ~ expfunc(beta0, beta1, beta2, treat),
                   data = resp.data,
                   fixed = list(beta0 ~ 1, beta1 ~ 1, beta2 ~ 1),
                   random = list(loc = pdDiag(beta0 + beta1 + beta2 ~ 1),
                                 block = pdDiag(beta0 + beta1 + beta2 ~ 1)),
                   start = c(12000, 3, -1),
                   na.action = na.omit,
                   verbose = F)

summary(exp.loc_block)

Теперь я хочу получить прогноз ответа на уровне 0 (фиксированный эффект).

#prediction at population level (fixed effect)
xvals <- with(resp.data,seq(min(treat),max(treat),length.out=100))
pframe <- data.frame(treat=xvals)
pframe$pred.resp.fix <- predict(exp.loc_block,newdata=pframe,level=0)

Пока все работает нормально.По ним я получаю прогноз на уровне блока и локации.

#prediction at block nested in location level (random effect)
Loc.names.vector <- unique(resp.data$loc)
block.names.vector <- unique(resp.data$block)
pframe2 <- with(resp.data,data.frame(treat=rep(xvals, 12)))
pframe2$loc <- as.factor(rep(Loc.names.vector, each = 400))
pframe2$block <- as.factor(rep(rep(block.names.vector, each = 100), 3))
pframe2$yield.exp <- predict(exp.loc_block,newdata=pframe2)

Итак, мой вопрос, есть ли способ предсказать ответ на уровне местоположения?

Действительно ценю любой намек.

Спасибо.

...