【发布时间】：2019-08-31 14:46:20
【问题描述】：
我有一个简单的实验（这里是一些虚构的数据），其中包含 3 个位置（“loc”）、位置内的块（“block”）、一组 8 个处理（“treat”）和一个响应变量（ “响应”）。根据该数据调整指数平台函数（响应作为处理的函数），并使用 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 的预测适用于预测固定效应。级别 = 2 的预测适用于块/位置随机级别的预测 level=1 的预测给了我 NAs.. 这将是 loc 级别的预测，这是我最感兴趣的。
关卡选项有错误吗？
这里是使用 nlme 函数进行预测的解释 https://rdrr.io/cran/nlme/man/predict.nlme.html
【问题讨论】：
-
我在
pframe中看不到块列。newdata参数需要与原始data=参数一致。 -
是的，没错..但我希望在 loc 级别进行预测。这可行吗？
-
当然可以，但是您需要为
block指定一些级别。可以都是同一个级别，但需要是原始数据范围内的值。 -
好的。完美的。最后一个问题。如果我使用块 1 而不是块级别的块，结果会改变吗？
-
我预计它会改变，假设您的意思是“Block1”与“Block2”，即现有级别..您为什么会期待其他情况？