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

#### Documented in CIrealCIreal.defaultCIreal.hierarchical

#' Confidence intervals for the natural (i.e., real) parameters
#'
#' Computes the standard errors and confidence intervals on the real (i.e., natural) scale of the data stream probability distribution parameters,
#' as well as for the transition probabilities parameters. If covariates are included in the probability distributions or TPM formula, the mean values
#' of non-factor covariates are used for calculating the natural parameters. For any covariate(s) of class 'factor', then the value(s) from the first observation
#' in the data are used.
#'
#' @param m A \code{momentuHMM}, \code{momentuHierHMM}, \code{miHMM}, or \code{miSum} object
#' @param alpha Significance level of the confidence intervals. Default: 0.95 (i.e. 95\% CIs).
#' @param covs Data frame consisting of a single row indicating the covariate values to be used in the calculations. By default, no covariates are specified.
#' @param parms Optional character vector indicating which groups of real parameters to calculate confidence intervals for (e.g., 'step', 'angle', 'gamma', 'delta', etc.). Default: NULL, in which case confidence intervals are calculated for all groups of parameters in the model.
#'
#' @details
#' For any covariates that are not specified using \code{covs}, the means of the covariate(s) are used
#' (unless the covariate is a factor, in which case the first factor in the data is used).
#'
#' @return A list of the following objects:
#' \item{...}{List(s) of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the natural parameters of the data streams}
#' \item{gamma}{List of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the transition probabilities}
#' \item{delta}{List of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the initial state probabilities}
#' \item{hierGamma}{A hierarchical data structure \code{\link[data.tree]{Node}} including a list of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the transition probabilities for each level of the hierarchy (only applies if \code{m} is a hierarchical model object)}
#' \item{hierDelta}{A hierarchical data structure \code{\link[data.tree]{Node}} including a list of estimates ('est'), standard errors ('se'), and confidence intervals ('lower', 'upper') for the initial state probabilities for each level of the hierarchy (only applies if \code{m} is a hierarchical model object)}
#'
#' @examples
#' # m is a momentuHMM object (as returned by fitHMM), automatically loaded with the package
#' m <- example$m #' #' ci1<-CIreal(m) #' #' # specify 'covs' #' ci2<-CIreal(m,covs=data.frame(cov1=mean(m$data$cov1),cov2=mean(m$data$cov2))) #' #' all.equal(ci1,ci2) #' #' @export CIreal <- function(m,alpha=0.95,covs=NULL,parms=NULL) { UseMethod("CIreal") } #' @rdname CIreal #' @method CIreal default #' @importFrom numDeriv grad #' @importFrom utils tail #' @importFrom Brobdingnag as.brob sum #' @importFrom CircStats circ.mean #' @export CIreal.default <- function(m,alpha=0.95,covs=NULL,parms=NULL) { if(is.miSum(m) | is.miHMM(m)){ if(is.miHMM(m)) m <- m$miSum
m <- formatmiSum(m)
m$mod$hessian <- NA
m$mod$Sigma <- matrix(0,length(m$mod$estimate),length(m$mod$estimate))
if(length(m$conditions$optInd)) m$mod$Sigma[-m$conditions$optInd,-m$conditions$optInd] <- m$MIcombine$variance
else m$mod$Sigma <- m$MIcombine$variance
m$CIreal <- m$Par$real m <- momentuHMM(m) } else if(!is.momentuHMM(m)) stop("'m' must be a momentuHMM, miSum, or miHMM object") if(!is.null(m$conditions$fit) && !m$conditions$fit) stop("The given model hasn't been fitted.") if(alpha<0 | alpha>1) stop("alpha needs to be between 0 and 1.") nbStates <- length(m$stateNames)
nbAnimals <- length(unique(m$data$ID))

dist <- m$conditions$dist
distnames <- names(dist)
DMind <- m$conditions$DMind

if(is.null(parms)) pparms <- c(distnames,"gamma","delta")
else pparms <- parms

m <- delta_bc(m)

tempCovs <- getCovs(m,covs,unique(m$data$ID),checkHier=FALSE)[1,]

# identify covariates
reForm <- formatRecharge(nbStates,m$conditions$formula,m$conditions$betaRef,m$data,covs=tempCovs,par=m$mle)
m$data <- cbind(m$data,reForm$newdata) recharge <- reForm$recharge
hierRecharge <- reForm$hierRecharge newformula <- reForm$newformula
tempCovs <- reForm$covs nbCovs <- reForm$nbCovs
if(!is.null(recharge) & is.null(parms)) pparms <- c(pparms,"g0","theta")

