如何在 NLME 中获得正确级别的预测而不是 NA？

Question

我有一个简单的实验（这里是一些虚构的数据），有 3 个位置（"loc"），位置内的块（"block"），一组 8 个处理（"treat), and a response variable ("response")。指数平台函数被调整到这个数据（响应作为治疗的函数）并且 NLME 函数被用来使用混合模型分析整个实验。指数平台函数的参数被认为是固定的部分而 loc 和 block 是随机的。

我的问题是模型的预测。我无法预测对 loc 级别的治疗的反应。我能够将它达到人口水平（固定效应），但是当尝试在 loc 水平上进行预测时，得到了所有 NA。我设置 "levels= " 选项的方式有问题吗？

这是我编造的数据

#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
library(nlme)
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)

# This provides evidence that levels loc and block within loc are having random effects calculated
ranef(exp.loc_block)

这是我做出预测的方式

#creation of the empty dataframe
#creation of the empty dataframe
xvals <- seq(min(resp.data$treat),max(resp.data$treat),length.out=100)
Loc.names.vector <- unique(resp.data$loc)
block <- c("Block_1","Block_2","Block_3","Block_4")
pframe.lp<-expand.grid(xvals, Loc.names.vector, block)
names(pframe.lp)[1]<-"treat"   
names(pframe.lp)[2]<-"loc"
names(pframe.lp)[3]<-"block"  
#prediction at location level (random effect)
pframe.lp$yield.exp <- predict(exp.loc_block,newdata=pframe.lp, level=0)
pframe.lp$yield.exp <- predict(exp.loc_block,newdata=pframe.lp, level=1)
pframe.lp$yield.exp <- predict(exp.loc_block,newdata=pframe.lp, level=2)

level=0 的预测适用于固定效应的预测。 level=2 的预测适用于 block/loc 随机级别的预测 level=1 的预测给了我 NAs.. 这将是 loc 级别的预测，这是我最感兴趣的级别。

是不是等级选项有误？

这里是使用nlme函数进行预测的解释 https://rdrr.io/cran/nlme/man/predict.nlme.html

Answer 1

没有看到任何 NA。你能展示代码和足够的输出来记录你的问题吗？（我改为使用 levels=0:2 并查看 str(pframe.lp) 然后从 summary(pframe.lp$yield.exp$predict.loc):

> pframe.lp$yield.exp <- predict(exp.loc_block,newdata=pframe.lp, level=0:2)
> 
> str(pframe.lp)
'data.frame':   1200 obs. of  4 variables:
 $ treat    : num  0 2.83 5.66 8.48 11.31 ...
 $ loc      : Factor w/ 3 levels "Loc1","Loc2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ block    : Factor w/ 4 levels "Block_1","Block_2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ yield.exp:'data.frame':  1200 obs. of  5 variables:
  ..$ loc          : Factor w/ 3 levels "Loc1","Loc2",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..$ block        : Factor w/ 12 levels "Loc1/Block_1",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..$ predict.fixed: num  5987 6188 6384 6575 6762 ...
  ..$ predict.loc  : num  7867 8125 8374 8613 8842 ...
  ..$ predict.block: num  7867 8121 8366 8601 8827 ...
 - attr(*, "out.attrs")=List of 2
  ..$ dim     : int  100 3 4
  ..$ dimnames:List of 3
  .. ..$ Var1: chr  "Var1=  0.000000" "Var1=  2.828283" "Var1=  5.656566" "Var1=  8.484848" ...
  .. ..$ Var2: chr  "Var2=Loc1" "Var2=Loc2" "Var2=Loc3"
  .. ..$ Var3: chr  "Var3=Block_1" "Var3=Block_2" "Var3=Block_3" "Var3=Block_4"
> summary(pframe.lp$yield.exp$predict.loc)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4380    9219   11426   10928   13103   14279

这种用法并不理想。将数据框添加到现有数据框是一种糟糕的 R 做法。更好的是：

pframe.lp$yield.exp <- data.matrix( predict( 
                                  exp.loc_block,newdata=pframe.lp, level=0:2) )

如何在 NLME 中获得正确级别的预测而不是 NA？

How to get predictions in NLME at the correct level instead of NAs?

r

nlme