Covariate transition probabilities and CIs for pooled model

[juliagulka]

I am working with some seabird satellite telemetry data, and so have been using the multiple simulation approach (MIfitHMM) to incorporate location error, and am in turn interested in the influence on different environmental covariates on transition probabilities. I am hoping to extract/predict the data needed to reproduce the covariate-specific transition probability plots produced by plot.momentuHMM (for plot customization) with both estimated transition probabilities and confidence intervals for a range of values for a specific covariate with mean values set for other covariates.

From my understanding, CIreal works to predict TPs and CIs for a singular model fit, but will not work on the pooled model. Looking at the documentation for plot.momentuHMM, I see that confidence intervals are calculated using finite-difference approximations, but is there an easy way to extract these values used to create these plots/predict TP estimates and CIs for a set of covariate values using a pooled model?

Thanks!

[bmcclintock]

Assuming the individual model fits are "good" (i.e., converged with no issues for the parameter or variance-covariance matrix estimates), then CIreal can be used for obtaining transition probability estimates and CIs from a pooled model for a given set of covariate values. CIreal does accept miSum object output from MIfitHMM, e.g., CIreal(mifit$miSum), where mifit is the object returned by MIfitHMM when poolEstimates=TRUE. The covariate values are specified using the covs argument, but CIreal requires covs to be a single row of values. To replicate the covariate values used by plot.momentuHMM, one would need to use a sequence of length 101 covering the range of the covariate of interest (e.g. seq(min(cov1),max(cov1),length=101)) and the mean value(s) for any other covariates in the model.

Here is some example code for the basic example that loads with momentuHMM (the only difference in your case is you'd be using mifit$miSum and your specific covariates instead of example$m):

covSeq <- data.frame(cov1=seq(min(example$m$data$cov1),max(example$m$data$cov1),length=101),cov2=mean(example$m$data$cov2))
estCI <- vector('list',4)
names(estCI) <- c("est","se","lower","upper")
estCI[names(estCI)] <- list(list())
for(i in 1:nrow(covSeq)){
  out <- CIreal(example$m, covs=covSeq[i,])$gamma
  estCI$est[[i]] <- out$est
  estCI$se[[i]] <- out$se
  estCI$lower[[i]] <- out$lower
  estCI$upper[[i]] <- out$upper
}
plot(example$m,plotCI=TRUE,ask=FALSE)

# plot transition probabilities and CIs for cov1
par(mfrow=c(2,2))
for(i in 1:2){
  for(j in 1:2){
    plot(covSeq$cov1,unlist(lapply(estCI$est,function(x) x[i,j])),ylim=c(0,1),type="l",ylab=paste(i,"->",j),xlab="cov1")
    lines(covSeq$cov1,unlist(lapply(estCI$lower,function(x) x[i,j])),lty=2)
    lines(covSeq$cov1,unlist(lapply(estCI$upper,function(x) x[i,j])),lty=2)
  }
}

[juliagulka]

I believe the individual model fits are "good" with no convergence issues. I have run the MIfitHMM model with poolEstimates=TRUE, however when I attempt to run the sample code that you provided, when I run the for loop to execute CIreal on each row of the covSeq dataframe using my pooled model$miSum object, I get the following error: "Error in CIreal.default(m$miSum, covs = covSeq[i, ]) : 'm' must be a momentuHMM object (as output by fitHMM)". If I run the code on a single model fit of the unpooled model (MIfitHMM with poolEstimates=FALSE), I have no problems in executing the code and creating the plots. I tried this with my own data/model outputs as well as the loggerhead turtle dataset in the package vignette and get the same error. Please let me know if I am missing something crucial to allow CIreal to work on the miSum object. Note: I am using the most up-to-date version of MomentuHMM (version 1.5.3).

Here is the example code that I ran with unpooled (MIfitHMM with poolEstimates = FALSE) and pooled (MIfitHMM with poolEstimates=TRUE) model outputs. Any additional insight would be greatly appreciated!

`#load data
unpooled <- readRDS("unpooled.rds")
pooled <- readRDS("pooled.rds")

#transition probability with single fit of model (unpooled) - this works
#the two spatial covariates in the model are chl a (logchl.scale) and sst (sst.scale)
covSeq <- data.frame(logchl.scale=seq(min(unpooled[[1]]$data$logchl.scale),max(unpooled[[1]]$data$logchl.scale),length=101),sst.scale=mean(unpooled[[1]]$data$sst.scale))
estCI <- vector('list',4)
names(estCI) <- c("est","se","lower","upper")
estCI[names(estCI)] <- list(list())

for(i in 1:nrow(covSeq)){
out <- CIreal(unpooled[[1]], covs=covSeq[i,])$gamma
estCI$est[[i]] <- out$est
estCI$se[[i]] <- out$se
estCI$lower[[i]] <- out$lower
estCI$upper[[i]] <- out$upper
}

par(mfrow=c(2,2))
for(i in 1:2){
for(j in 1:2){
plot(covSeq$logchl.scale,unlist(lapply(estCI$est,function(x) x[i,j])),ylim=c(0,1),type="l",ylab=paste(i,"->",j),xlab="chl a")
lines(covSeq$logchl.scale,unlist(lapply(estCI$lower,function(x) x[i,j])),lty=2)
lines(covSeq$logchl.scale,unlist(lapply(estCI$upper,function(x) x[i,j])),lty=2)
}
}

#running the same script on the pooled model produces an error
covSeq <- data.frame(logchl.scale=seq(min(pooled$miSum$data$logchl.scale),max(pooled$miSum$data$logchl.scale),length=101),sst.scale=mean(pooled$miSum$data$sst.scale))
estCI <- vector('list',4)
names(estCI) <- c("est","se","lower","upper")
estCI[names(estCI)] <- list(list())

for(i in 1:nrow(covSeq)){
out <- CIreal(pooled$miSum, covs=covSeq[i,])$gamma
estCI$est[[i]] <- out$est
estCI$se[[i]] <- out$se
estCI$lower[[i]] <- out$lower
estCI$upper[[i]] <- out$upper
}
#get error: Error in CIreal.default(pooled$miSum, covs = covSeq[i, ]) :
#'m' must be a momentuHMM object (as output by fitHMM)
`

[bmcclintock]

Ah, CIreal wasn't correctly checking for miSum object inputs. Thank you for bringing this to my attention. I've fixed this for the next release (v1.5.4), which will be submitted when CRAN returns from vacation next month. CIreal will also now accept miHMM objects, so no need for the $miSum part (i.e., either CIreal(mifit$miSum) or CIreal(mifit) will work). In the meantime, you can install it from here:

remotes::install_github('bmcclintock/momentuHMM')