# vignettes/examples/hornSharkExample.R In momentuHMM: Maximum Likelihood Analysis of Animal Movement Behavior Using Multivariate Hidden Markov Models

```library(momentuHMM)
library(data.tree)

# coarse-scale data
data <- data.frame(ID=unlist(mapply(function(x) rep(paste0("seg",x),nrow(steps[[x]])-1),1:length(steps))),
level="1",
steps=unlist(lapply(steps,function(x) x\$steps[-nrow(x)])),
stepCat=unlist(lapply(steps,function(x) x\$cats[-nrow(x)])),
odbas=NA)
data\$stepCat[which(is.na(data\$steps))] <- NA

### add extra rows for fine-scale data
# level=1  covariate indicates when coarse scale behavior switching can occur (i.e., coarse scale t.p.m. used when level=1)
# level=2i covariate indicates start of each fine scale interval (i.e., fine scale initial distribution used when level=2i)
# level=2  otherwise (i.e., fine scale t.p.m used when level=2)
odbas <- unlist(odbas)
sharkData <- NULL
for(i in 1:nrow(data)){
fineInd <- data.frame(ID=data\$ID[i],level="2",steps=NA,stepCat=NA,odbas=odbas[(i-1)*50+1:50])
tmp <- rbind(data[i,,drop=FALSE],data.frame(ID=data\$ID[i],level="2i",steps=NA,stepCat=NA,odbas=NA),fineInd)
sharkData <- rbind(sharkData,tmp)
}

# prepare hierarchical data
sharkData <- prepData(sharkData,coordNames=NULL,hierLevels=c("1","2i","2"))

# data summary
summary(sharkData,dataNames=names(sharkData)[-1])

### define hierarchical HMM: states 1-3 = coarse state 1 (high activity); states 4-6 = coarse state 2 (resting); states 7-9 = coarse state 3 (travelling)
hierStates <- data.tree::Node\$new("shark HHMM states")
hierStates\$AddChild("activity")   # zero distance travelled, high activity
hierStates\$AddChild("resting")  # zero distance travelled, low activity

print(hierStates,"state")

nbStates <- length(hierStates\$Get("state",filterFun=data.tree::isLeaf))
nCat <- 8 # number of stepCat categories

# data stream distributions: level 1 = coarse level (stepCat="cat8"); level 2 = fine level (odbas="gamma")
hierDist <- data.tree::Node\$new("shark HHMM dist")
hierDist\$level1\$AddChild("stepCat", dist=paste0("cat",nCat)) # categorical distribution with 8 categories (1, 2, 3, 4, 5, 6, 7, or 8)
# equivalent
#hierDist <- data.tree::as.Node(list(name="shark HHMM dist",
#                           level1=list(stepCat=list(dist="cat8")),
#                           level2=list(odbas=list(dist="gamma"))))

print(hierDist,"dist")

### defining start values based on Adam et al
mu0 <- c(0.191, 0.323, 0.721, 0.084, 0.15, 0.228, 0.094, 0.191, 0.39) # real scale
sd0 <- c(0.047, 0.051, 0.248, 0.021, 0.025, 0.033, 0.026, 0.039, 0.159) # real scale
probs0 <- c(1e+10, -1e+10, 0.958, 5.064, 3.261, 4.021, 0.473, 2.568) # working scale

# constrain coarse-scale data stream parameters for corresponding fine-scale states
DM <- list(stepCat=matrix(cbind(c(rep(c(1,1,1),2),rep(0,(nCat-1)*nbStates-6)),kronecker(diag(nCat-1),c(0,0,0,0,0,0,1,1,1))),
nrow=(nCat-1)*nbStates,
ncol=8,
dimnames=list(paste0(rep(paste0("prob",1:(nCat-1),"_"),each=nbStates),1:nbStates),
c(paste0(c("prob1_12",paste0("prob",1:(nCat-1),"_3")),":(Intercept)")))))

Par0 <- list(stepCat=probs0,
odbas=c(mu0,sd0))
fixPar <- list(stepCat=c(1.e+10,-1.e+10,rep(NA,6)))

hierBeta <- data.tree::Node\$new("shark beta")
hierBeta\$AddChild("level1",beta=matrix(c(-0.651, -0.169, -1.884, -0.369, -1.596, -0.737),1,length(hierStates\$children)*(length(hierStates\$children)-1)))
hierBeta\$level2\$AddChild("activity",beta=matrix(c(-2.902, -14.032, 3.059, -1.37, 8.098, 11.66),1,length(hierStates\$activity\$children)*(length(hierStates\$activity\$children)-1)))
hierBeta\$level2\$AddChild("resting",beta=matrix(c(-3.264, -14.279, 3.107, 0.252, 13.861, 16.468),1,length(hierStates\$resting\$children)*(length(hierStates\$resting\$children)-1)))
hierBeta\$level2\$AddChild("transit",beta=matrix(c(-3.21, -21.32, 3.463, -0.598, 14.636, 17.811),1,length(hierStates\$transit\$children)*(length(hierStates\$transit\$children)-1)))

hierDelta <- data.tree::Node\$new("shark delta")

# check hierarchical model specification and parameters
checkPar0(sharkData,hierStates=hierStates,hierDist=hierDist,Par0=Par0,DM=DM,hierBeta=hierBeta,hierDelta=hierDelta,fixPar=fixPar)

# compute hierarchical state transition probabilities based on initial values
iTrProbs <- getTrProbs(sharkData,hierStates=hierStates,hierBeta=hierBeta,hierDist=hierDist)
iTrProbs\$level1\$gamma[,,1] # tpm at first time step for level1
lapply(iTrProbs\$level2,function(x) x\$gamma[,,1]) # tpm at first time step for level2

hhmm <- fitHMM(sharkData,hierStates=hierStates,hierDist=hierDist,Par0=Par0,hierBeta=hierBeta,hierDelta=hierDelta,DM=DM,fixPar=fixPar,
nlmPar=list(print.level=2,iterlim=1000,print.level=2,stepmax=25,hessian=FALSE))
hhmm
plotPR(hhmm)

### most likely state sequence #######################
states <- viterbi(hhmm)
hStates <- viterbi(hhmm,hierarchical=TRUE)
# coarse scale sequence
coarseStates <- hStates\$level1
# fine scale sequence
fineStates <- hStates\$level2
######################################################

# transition probabilities
trProbs12 <- getTrProbs(hhmm,covIndex=c(1,3))

# stationary distributions
stats12 <- stationary(hhmm,covIndex=c(1,3))

### coarse scale #####################################
# t.p.m.
round(trProbs12\$level1\$gamma[,,1],3)

# stationary distribution
round(stats12[[1]]\$level1[1,],3)
######################################################

### fine scale #######################################
# t.p.m.
lapply(trProbs12\$level2,function(x) round(x\$gamma[,,1],3))

# stationary distribution
lapply(stats12[[1]]\$level2,function(x) round(x[1,],3))
######################################################

### Simulate from fitted model ##############################
obsPerLevel <- data.tree::Node\$new("simHierData")
obsPerLevel\$AddChild("level1",obs=7) # number of level 1 observations
obsPerLevel\$AddChild("level2",obs=10)  # number of level 2 observations that follow each level 1 observation

simHHMM <- simHierData(model=hhmm, nbAnimals = length(unique(sharkData\$ID)), obsPerLevel = obsPerLevel, states = TRUE)