# inverse of Hessian
if(!is.null(m$mod$hessian) && !inherits(m$mod$Sigma,"error")) Sigma <- m$mod$Sigma
else Sigma <- NULL

ncmean <- get_ncmean(distnames,m$conditions$fullDM,m$conditions$circularAngleMean,nbStates)
nc <- ncmean$nc meanind <- ncmean$meanind

tmPar <- lapply(m$mle[distnames],function(x) c(t(x))) parCount<- lapply(m$conditions$fullDM,ncol) for(i in distnames[!unlist(lapply(m$conditions$circularAngleMean,isFALSE))]){ parCount[[i]] <- length(unique(gsub("cos","",gsub("sin","",colnames(m$conditions$fullDM[[i]]))))) } parindex <- c(0,cumsum(unlist(parCount))[-length(m$conditions$fullDM)]) names(parindex) <- distnames for(i in distnames){ if(!is.null(m$conditions$DM[[i]])){# & m$conditions$DMind[[i]]){ tmPar[[i]] <- m$mod$estimate[parindex[[i]]+1:parCount[[i]]] if(!isFALSE(m$conditions$circularAngleMean[[i]])){ names(tmPar[[i]]) <- unique(gsub("cos","",gsub("sin","",colnames(m$conditions$fullDM[[i]])))) } else names(tmPar[[i]])<-colnames(m$conditions$fullDM[[i]]) } else{ if(dist[[i]] %in% angledists) if(!m$conditions$estAngleMean[[i]]) tmPar[[i]] <- tmPar[[i]][-(1:nbStates)] } } Par <- list() lower<-list() upper<-list() se<-list() inputs <- checkInputs(nbStates,dist,tmPar,m$conditions$estAngleMean,m$conditions$circularAngleMean,m$conditions$zeroInflation,m$conditions$oneInflation,m$conditions$DM,m$conditions$userBounds,m$stateNames)
p<-inputs$p splineInputs<-getSplineDM(distnames,inputs$DM,m,tempCovs)
covs<-splineInputs$covs DMinputs<-getDM(covs,splineInputs$DM,inputs$dist,nbStates,p$parNames,p$bounds,tmPar,m$conditions$zeroInflation,m$conditions$oneInflation,m$conditions$circularAngleMean) fullDM <- DMinputs$fullDM
#DMind <- DMinputs$DMind #DMinputs<-getDM(tempCovs,inputs$DM,inputs$dist,nbStates,p$parNames,p$bounds,tmPar,m$conditions$zeroInflation,m$conditions$oneInflation,m$conditions$circularAngleMean) #fullDM<-DMinputs$fullDM

for(i in distnames[which(distnames %in% pparms)]){
tmpParNames <- p$parNames[[i]] tmpParNames[which(p$parNames[[i]]=="kappa")] <- "concentration"
if(!m$conditions$DMind[[i]]){
par <- c(w2n(m$mod$estimate,p$bounds,p$parSize,nbStates,nbCovs,m$conditions$estAngleMean,m$conditions$circularAngleMean[i],inputs$consensus[i],m$conditions$stationary,fullDM,m$conditions$DMind,1,inputs$dist[i],m$conditions$Bndind,nc,meanind,m$covsDelta,m$conditions$workBounds,m$covsPi)[[i]])
} else {
par <- as.vector(t(m$mle[[i]])) } if(!(inputs$dist[[i]] %in% angledists) | (inputs$dist[[i]] %in% angledists & m$conditions$estAngleMean[[i]] & !m$conditions$Bndind[[i]])) { Par[[i]] <- get_CI(m$mod$estimate,par,m,parindex[[i]]+1:parCount[[i]],fullDM[[i]],DMind[[i]],p$bounds[[i]],Sigma,m$conditions$circularAngleMean[[i]],inputs$consensus[[i]],nbStates,alpha,tmpParNames,m$stateNames,nc[[i]],meanind[[i]],m$conditions$workBounds[[i]],inputs$dist[[i]]) } else { if(!m$conditions$estAngleMean[[i]]) Par[[i]] <- get_CI(m$mod$estimate,par[-c(1:nbStates)],m,parindex[[i]]+1:parCount[[i]],fullDM[[i]],DMind[[i]],p$bounds[[i]],Sigma,m$conditions$circularAngleMean[[i]],inputs$consensus[[i]],nbStates,alpha,tmpParNames,m$stateNames,nc[[i]],meanind[[i]],m$conditions$workBounds[[i]],inputs$dist[[i]]) else { if(m$conditions$Bndind[[i]]){ Par[[i]] <- CI_angle(m$mod$estimate,par,m,parindex[[i]]+1:parCount[[i]],fullDM[[i]],DMind[[i]],p$bounds[[i]],Sigma,m$conditions$circularAngleMean[[i]],inputs$consensus[[i]],nbStates,alpha,tmpParNames,m$stateNames,nc[[i]],meanind[[i]],m$conditions$workBounds[[i]],inputs$dist[[i]]) } } } } mixtures <- m$conditions$mixtures if(mixtures>1 & is.null(parms)) pparms <- c(pparms,"pi") if(nbStates>1) { if("gamma" %in% pparms){ # identify parameters of interest i2 <- tail(cumsum(unlist(parCount)),1)+1 i3 <- i2+nbStates*(nbStates-1)*(nbCovs+1)*mixtures-1 quantSup <- qnorm(1-(1-alpha)/2) tmpSplineInputs<-getSplineFormula(newformula,m$data,tempCovs)
tempCovMat <- stats::model.matrix(tmpSplineInputs$formula,data=tmpSplineInputs$covs)

est<-lower<-upper<-se<-matrix(NA,nbStates*mixtures,nbStates)
for(mix in 1:mixtures){

if(is.null(recharge)){
wpar <- m$mod$estimate[i2:i3][unique(c(m$conditions$betaCons))]
est[(mix-1)*nbStates+1:nbStates,] <- get_gamma(wpar,tempCovMat,nbStates,betaRef=m$conditions$betaRef,betaCons=m$conditions$betaCons,workBounds=m$conditions$workBounds$beta,mixture=mix) tmpSig <- Sigma[(i2:i3)[unique(c(m$conditions$betaCons))],(i2:i3)[unique(c(m$conditions$betaCons))]] } else { wpar <- c(m$mod$estimate[i2:i3][unique(c(m$conditions$betaCons))],m$mod$estimate[length(m$mod$estimate)-reForm$nbRecovs:0])
est[(mix-1)*nbStates+1:nbStates,] <- get_gamma_recharge(wpar,tmpSplineInputs$covs,tmpSplineInputs$formula,hierRecharge,nbStates,betaRef=m$conditions$betaRef,betaCons=m$conditions$betaCons,workBounds=rbind(m$conditions$workBounds$beta,m$conditions$workBounds$theta),mixture=mix)
tmpSig <- Sigma[c((i2:i3)[unique(c(m$conditions$betaCons))],length(m$mod$estimate)-reForm$nbRecovs:0),c((i2:i3)[unique(c(m$conditions$betaCons))],length(m$mod$estimate)-reForm$nbRecovs:0)]
}

if(!is.null(Sigma)){
for(i in 1:nbStates){
for(j in 1:nbStates){
if(is.null(recharge)){
dN<-numDeriv::grad(get_gamma,wpar,covs=tempCovMat,nbStates=nbStates,i=i,j=j,betaRef=m$conditions$betaRef,betaCons=m$conditions$betaCons,workBounds=m$conditions$workBounds$beta,mixture=mix) } else { dN<-numDeriv::grad(get_gamma_recharge,wpar,covs=tmpSplineInputs$covs,formula=tmpSplineInputs$formula,hierRecharge=hierRecharge,nbStates=nbStates,i=i,j=j,betaRef=m$conditions$betaRef,betaCons=m$conditions$betaCons,workBounds=rbind(m$conditions$workBounds$beta,m$conditions$workBounds$theta),mixture=mix) } se[(mix-1)*nbStates+i,j]<-suppressWarnings(sqrt(dN%*%tmpSig%*%dN)) lower[(mix-1)*nbStates+i,j]<-1/(1+exp(-(log(est[(mix-1)*nbStates+i,j]/(1-est[(mix-1)*nbStates+i,j]))-quantSup*(1/(est[(mix-1)*nbStates+i,j]-est[(mix-1)*nbStates+i,j]^2))*se[(mix-1)*nbStates+i,j])))#est[i,j]-quantSup*se[i,j] upper[(mix-1)*nbStates+i,j]<-1/(1+exp(-(log(est[(mix-1)*nbStates+i,j]/(1-est[(mix-1)*nbStates+i,j]))+quantSup*(1/(est[(mix-1)*nbStates+i,j]-est[(mix-1)*nbStates+i,j]^2))*se[(mix-1)*nbStates+i,j])))#m$mle$gamma[i,j]+quantSup*se[i,j] } } } } Par$gamma <- list(est=est,se=se,lower=lower,upper=upper)
dimnames(Par$gamma$est) <- dimnames(Par$gamma$se) <- dimnames(Par$gamma$lower) <- dimnames(Par$gamma$upper) <- list(rep(m$stateNames,mixtures),m$stateNames)
if(mixtures>1) dimnames(Par$gamma$est) <- dimnames(Par$gamma$se) <- dimnames(Par$gamma$lower) <- dimnames(Par$gamma$upper) <- list(paste0(rep(m$stateNames,mixtures),"_mix",rep(1:mixtures,each=nbStates)),m$stateNames)
}

# pi
if(mixtures>1 & ("pi" %in% pparms)){
wpar<-m$mod$estimate
pie <- matrix(wpar[i3+1:(ncol(m$covsPi)*(mixtures-1))],nrow=ncol(m$covsPi),ncol=mixtures-1)
tmpSig <- Sigma[i3+1:(ncol(m$covsPi)*(mixtures-1)),i3+1:(ncol(m$covsPi)*(mixtures-1))]
quantSup <- qnorm(1-(1-alpha)/2)
lower<-upper<-se<-matrix(NA,nrow=nrow(m$covsPi),ncol=mixtures) est<-matrix(m$mle[["pi"]],nrow=nrow(m$covsPi),ncol=mixtures) if(!is.null(Sigma)){ for(j in 1:nrow(m$covsPi)){
for(i in 1:mixtures){
dN<-numDeriv::grad(get_delta,pie,covsDelta=m$covsPi[j,,drop=FALSE],i=i,workBounds=m$conditions$workBounds[["pi"]]) se[j,i]<-suppressWarnings(sqrt(dN%*%tmpSig%*%dN)) lower[j,i]<-1/(1+exp(-(log(est[j,i]/(1-est[j,i]))-quantSup*(1/(est[j,i]-est[j,i]^2))*se[j,i]))) upper[j,i]<-1/(1+exp(-(log(est[j,i]/(1-est[j,i]))+quantSup*(1/(est[j,i]-est[j,i]^2))*se[j,i]))) } } } Par[["pi"]] <- list(est=est,se=se,lower=lower,upper=upper) colnames(Par[["pi"]]$est) <- colnames(Par[["pi"]]$se) <- colnames(Par[["pi"]]$lower) <- colnames(Par[["pi"]]$upper) <- paste0("mix",1:mixtures) rownames(Par[["pi"]]$est) <- rownames(Par[["pi"]]$se) <- rownames(Par[["pi"]]$lower) <- rownames(Par[["pi"]]$upper) <- paste0("ID:",unique(m$data$ID)) } } # delta if("delta" %in% pparms){ if(nbStates>1){ if(!m$conditions$stationary){ wpar<-m$mod$estimate nbCovsDelta <- ncol(m$covsDelta)-1
foo <- length(wpar)-ifelse(reForm$nbRecovs,(reForm$nbRecovs+1)+(reForm$nbG0covs+1),0)-(nbCovsDelta+1)*(nbStates-1)*mixtures delta <- matrix(wpar[foo+1:((nbCovsDelta+1)*(nbStates-1)*mixtures)],nrow=(nbCovsDelta+1)*mixtures,ncol=nbStates-1) tmpSig <- Sigma[foo+1:((nbCovsDelta+1)*(nbStates-1)*mixtures),foo+1:((nbCovsDelta+1)*(nbStates-1)*mixtures)] quantSup <- qnorm(1-(1-alpha)/2) lower<-upper<-se<-matrix(NA,nrow=nrow(m$covsDelta)*mixtures,ncol=nbStates)
est<-matrix(m$mle$delta,nrow=nrow(m$covsDelta)*mixtures,ncol=nbStates) for(mix in 1:mixtures){ if(!is.null(Sigma)){ for(j in 1:nrow(m$covsDelta)){
for(i in 1:nbStates){
dN<-numDeriv::grad(get_delta,delta,covsDelta=m$covsDelta[j,,drop=FALSE],i=i,workBounds=m$conditions$workBounds$delta,mixture=mix)
se[(mix-1)*nrow(m$covsDelta)+j,i]<-suppressWarnings(sqrt(dN%*%tmpSig%*%dN)) lower[(mix-1)*nrow(m$covsDelta)+j,i]<-1/(1+exp(-(log(est[(mix-1)*nrow(m$covsDelta)+j,i]/(1-est[(mix-1)*nrow(m$covsDelta)+j,i]))-quantSup*(1/(est[(mix-1)*nrow(m$covsDelta)+j,i]-est[(mix-1)*nrow(m$covsDelta)+j,i]^2))*se[(mix-1)*nrow(m$covsDelta)+j,i])))#est[(mix-1)*nrow(m$covsDelta)+j,i]-quantSup*se[i]
upper[(mix-1)*nrow(m$covsDelta)+j,i]<-1/(1+exp(-(log(est[(mix-1)*nrow(m$covsDelta)+j,i]/(1-est[(mix-1)*nrow(m$covsDelta)+j,i]))+quantSup*(1/(est[(mix-1)*nrow(m$covsDelta)+j,i]-est[(mix-1)*nrow(m$covsDelta)+j,i]^2))*se[(mix-1)*nrow(m$covsDelta)+j,i])))#est[j,i]+quantSup*se[i]
}
}
}
}
} else {
covs<-tempCovMat
statFun<-function(beta,nbStates,covs,i,mixture=1){
gamma <- trMatrix_rcpp(nbStates,beta[(mixture-1)*ncol(covs)+1:ncol(covs),,drop=FALSE],covs,m$conditions$betaRef)[,,1]
tryCatch(solve(t(diag(nbStates)-gamma+1),rep(1,nbStates))[i],error = function(e) {
"A problem occurred in the calculation of the stationary distribution."})
}
est <- lower <- upper <- se <- matrix(NA,nbAnimals*mixtures,nbStates)
wpar <- m$mod$estimate[i2:i3][unique(c(m$conditions$betaCons))]
tmpSig <- Sigma[(i2:i3)[unique(c(m$conditions$betaCons))],(i2:i3)[unique(c(m$conditions$betaCons))]]

for(mix in 1:mixtures){
delta <- statFun(matrix(wpar,nrow=(nbCovs+1)*mixtures),nbStates,covs,1:nbStates,mixture=mix)
est[nbAnimals*(mix-1)+1:nbAnimals,] <- matrix(delta,nrow=nbAnimals,ncol=nbStates,byrow=TRUE)
if(!is.null(Sigma)){
for(k in 1:nbStates){
se[nbAnimals*(mix-1)+1:nbAnimals,k]<-suppressWarnings(sqrt(dN%*%tmpSig%*%dN))
lower[nbAnimals*(mix-1)+1:nbAnimals,k] <- probCI(est[nbAnimals*(mix-1)+1:nbAnimals,k],se[nbAnimals*(mix-1)+1:nbAnimals,k],quantSup,bound="lower")
upper[nbAnimals*(mix-1)+1:nbAnimals,k] <- probCI(est[nbAnimals*(mix-1)+1:nbAnimals,k],se[nbAnimals*(mix-1)+1:nbAnimals,k],quantSup,bound="upper")
}
}
}
}
Par$delta <- list(est=est,se=se,lower=lower,upper=upper) } else { Par$delta <- list(est=matrix(1,nrow(m$covsDelta)),se=matrix(NA,nrow(m$covsDelta)),lower=matrix(NA,nrow(m$covsDelta)),upper=matrix(NA,nrow(m$covsDelta)))
}
colnames(Par$delta$est) <- colnames(Par$delta$se) <- colnames(Par$delta$lower) <- colnames(Par$delta$upper) <- m$stateNames rownames(Par$delta$est) <- rownames(Par$delta$se) <- rownames(Par$delta$lower) <- rownames(Par$delta$upper) <- paste0("ID:",rep(unique(m$data$ID),mixtures)) if(mixtures>1) rownames(Par$delta$est) <- rownames(Par$delta$se) <- rownames(Par$delta$lower) <- rownames(Par$delta$upper) <- paste0("ID:",rep(unique(m$data$ID),mixtures),"_mix",rep(1:mixtures,each=nbAnimals)) } return(Par) } get_gamma <- function(beta,covs,nbStates,i,j,betaRef,betaCons,workBounds=NULL,mixture=1){ tmpBeta <- rep(NA,length(betaCons)) tmpBeta[unique(c(betaCons))] <- beta beta <- w2wn(matrix(tmpBeta[betaCons],nrow(betaCons),ncol(betaCons)),workBounds) gamma <- trMatrix_rcpp(nbStates,beta[(mixture-1)*ncol(covs)+1:ncol(covs),,drop=FALSE],covs,betaRef)[,,1] gamma[i,j] } get_recharge <- function(g0theta,recovs,g0covs,recharge,hierRecharge,rechargeName="recharge",covs,workBounds=NULL,k=0){ g0 <- w2wn(g0theta[1:ncol(g0covs)],workBounds$g0)
theta <- w2wn(g0theta[-(1:ncol(g0covs))],workBounds$theta) if(inherits(covs,"hierarchical")){ recLevelNames <- names(hierRecharge) rechargeNames <- paste0("recharge",gsub("level","",recLevelNames)) colInd <- which(grepl(paste0("I((level == \"",gsub("level","",recLevelNames[which(rechargeName==rechargeNames)]),"\")"),colnames(recovs),fixed=TRUE)) g0covs <- g0covs[which(covs$level==gsub("level","",recLevelNames[which(rechargeName==rechargeNames)]))[1],,drop=FALSE]
} else {
colInd <- 1:ncol(recovs)
g0covs <- g0covs[1,,drop=FALSE]
}
g <- g0 %*% t(g0covs)
rec <- cumsum(c(g,theta[colInd]%*%t(recovs[-nrow(recovs),colInd])))

if(k) rec <- rec[k]
return(rec)
}

get_gamma_recharge <- function(beta,covs,formula,hierRecharge,nbStates,i,j,betaRef,betaCons,workBounds=NULL,mixture=1){

recharge <- expandRechargeFormulas(hierRecharge)

recovs <- stats::model.matrix(recharge$theta,covs) g0covs <- stats::model.matrix(recharge$g0,covs)

tmpBeta <- rep(NA,length(betaCons))
tmpBeta[unique(c(betaCons))] <- beta[1:(length(beta)-(ncol(recovs)))]
tmpBeta <- matrix(tmpBeta[betaCons],nrow(betaCons),ncol(betaCons))
beta <- w2wn(c(tmpBeta,beta[length(beta)-(ncol(recovs)-1):0]),workBounds)

#g0 <- beta[length(beta)-ncol(recovs)-(ncol(g0covs)-1):0] %*% t(g0covs)
theta <- beta[length(beta)-(ncol(recovs)-1):0]

if(inherits(covs,"hierarchical")){
recLevels <- length(hierRecharge)
recLevelNames <- names(hierRecharge)
rechargeNames <- paste0("recharge",gsub("level","",recLevelNames))
colInd <- lapply(recLevelNames,function(x) which(grepl(paste0("I((level == \"",gsub("level","",x),"\")"),colnames(recovs),fixed=TRUE)))
} else {
recLevels <- 1
rechargeNames <- "recharge"
colInd <- list(1:ncol(recovs))
}
for(iLevel in 1:recLevels){
covs[,rechargeNames[iLevel]] <- covs[,rechargeNames[iLevel],drop=FALSE] + theta[colInd[[iLevel]]]%*%t(recovs[,colInd[[iLevel]],drop=FALSE]) # g0  + theta%*%t(recovs)
}

newcovs <- stats::model.matrix(formula,covs)
beta <- matrix(beta[1:(length(beta)-(ncol(recovs)))],ncol=nbStates*(nbStates-1))
gamma <- trMatrix_rcpp(nbStates,beta[(mixture-1)*ncol(newcovs)+1:ncol(newcovs),,drop=FALSE],newcovs,betaRef)[,,1]
gamma[i,j]
}

get_delta <- function(delta,covsDelta,i,workBounds=matrix(c(-Inf,Inf),length(delta),2,byrow=TRUE),mixture=1){
delta <- w2wn(delta,workBounds)
nbCovsDelta <- ncol(covsDelta)-1
delta <- delta[(mixture-1)*(nbCovsDelta+1)+1:(nbCovsDelta+1),]
delta <- c(rep(0,nbCovsDelta+1),delta)
deltaXB <- covsDelta%*%matrix(delta,nrow=nbCovsDelta+1)
expdelta <- exp(deltaXB)
if(!is.finite(sum(expdelta))){
tmp <- exp(Brobdingnag::as.brob(as.vector(deltaXB)))
delta <- as.numeric(tmp/Brobdingnag::sum(tmp))
} else {
delta <- expdelta/sum(expdelta)
}
delta[i]
}

parm_list<-function(est,se,lower,upper,rnames,cnames){
Par <- list(est=est,se=se,lower=lower,upper=upper)
rownames(Par$est) <- rnames rownames(Par$se) <- rnames
rownames(Par$lower) <- rnames rownames(Par$upper) <- rnames
colnames(Par$est) <- cnames colnames(Par$se) <- cnames
colnames(Par$lower) <- cnames colnames(Par$upper) <- cnames
Par
}

get_CI<-function(wpar,Par,m,ind,DM,DMind,Bounds,Sigma,circularAngleMean,consensus,nbStates,alpha,rnames,cnames,nc,meanind,workBounds,dist){

w<-wpar[ind]
lower<-upper<-se<-rep(NA,nrow(DM))
if(!is.null(Sigma)){
for(k in 1:nrow(DM)){
if(inherits(dN,"error")){
warning(dN)
dN <- rep(NaN,length(w))
}
se[k]<-suppressWarnings(sqrt(dN%*%Sigma[ind,ind]%*%dN))
lower[k] <- Par[k] - qnorm(1-(1-alpha)/2) * se[k]
upper[k] <- Par[k] + qnorm(1-(1-alpha)/2) * se[k]
#cn<-exp(qnorm(1-(1-alpha)/2)*sqrt(log(1+(se[k]/Par[k])^2)))
#lower[k]<-Par[k]/cn
#upper[k]<-Par[k]*cn
}
}
est<-matrix(Par,ncol=nbStates,byrow=TRUE)
l<-matrix(lower,ncol=nbStates,byrow=TRUE)
u<-matrix(upper,ncol=nbStates,byrow=TRUE)
s<-matrix(se,ncol=nbStates,byrow=TRUE)
out <- parm_list(est[1:length(rnames),,drop=FALSE],s[1:length(rnames),,drop=FALSE],l[1:length(rnames),,drop=FALSE],u[1:length(rnames),,drop=FALSE],rnames,cnames)
out
}

CI_angle<-function(wpar,Par,m,ind,DM,DMind,Bounds,Sigma,circularAngleMean,consensus,nbStates,alpha,rnames,cnames,nc,meanind,workBounds,dist){

w<-wpar[ind]
lower<-upper<-se<-rep(NA,nrow(DM))
if(!is.null(Sigma)){
for(k in 1:nrow(DM)){
se[k]<-suppressWarnings(sqrt(dN%*%Sigma[ind,ind]%*%dN))
lower[k] <- Par[k] - qnorm(1-(1-alpha)/2) * se[k]
upper[k] <- Par[k] + qnorm(1-(1-alpha)/2) * se[k]
#cn<-exp(qnorm(1-(1-alpha)/2)*sqrt(log(1+(se[k]/Par[k])^2)))
#lower[k]<-Par[k]/cn
#upper[k]<-Par[k]*cn
}
}
est<-matrix(Par,ncol=nbStates,byrow=TRUE)
l<-matrix(lower,ncol=nbStates,byrow=TRUE)
u<-matrix(upper,ncol=nbStates,byrow=TRUE)
s<-matrix(se,ncol=nbStates,byrow=TRUE)
out <- parm_list(est,s,l,u,rnames,cnames)
out
}

w2nDMangle<-function(w,bounds,DM,DMind,nbObs,circularAngleMean,consensus,nbStates,k,nc,meanind,workBounds,dist){

bounds[,1] <- -Inf
bounds[which(bounds[,2]!=1),2] <- Inf

foo <- length(w) - nbStates + 1
x <- w[(foo - nbStates):(foo - 1)]
y <- w[foo:length(w)]
angleMean <- Arg(x + (0+1i) * y)
kappa <- sqrt(x^2 + y^2)
w[(foo - nbStates):(foo - 1)] <- angleMean
w[foo:length(w)] <- kappa

w2nDM(w,bounds,DM,DMind,nbObs,circularAngleMean,consensus,nbStates,k,nc,meanind,workBounds,dist)
}

#' @rdname CIreal
#' @method CIreal hierarchical
#' @importFrom utils tail
#' @importFrom Brobdingnag as.brob sum
#' @importFrom CircStats circ.mean
#' @export
CIreal.hierarchical <- function(m,alpha=0.95,covs=NULL,parms=NULL){

if(is.miSum(m) | is.miHMM(m)){
if(is.miHMM(m)) m <- m$miSum m <- formatmiSum(m) m$mod$hessian <- NA m$mod$Sigma <- matrix(0,length(m$mod$estimate),length(m$mod$estimate)) if(length(m$conditions$optInd)) m$mod$Sigma[-m$conditions$optInd,-m$conditions$optInd] <- m$MIcombine$variance else m$mod$Sigma <- m$MIcombine$variance m$CIreal <- m$Par$real
m <- momentuHMM(m)
} else if(!inherits(m,"momentuHierHMM")) stop("m must be a momentuHierHMM, hierarchical miSum, or hierarchical miHMM object")

installDataTree()

if(!is.null(covs)){
ci <- CIreal.default(m=m,alpha=alpha,covs=covs,parms=parms)
tmpcovs <- covs
tmpcovs$level <- "1" } else { ci <- m$CIreal
tmpcovs <- data.frame(level="1")
}

class(tmpcovs) <- append("hierarchical",class(tmpcovs))

CIgamma <- data.tree::Node$new("hierGamma") CIdelta <- data.tree::Node$new("hierDelta")

hierStates <- m$conditions$hierStates

mixtures <- m$conditions$mixtures

ref <- hierStates$Get(function(x) data.tree::Aggregate(x,"state",min),filterFun=function(x) x$level==2)
mixref <- rep(seq(0,mixtures*length(m$stateNames)-1,length(m$stateNames)),each=length(ref))+ref
if(mixtures>1) nameref <- paste0(rep(names(ref),mixtures),"_mix",rep(1:mixtures,each=length(ref)))
else nameref <- names(ref)

CIgamma$AddChild("level1",gamma=lapply(m$CIreal$gamma,function(x) matrix(x[mixref,ref],length(mixref),length(ref),dimnames=list(nameref,names(ref))))) CIdelta$AddChild("level1",delta=lapply(m$CIreal$delta,function(x) matrix(x[,ref],nrow(x),dimnames=list(rownames(x),names(ref)))))

for(j in 2:(hierStates$height-1)){ t <- data.tree::Traverse(hierStates,filterFun=function(x) x$level==j)
names(t) <- hierStates$Get("name",filterFun=function(x) x$level==j)

tmpcovs$level <- j tmpGamma <- CIreal.default(m=m,alpha=alpha,covs=tmpcovs,parms=c("gamma"))$gamma
tmpcovs$level <- paste0(j,"i") tmpDelta <- CIreal.default(m=m,alpha=alpha,covs=tmpcovs,parms=c("gamma"))$gamma
CIgamma$AddChild(paste0("level",j),gamma=list()) CIdelta$AddChild(paste0("level",j),delta=list())

ref <- hierStates$Get(function(x) data.tree::Aggregate(x,"state",min),filterFun=function(x) x$level==j)

for(k in names(t)){
levelStates <- t[[k]]$Get(function(x) data.tree::Aggregate(x,"state",min),filterFun=function(x) x$level==j+1)#t[[k]]$Get("state",filterFun = data.tree::isLeaf) if(!is.null(levelStates)){ mixref <- rep(seq(0,mixtures*length(m$stateNames)-1,length(m$stateNames)),each=length(levelStates))+levelStates if(mixtures>1) nameref <- paste0(rep(names(levelStates),mixtures),"_mix",rep(1:mixtures,each=length(levelStates))) else nameref <- names(levelStates) mixrefk <- rep(seq(0,mixtures*length(m$stateNames)-1,length(m$stateNames)),each=length(ref[[k]]))+ref[[k]] if(mixtures>1) namerefk <- paste0(rep(names(ref[k]),mixtures),"_mix",1:mixtures) else namerefk <- k CIgamma[[paste0("level",j)]]$gamma[[k]] <- lapply(tmpGamma,function(x) matrix(x[mixref,levelStates],length(mixref),length(levelStates),dimnames=list(nameref,names(levelStates))))
CIdelta[[paste0("level",j)]]$delta[[k]] <- lapply(tmpDelta,function(x) matrix(x[mixrefk,levelStates],mixtures,dimnames=list(namerefk,names(levelStates)))) } } } ci$hierDelta <- CIdelta
ci\$hierGamma <- CIgamma
ci
}


## Try the momentuHMM package in your browser

Any scripts or data that you put into this service are public.

momentuHMM documentation built on Sept. 5, 2021, 5:17 p.m.