if (FALSE) {
data(Soybean)
attach(Soybean)
#################################
#A full model (no covariate)
y = weight #response
x = Time #time
#Expression for the three parameter logistic curve
exprNL = expression((fixed[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))
#A median regression (by default)
median_reg = QRNLMM(y,x,Plot,initial,exprNL)
#Assing the fit
fxd = median_reg$res$beta
nlmodel = median_reg$res$nlmodel
seqc = seq(min(x),max(x),length.out = 500)
group.plot(x = Time,y = weight,groups = Plot,type="l",
main="Soybean profiles",xlab="time (days)",
ylab="mean leaf weight (gr)",col="gray")
lines(seqc,nlmodel(x = seqc,fixed = fxd,random = rep(0,3)),
lwd=2,col="blue")
#Histogram for residuals
hist(median_reg$res$residuals)
#########################################
#A model for comparing the two genotypes (with covariates)
y = weight #response
x = Time #time
covar = c(Variety)-1 #factor genotype (0=Forrest, 1=Plan Introduction)
#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 the three quartiles
box_reg = QRNLMM(y,x,Plot,initial,exprNL,covar,p=c(0.25,0.50,0.75))
#Assing the fit for the median (second quartile)
fxd = box_reg[[2]]$res$beta
nlmodel = box_reg[[2]]$res$nlmodel
seqc = seq(min(x),max(x),length.out = 500)
group.plot(x = Time[Variety=="P"],y = weight[Variety=="P"],
groups = Plot[Variety=="P"],type="l",col="light blue",
main="Soybean profiles by genotype",xlab="time (days)",
ylab="mean leaf weight (gr)")
group.lines(x = Time[Variety=="F"],y = weight[Variety=="F"],
groups = Plot[Variety=="F"],col="gray")
lines(seqc,nlmodel(x = seqc,fixed = fxd,random = rep(0,3),covar=1),
lwd=2,col="blue")
lines(seqc,nlmodel(x = seqc,fixed = fxd,random = rep(0,3),covar=0),
lwd=2,col="black")
}
Run the code above in your browser using DataLab