Plotting model-averaged SIF prediction with CI

[Hsin Cheng]

Hi,

I am building single-season 2 species occupancy models, and I try to predict psi and SIF against a target covariate using the modAvg() output. To see whether the trends in psi and SIF are significant, I need to plot the predicted values and confidence intervals.

I have successfully done the plot for psi with predicted values and CIs. Here's the link to the plot:
https://drive.google.com/file/d/1wProvy ... sp=sharing
The X-axis is the target covariate I hypothesized to influence the species' interaction, and the y-axis is the psi value.

However, as the right-hand side plot shows, I can only plot the predicted values of model-averaged SIF, but I can't find a way to plot its predicted CI.
Below is the code:

Code: Select all

newcov = seq(-5,5,.1)
nu=length(newcov)

newdata = data.frame(
  SP=factor(c(rep('A',nu),rep('B',2*nu))),
  INT=c(rep(0,2*nu),rep(1,nu)),
  logDog=0,
  rdist=newcov, # my target covariate
  elevation=0,
  YearPerc=0,
  forest=0,
  slope=0)

avgpsi <- modAvg(results, param = 'psi', index = 1:nrow(results$table),predict = T, newdata = newdata)
avgpsi$param = c(rep('psiA',nu),rep('psiBA',nu),rep('psiBa',nu))
avgpsi$rdist = newdata$rdist

#  Plot psi by covariate
library(ggplot2)
ggplot(avgpsi, aes(x=rdist,y=est)) + xlab('Scaled Distance to Road') +
  geom_line(aes(color=param),size=0.8)+
  geom_ribbon(aes(x =rdist, ymin=lower_0.95, ymax=upper_0.95, group=param, color=param), alpha=0.2)


# Calculate SIF
psiBA <- avgpsi[avgpsi$param == 'psiBA',]
psiBa <- avgpsi[avgpsi$param == 'psiBa',]
psiA <- avgpsi[avgpsi$param == 'psiA',]
SIF = (psiA$est*psiBA$est) / (psiA$est * (psiA$est*psiBA$est + (1 - psiA$est) * psiBa$est))
sif = data.frame(SIF, rdist = psiBA$rdist)

# Plot SIF by covariate
ggplot(sif,aes(x=rdist,y=SIF)) + xlab('Scaled Distance to Road') +
  geom_line()


If there's any solution to plotting the predicted CI of model-averaged SIF, I would be grateful if you could provide the code or any guidance.

Hsin Cheng

[jhines]

Hi,

You can compute model-averaged phi's with the command:

phi <- modAvg(aic.tab=x, param='---phi--', parmgrp='derived')

You'll need to update RPresence to the latest version (2.x.50). This should return the estimates, standard errors and 95% confidence intervals.

Jim

[Hsin Cheng]

Thank you.

This works, but it only shows the SIF of each site, and only work when the argument "predict = F". I would like to plot SIF against a single target covariate, and also get predicted model averaged values.

However, when I try to predict SIF through the newdata I build (containing all other covariates equal to 0, only the target covariate varies between -5 and 5), it generates an error.

Code: Select all

newcov = seq(-5,5,.1)
nu=length(newcov)

newdata = data.frame(
   SP=factor(c(rep('A',nu),rep('B',2*nu))),
   INT=c(rep(0,2*nu),rep(1,nu)),
   rdist=newcov,  # target covariate
   elevation=0,
   YearPerc=0,
   forest=0,
   slope=0)

phi <- modAvg (results, param='phi', parmgrp='derived', predict = T, newdata = newdata)


Error in get(param, occmod$beta) : object 'phi' not found
I've tried "---phi--" but it still doesn't work.