cow.saemix | R Documentation |
The cow.saemix
data contains records of the weight of 560 cows on 9 or 10 occasions.
cow.saemix
This data frame contains the following columns:
the unique identifier for each cow
time (days)
a numeric vector giving the weight of the cow (kg)
year of birth (between 1988 and 1998)
existence of a twin (no=1, yes=2)
the rank of birth (beetween 3 and 7)
An exponential model was assumed to describe the weight gain with time: y_ij = A_i (1- B_i exp( - K_i t_ij)) +epsilon_ij
JC Pinheiro, DM Bates (2000) Mixed-effects Models in S and S-PLUS, Springer, New York (Appendix A.19)
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)")
saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.