# NOT RUN {
nbStates <- 2
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
# extract data from momentuHMM example
data <- example$m$data
### 1. fit the model to the simulated data
# define initial values for the parameters
mu0 <- c(20,70)
sigma0 <- c(10,30)
kappa0 <- c(1,1)
stepPar <- c(mu0,sigma0) # no zero-inflation, so no zero-mass included
anglePar <- kappa0 # not estimating angle mean, so not included
formula <- ~cov1+cos(cov2)
m <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar,angle=anglePar),formula=formula)
print(m)
# }
# NOT RUN {
### 2. fit the exact same model to the simulated data using DM formulas
# Get initial values for the parameters on working scale
Par0 <- getParDM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par=list(step=stepPar,angle=anglePar),
DM=list(step=list(mean=~1,sd=~1),angle=list(concentration=~1)))
mDMf <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par0=Par0,formula=formula,
DM=list(step=list(mean=~1,sd=~1),angle=list(concentration=~1)))
print(mDMf)
### 3. fit the exact same model to the simulated data using DM matrices
# define DM
DMm <- list(step=diag(4),angle=diag(2))
# user-specified dimnames not required but are recommended
dimnames(DMm$step) <- list(c("mean_1","mean_2","sd_1","sd_2"),
c("mean_1:(Intercept)","mean_2:(Intercept)",
"sd_1:(Intercept)","sd_2:(Intercept)"))
dimnames(DMm$angle) <- list(c("concentration_1","concentration_2"),
c("concentration_1:(Intercept)","concentration_2:(Intercept)"))
mDMm <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par0=Par0,formula=formula,
DM=DMm)
print(mDMm)
### 4. fit step mean parameter covariate model to the simulated data using DM
stepDMf <- list(mean=~cov1,sd=~1)
Par0 <- getParDM(data,nbStates,list(step=stepDist,angle=angleDist),
Par=list(step=stepPar,angle=anglePar),
DM=list(step=stepDMf,angle=DMm$angle))
mDMfcov <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par0=Par0,
formula=formula,
DM=list(step=stepDMf,angle=DMm$angle))
print(mDMfcov)
### 5. fit the exact same step mean parameter covariate model using DM matrix
stepDMm <- matrix(c(1,0,0,0,"cov1",0,0,0,0,1,0,0,0,"cov1",0,0,
0,0,1,0,0,0,0,1),4,6,dimnames=list(c("mean_1","mean_2","sd_1","sd_2"),
c("mean_1:(Intercept)","mean_1:cov1","mean_2:(Intercept)","mean_2:cov1",
"sd_1:(Intercept)","sd_2:(Intercept)")))
Par0 <- getParDM(data,nbStates,list(step=stepDist,angle=angleDist),
Par=list(step=stepPar,angle=anglePar),
DM=list(step=stepDMm,angle=DMm$angle))
mDMmcov <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par0=Par0,
formula=formula,
DM=list(step=stepDMm,angle=DMm$angle))
print(mDMmcov)
### 6. fit circular-circular angle mean covariate model to the simulated data using DM
# Generate fake circular covariate using circAngles
data$cov3 <- circAngles(refAngle=2*atan(rnorm(nrow(data))),data)
# Fit circular-circular regression model for angle mean
# Note no intercepts are estimated for angle means because these are by default
# the previous movement direction (i.e., a turning angle of zero)
mDMcircf <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar,angle=c(0,0,Par0$angle)),
formula=formula,
estAngleMean=list(angle=TRUE),
circularAngleMean=list(angle=TRUE),
DM=list(angle=list(mean=~cov3,concentration=~1)))
print(mDMcircf)
### 7. fit the exact same circular-circular angle mean model using DM matrices
# Note no intercept terms are included in DM for angle means because the intercept is
# by default the previous movement direction (i.e., a turning angle of zero)
mDMcircm <- fitHMM(data=data,nbStates=nbStates,dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar,angle=c(0,0,Par0$angle)),
formula=formula,
estAngleMean=list(angle=TRUE),
circularAngleMean=list(angle=TRUE),
DM=list(angle=matrix(c("cov3",0,0,0,0,"cov3",0,0,0,0,1,0,0,0,0,1),4,4)))
print(mDMcircm)
### 8. Cosinor and state-dependent formulas
nbStates<-2
dist<-list(step="gamma")
Par<-list(step=c(100,1000,50,100))
# include 24-hour cycle on all transition probabilities
# include 12-hour cycle on transitions from state 2
formula=~cosinor(hour24,24)+state2(cosinor(hour12,12))
# specify appropriate covariates
covs<-data.frame(hour24=0:23,hour12=0:11)
beta<-matrix(c(-1.5,1,1,NA,NA,-1.5,-1,-1,1,1),5,2)
# row names for beta not required but can be helpful
rownames(beta)<-c("(Intercept)",
"cosinorCos(hour24, 24)",
"cosinorSin(hour24, 24)",
"cosinorCos(hour12, 12)",
"cosinorSin(hour12, 12)")
data.cos<-simData(nbStates=nbStates,dist=dist,Par=Par,
beta=beta,formula=formula,covs=covs)
m.cosinor<-fitHMM(data.cos,nbStates=nbStates,dist=dist,Par0=Par,formula=formula)
m.cosinor
### 9. Piecewise constant B-spline on step length mean and angle concentration
nObs <- 1000 # length of simulated track
cov <- data.frame(time=1:nObs) # time covariate for splines
dist <- list(step="gamma",angle="vm")
stepDM <- list(mean=~bSpline(time,df=3,degree=0),sd=~1)
angleDM <- list(mean=~1,concentration=~bSpline(time,df=3,degree=0))
DM <- list(step=stepDM,angle=angleDM)
Par <- list(step=c(log(1000),1,-1,log(100)),angle=c(0,log(10),2,-5))
data.spline<-simData(obsPerAnimal=nObs,nbStates=1,dist=dist,Par=Par,DM=DM,covs=cov)
Par0 <- list(step=Par$step,angle=Par$angle[-1])
m.spline<-fitHMM(data.spline,nbStates=1,dist=dist,Par0=Par0,
DM=list(step=stepDM,
angle=list(concentration=~bSpline(time,df=3,degree=0))))
### 10. Initial state (delta) based on covariate
nObs <- 100
dist <- list(step="gamma",angle="vm")
Par <- list(step=c(100,1000,50,100),angle=c(0,0,0.01,0.75))
# create sex covariate
cov <- data.frame(sex=factor(rep(c("F","M"),each=nObs))) # sex covariate
formulaDelta <- ~ sex + 0
# Female begins in state 1, male begins in state 2
delta <- matrix(c(-100,100),2,1,dimnames=list(c("sexF","sexM"),"state 2"))
data.delta<-simData(nbAnimals=2,obsPerAnimal=nObs,nbStates=2,dist=dist,Par=Par,
delta=delta,formulaDelta=formulaDelta,covs=cov)
Par0 <- list(step=Par$step, angle=Par$angle[3:4])
m.delta <- fitHMM(data.delta, nbStates=2, dist=dist, Par0 = Par0,
formulaDelta=formulaDelta)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab