if (FALSE) {
# simulated data with normal random effects
# and binary individual covariate
nbAnimals <- 5 # should be larger for random effects estimation
obsPerAnimal <- 110
indCov <- rbinom(nbAnimals,1,0.5) # individual covariate
betaCov <- c(-0.5,0.5) # covariate effects
mu <- c(-0.1,0.1) # mean for random effects
sigma <- c(0.2,0.4) # sigma for random effects
beta0 <- cbind(rnorm(nbAnimals,mu[1],sigma[1]),
rnorm(nbAnimals,mu[2],sigma[2]))
reData <- simData(nbAnimals=nbAnimals,obsPerAnimal=obsPerAnimal,nbStates=2,
dist=list(step="gamma"),formula=~0+ID+indCov,
Par=list(step=c(1,10,1,2)),
beta=rbind(beta0,betaCov),
covs=data.frame(indCov=rep(indCov,each=obsPerAnimal)))
# fit null model
nullFit <- fitHMM(reData,nbStates=2,
dist=list(step="gamma"),
Par0=list(step=c(1,10,1,2)))
# fit covariate model
covFit <- fitHMM(reData,nbStates=2,
dist=list(step="gamma"),formula=~indCov,
Par0=list(step=c(1,10,1,2)),
beta0=rbind(mu,betaCov))
# fit fixed effects model
fixFit <- fitHMM(reData,nbStates=2,
dist=list(step="gamma"),formula=~0+ID,
Par0=list(step=c(1,10,1,2)),
beta0=beta0)
# fit random effect model
reFit <- randomEffects(fixFit)
# fit random effect model with individual covariate
reCovFit <- randomEffects(fixFit, Xformula=~indCov)
# compare by AICc
AIC(nullFit,covFit,fixFit,reFit,reCovFit, n=nrow(reData))
}
Run the code above in your browser using DataLab