data(cow.saemix)
saemix.data<-saemixData(name.data=cow.saemix,header=TRUE,name.group=c("cow"),
name.predictors=c("time"),name.response=c("weight"),
name.covariates=c("birthyear","twin","birthrank"),
units=list(x="days",y="kg",covariates=c("yr","-","-")))
growthcow<-function(psi,id,xidep) {
x<-xidep[,1]
a<-psi[id,1]
b<-psi[id,2]
k<-psi[id,3]
f<-a*(1-b*exp(-k*x))
return(f)
}
saemix.model<-saemixModel(model=growthcow,
description="Exponential growth model",
psi0=matrix(c(700,0.9,0.02,0,0,0),ncol=3,byrow=TRUE,
dimnames=list(NULL,c("A","B","k"))),transform.par=c(1,1,1),fixed.estim=c(1,1,1),
covariate.model=matrix(c(0,0,0),ncol=3,byrow=TRUE),
covariance.model=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE),
omega.init=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE),error.model="constant")
saemix.options<-list(algorithms=c(1,1,1),nb.chains=1,nbiter.saemix=c(200,100),
seed=4526,save=FALSE,save.graphs=FALSE,displayProgress=FALSE)
# Plotting the data
plot(saemix.data,xlab="Time (day)",ylab="Weight of the cow (kg)")
# \donttest{
saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)
# }
Run the code above in your browser using DataLab