if (FALSE) {
#A model for comparing the two genotypes (with covariates)
data(Soybean)
attach(Soybean)
y = weight #response
x = Time #time
covar = c(Variety)-1 #factor genotype (0=Forrest, 1=Plan Intro)
#Expression for the three parameter logistic curve with a covariate
exprNL = expression((fixed[1]+(fixed[4]*covar[1])+random[1])/
(1 + exp(((fixed[2]+random[2])- x)/
(fixed[3]+random[3]))))
#Initial values for fixed effects
initial = c(max(y),0.6*max(y),0.73*max(y),3)
# A quantile regression for percentiles p = c(0.05,0.50,0.95)
#Take your sit and some popcorn
results = qrNLMM::QRNLMM(
y = y,
x = x,
groups = Plot,
initial = initial,
exprNL = exprNL,
covar = covar,
p= c(0.05,0.50,0.95),# quantiles to estimate
MaxIter = 50,M = 15, # to accelerate
verbose = FALSE # show no output
)
######################################################################
# Predicting
######################################################################
# now we select two random subjects from original data
set.seed(19)
index = Plot %in% sample(Plot,size = 2)
index2 = c(1,diff(as.numeric(Plot)))>0 & index
# 1. Original dataset
#####################
prediction = predict(object = results)
head(prediction) # fitted values
if(TRUE){
group.plot(x = Time[index],
y = weight[index],
groups = Plot[index],
type="b",
main="Soybean profiles",
xlab="time (days)",
ylab="mean leaf weight (gr)",
col= ifelse(covar[index2],"gray70","gray90"),
ylim = range(prediction[index,]),
lty = 1
)
legend("bottomright",
legend = c("Forrest","Plan Intro"),
bty = "n",col = c(4,2),lty = 1)
# predictions for these two plots
group.lines(x = Time[index], # percentile 5
y = prediction[index,1],
groups = Plot[index],
type = "l",
col=ifelse(covar[index2],"red","blue"),
lty = 2
)
group.lines(x = Time[index], # median
y = prediction[index,2],
groups = Plot[index],
type = "l",
col="black",
lty = 2)
group.lines(x = Time[index], # percentile 95
y = prediction[index,3],
groups = Plot[index],
type = "l",
col=ifelse(covar[index2],"red","blue"),
lty = 2)
legend("topleft",
legend = c("p = 5","p = 50","p = 95"),
col = c(4,1,4),lty = c(2,2,2),bty = "n")
}
# 2. New covariates with no response
####################################
# For the two randomly selected plots (index == TRUE)
# newdata
newdata = data.frame(new.groups = Plot[index],
new.x = Time[index],
new.covar = covar[index])
newdata
attach(newdata)
prediction2 = predict(object = results,
groups = new.groups,
x = new.x,
covar = new.covar)
# population curves
if(TRUE){
group.plot(x = new.x, # percentile 5
y = prediction2[,1],
groups = new.groups,
type = "l",
col=ifelse(covar[index2],"red","blue"),
lty = 2,
main="Soybean profiles",
xlab="time (days)",
ylab="mean leaf weight (gr)",
ylim = range(prediction2)
)
legend("bottomright",
legend = c("Forrest","Plan Intro"),
bty = "n",col = c(4,2),lty = 1)
# predictions for these two plots
group.lines(x = new.x, # median
y = prediction2[,2],
groups = new.groups,
type = "l",
col="black",
lty = 1)
group.lines(x = new.x, # percentile 95
y = prediction2[,3],
groups = new.groups,
type = "l",
col=ifelse(covar[index2],"red","blue"),
lty = 2)
legend("topleft",
legend = c("p = 5","p = 50","p = 95"),
col = c(4,1,4),lty = c(2,1,2),bty = "n")
segments(x0 = new.x[new.covar==1],
y0 = prediction2[new.covar==1,1],
y1 = prediction2[new.covar==1,3],
lty=2,col = "red")
segments(x0 = new.x[new.covar==0],
y0 = prediction2[new.covar==0,1],
y1 = prediction2[new.covar==0,3],
lty=2,col = "blue")
}
# 3. New covariates + response
####################################
# newdata
newdata2 = data.frame(new.groups = Plot[index],
new.x = Time[index],
new.covar = covar[index],
new.y = weight[index])
newdata2
attach(newdata2)
prediction2 = predict(object = results,
groups = new.groups,
x = new.x,
covar = new.covar,
y = new.y)
# individual curves (random-effects to be computed)
if(TRUE){
group.plot(x = Time[index],
y = weight[index],
groups = Plot[index],
type="b",
main="Soybean profiles",
xlab="time (days)",
ylab="mean leaf weight (gr)",
col= ifelse(covar[index2],"gray70","gray90"),
ylim = range(prediction[index,]),
lty = 1
)
legend("bottomright",
legend = c("Forrest","Plan Intro"),
bty = "n",col = c(4,2),lty = 1)
# predictions for these two plots
group.lines(x = new.x, # percentile 5
y = prediction2[,1],
groups = new.groups,
type = "l",
col=ifelse(covar[index2],"red","blue"),
lty = 2)
group.lines(x = new.x, # median
y = prediction2[,2],
groups = new.groups,
type = "l",
col="black",
lty = 1)
group.lines(x = new.x, # percentile 95
y = prediction2[,3],
groups = new.groups,
type = "l",
col=ifelse(covar[index2],"red","blue"),
lty = 2)
legend("topleft",
legend = c("p = 5","p = 50","p = 95"),
col = c(4,1,4),lty = c(2,1,2),bty = "n")
}
}
Run the code above in your browser using DataLab