occuMS and predict

[Benedikt Schmidt]

Dear all,

How does one compute the confidence intervals using predict() and occuMS()?

 

I fitted a model using occuMS() and then used predict() to calculaze predicted values and their confidence intervals on the probability scale. The problem is that some upper limits of the CI were greater than 1. What’s wrong?

 

The code is shown below.

 

Thanks for your help.

 

Best wishes, Beni

 

###

library(unmarked)

 

y.data = matrix(
c(1,0,NA,NA,
0,2,2,2,
1,0,0,0,
1,1,1,0,
0,1,1,1,
0,1,1,1,
1,2,0,0,
1,1,2,0,
1,1,2,0,
0,0,0,0,
1,1,2,0,
1,2,2,0,
1,1,2,NA,
0,1,0,2,
1,2,0,0,
1,1,1,0,
0,0,0,0,
1,1,1,0,
1,1,2,0,
0,1,2,0,
1,1,2,0,
1,1,1,0,
1,1,2,0,
1,1,2,NA,
1,1,0,0,
1,1,2,0,
1,1,1,0,
0,0,0,0,
0,0,0,0,
1,1,2,2,
1,1,2,2,
1,2,2,2,
1,2,2,2,
1,2,0,0,
0,0,0,NA,
0,0,0,0,
1,0,0,0,
1,1,1,0),
nrow=38, ncol=4, byrow=T)

 

area = c(2.16743982, 3.22957599,-0.7722639,0.6139139,-0.88246398,-0.38416787,-0.90461981,-0.24423656,1.22692548,-0.43040444,-0.66019999,-0.85714662,-0.25446389,-0.92644011,-0.59230537,0.61632236,-0.08612846,-0.02361069,0.22154757,-0.64852993,-1.03235202,-0.74222298,-0.87926839,0.25986109,-0.32464664,-0.46302077,-0.72726601,2.74550461,-0.01893848,0.28137424,0.19222363,-0.04268901,1.41030859,0.89018083,-0.24518774,-0.09871367,-0.35404597, -0.20002436)

 

umf = unmarkedFrameOccuMS(y = y.data, siteCovs = data.frame(area = area))

 

detformulas <- c('~1','~1','~1')
stateformulas <- c('~ area','~1')

 

(m1 = occuMS(detformulas, stateformulas, data = umf, parameterization='condbinom'))

 

newdata = data.frame(area = area)

 

(p1 = predict(m1, newdata=newdata, type='psi'))

 

max(p1$psi$upper) > 1
print(max(p1$psi$upper))

 

[Ken Kellner]

Hi Beni,

You're not doing anything wrong. For the conditional binomial parameterization of this model unmarked calculates SEs for predicted probabilities using the delta method (see the last section of this post). Unfortunately the delta method does not take into account the 0-1 bounds of a probability. So if you have a predicted probability close to 0 or 1 and a fairly large SE, which you do in your example, you can exceed the boundaries.

The alternative is to bootstrap the confidence intervals using the model output. This is what the other (multinomial) parameterization does. See example code below. Possibly this should just be the default approach for both parameterizations.

#Number of bootstrap replicates

nsim <- 100

#Coefficient vector

Mu <- coef(m1)

#Variance-covariance matrix

Sigma <- vcov(m1)

#Sample new vectors of coefficients from multivariate random normal

#each row is a sample

samps <- MASS::mvrnorm(nsim, Mu, Sigma)

#Copy model

m1_star <- m1

#Matrices to hold output

out_psi <- out_R <- matrix(NA, nrow=numSites(getData(m1)), ncol=nsim)

for (i in 1:nsim){

  #Replace coefficient estimates with bootstrap sample

  m1_star@estimates@estimates$state@estimates <- samps[i,1:3]

  m1_star@estimates@estimates$det@estimates <- samps[i,4:6]

  

  #Run predict with modified model

  pr <- predict(m1_star, "psi", se.fit=FALSE)

  

  #Save to output object

  out_psi[,i] <- pr$psi$Predicted

  out_R[,i] <- pr$R$Predicted

}

#calculate estimate and CI

est <- predict(m1, "psi", se.fit=FALSE)

pr_psi <- data.frame(Predicted=est$psi$Predicted,

                     SE=apply(out_psi, 1, sd),

                     lower=apply(out_psi, 1, quantile, 0.025),

                     upper=apply(out_psi, 1, quantile, 0.975))

pr_R <- data.frame(Predicted=est$R$Predicted,

                     SE=apply(out_R, 1, sd),

                     lower=apply(out_R, 1, quantile, 0.025),

                     upper=apply(out_R, 1, quantile, 0.975))

#Replicate unmarked output

pr_final <- list(pr_psi, pr_R)

[Benedikt Schmidt]

Hi Ken,

Thank you very much for this helpful answer!

Best wishes, Beni

 

 

 

Gesendet: Mittwoch, 21. Oktober 2020 um 22:57 Uhr
Von: "Ken Kellner" <ken.k...@gmail.com>
An: "unmarked" <unma...@googlegroups.com>
Betreff: [unmarked] Re: occuMS and predict

--
You received this message because you are subscribed to the Google Groups "unmarked" group.
To unsubscribe from this group and stop receiving emails from it, send an email to unmarked+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/869183ab-98ec-4aea-a480-c73e0505439an%40googlegroups.com.

[Ken Kellner]

I wrote a patch to improve calculation of the confidence intervals so that they would be properly bounded, without needing to bootstrap (attached).

Just source it in after loading unmarked and predict() will be updated. I'll get the fix into the next version of the package.

Ken

To unsubscribe from this group and stop receiving emails from it, send an email to unma...@googlegroups.com.

[Ken Kellner]

I updated the patch so it would work properly with the CRAN version of unmarked. CI calculation is also now fixed in the dev version and will be in next CRAN release. Thanks Beni and Evan for help with this.

Ken