I have a simple experiment (here is some made-up data) with 3 locations ("loc"), blocks within locations ("block"), a set of 8 treatments ("treat), and a response variable ("response"). An exponential plateau function was adjusted to this data (response as a function of treat) and the NLME function was used to analyze the whole experiment using mixed models. The parameters of the exponential plateau function is considered the fixed part while loc and block are random.
My problem is with the prediction of the model. I am not able to get the prediction for the response to the treatments at the loc level. I am able to get it to the population level (fixed effect), but got all NAs for the prediction when try to do it a the loc level. Is there something wrong with the way I set the "levels= " option?
Here is my made up data
#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")
Here is the function and the model fit
#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)
And here is how I made prediction
#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)
Prediction with level=0 works fine for the prediction of fixed effects. Prediction with level=2 works fine for the prediction at block/loc random level Prediction with level=1 gives me NAs.. this would be the prediction at the loc level, which is the one I am most interested in.
Is there any mistake in the level option?
Here is the explanation of prediction using the nlme function https://rdrr.io/cran/nlme/man/predict.nlme.html
