Confidence intervals for parameter estimates and mean p/psi

[Maresa]

Hello you all, 

I am getting a little confused about the different ways to get confidence intervals for detection probability and occupancy and the cov-coefficients.

 

I want two things: 

1) get an estimate for (mean?) occupancy/detection probablity with CIs for this estimate based on my top model. I would use predict() on my topmodel with the actual covariate values belonging to the covs included in the topmodel (not some made-up new data frame). Then I would just take the mean predicted psi/p + mean lower/upper CI values. If this is the way to got - do I need to backtransform anything or is predict doing this automatically? Or when would I use the backtransform() functionality?

2) get confidence intervals for the influence of every covariate of my top model on p and psi -> i would use confint(), but I am not really sure about it. 

I am reading a loooooot of different tutorials and everyone seems to do this differently. 

Then there is also ranef() for POA, but I won't even start to get into that...

Thanks in advance!

Maresa

[Josh Twining]

Hi Maresa,

For 1) predict is a really handy function - it will provide you the mean, and upper and lower confidence intervals for each set of values you provide it in newdata (what you want the mean value of isn't clear to me, your sampled data? the landscape you sampled? below I show an example of getting marginal occupancy values for sitecov1, but if you just want mean occupancy vs det of the data you modelled, you could do similar by providing your data as the newdata). 

#create single season occupancy model

top <- occu(formula = ~ scale(observationcov1) + scale(observationcov2 ) 
             ~ scale(sitecov1) +scale(sitecov2) +scale(sitecov3), data = umf)

#create newdata to predict out to 

range(siteCovs$sitecov1) #0 - 0.66

sitecov1_newdata <- data.frame(sitecov1= seq(0, 0.66, by = 0.0025),

                          sitecov2 = mean(umf@siteCovs$sitecov2),

                           sitecov3 = mean(umf@siteCovs$sitecov3))

#predict out to new data using top model

sitecov1_preds <- predict(top, newdata = sitecov1_newdata, type = "state", level = 0.85)  

#use type = "det" for detection probabilities

#you can set your confidence intervals using "level = ...."  i.e. 0.85 = 85% confidence intervals

#combine predictions (which within it has mean, lower, and upper confidence intervals - these are 95% CI by default, or whatever you set level too)

sitecov1_preds_df <- cbind(sitecov1_newdata, sitecov1_preds)

#plot it
ggplot(sitecov1_preds_df, aes(x = scale(sitecov1), y = Predicted)) + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.1, linetype = "dashed")+ geom_path(size=1)+  theme_classic()

For 2)  If you are after beta coefficients and CIs for covariates in an unmarked model you can use the functions coef and confint as you suggest - see below: 

#create single season occupancy model

top <- occu(formula = ~ scale(observationcov1) + scale(observationcov2 ) 
             ~ scale(sitecov1) +scale(sitecov2) +scale(sitecov3), data = umf)

#get beta coefficients of top model 

coef(top)

#get confidence intervals of top model

#Note you run once for state covariates and once for detection covariates.

confint(type = "state", level = 0.85, top) #level lets you set CI, which by default are 95% CI)
confint(type = "det", level = 0.85, top),

I hope this helps and intrigued to see other approaches! 

All the best

Josh

--
*** Three hierarchical modeling email lists ***
(1) unmarked (this list): for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology: for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
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/4403e292-2b5b-4af7-ae93-f5466bd0d651n%40googlegroups.com.

--

Postdoctoral Research Scientist,

321 Fernow Hall,

Cornell University,

Ithaca, New York, 

14850

Email: joshuap...@gmail.com / jpt...@cornell.edu 

Mob: (607) 262 2918

[Maresa]

Hi Josh, 

Thanks a lot, for 2) Great, so I was right there. 

For 1) Well, I wanted to get something like one overall value for an "estimated rate of detection" and "estimated occupancy" + associated CIs based on my top model and the landscape I sampled (as provided e.g. in Kupferman 2021). What you are doing is to predict along a sequence of one covariate while holding the rest at their mean to visualize the influence of this covariate on psi, right?

Well, my idea was to provide simply my sampled covs contained in my topmodel as "newhab", predict on this "old" landscape (both for psi an p) and compute the mean, like this: 

#create new data frame with the values of the covs contained in the topmodel
newhab <- data.frame(centr_d = umf@siteCovs$centr_d,
                     FRAC_patch = umf@siteCovs  $FRAC_patch,
                     LAI=  umf@siteCovs$LAI
                     )

psi_predict <- predict(m, type = "state", newdata = newhab, appendData = T, 1000)

#get estimate for occupancy 
summary( psi_predict   ) #confidence intervals of mean predicted psi
apply( psi_predict  , 2,  mean)

same for p...

Best wishes

Maresa

[Josh Twining]

Hi Maresa,

An alternate slightly quicker way to get a single mean value (and confidence intervals) from your modelled data would be to create a dataframe with a single row - all covariates at their mean values and predict out to that:

#create newdata to predict out to 

mean_pred  <- data.frame(sitecov1= mean(umf@siteCovs$sitecov1),

                          sitecov2 = mean(umf@siteCovs$sitecov2),

                           sitecov3 = mean(umf@siteCovs$sitecov3))

psi_predict <- predict(top, type = "state", newdata = mean_pred)

All the best

Josh

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/59537b49-1e12-40a9-bfdc-c909d9811733n%40googlegroups.com.

[Maresa]

Thanks, it's always good to save some lines in your code ;) 

So, I will not have to deal with backtransformation, right? seems to me, that predict() does this internally?

All the best

Maresa

[Coral Gardner]

Hi, did you get an answer to your question,  "So, I will not have to deal with backtransformation, right? seems to me, that predict() does this internally"?

I have the same question, so if you figured this part out, do you mind sharing?  I am not sure if I need to backtransform the predict() values.

My covariates are were standardized with scale() and then I used a negative binomial for my model.

Coral

[Marc Kery]

Dear Coral,

 

you don’t need to backtransform the predictions: predict() in unmarked yields predicted values on the natural, rather than the link, scale. Attached is an example with a static occupancy model from Chapter 10 in the AHM1 book; see also https://github.com/mikemeredith/AHM_code

 

But note that you will still want to backtransform the covariate values, since nobody is interested in a plot with, say, scaled rainfall on the x axis.

 

Best regards  --- Marc

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/b900be02-a0ee-4936-9211-4312a0d9942an%40googlegroups.com.

[Ken Kellner]

To add to Marc's answer, if you used scale() on the covariates before you put them in the unmarkedFrame, you should make sure they are also scaled the same way when provided to predict() in the newdata. And then as Marc says when plotting you should convert the covariates back to the original scale so they are easier to interpret.

However if you used scale() inside the formulas provided to unmarked you don't need to do any extra steps. In this case unmarked will scale the covariates for you in all steps including automatically with predict().

Ken

[Marc Kery]

Dear Ken,

 

thanks for the important addition. And sometimes, when we have some time, we ought perhaps to collect such tipps on some website for unmarked 😊

To view this discussion on the web visit https://groups.google.com/d/msgid/unmarked/b2ae8f8e-6beb-492f-9206-22d81110e7c0n%40googlegroups.com.