www.phidot.org

Hello,

I am running a single-season two-species model in RPresence using PsiBA parameterization. I am primarily interested in determining the effect of an apex predator on smaller carnivores. I am little bit caught up in formalizing the codes. For species-specific occupancy without any interaction where psiA, psiBA=psiBa. the code is :

mod1<-occMod(model=list(psi~SP,p~SP),data=data,type="so.2sp.1",parameterization="PsiBA")

The matrix of the model is:

a1 a2
psiA "1" "0"
psiBA "1" "1"
psiBa "1" "1"

Is it correct for estimating species-specific occupancy without any interactions? or will it should be like this:
a1 a2
psiA 1 0
psiBA 0 1
psiBa 0 1

How should I enter my codes for getting the above matrix where the a1 parameter of psiBA=psiBa=0. Does both the matrices signify the same thing that is occupancy without any interactions? if not what is the difference between the two.

Moreover, if anyone can tell me how will I predict my beta estimates generate through RPresence for plotting purposes.

Thank You

Both design matrices are correct and will give the same real parameter estimates. The difference is that the first design matrix estimates the logit of psiA as the beta parameter, a1, and the logit of the difference between species A and B, as the beta parameter, a2. In the second design matrix, the parameter, a2 is estimating the logit of psiB (psiBA=psiBa) independently of psiA. The first design matrix is better suited for testing for a significant difference between psiA and psiB, as you only have to test if a2 is significantly different from zero. With the second design matrix, you would have to test if a2 not equal a1, where a1 and a2 both have standard errors to deal with.

To get RPresence to produce the 2nd matrix, the formula is:

Code: Select all

psi ~ -1 + SP

The "-1" tells R to remove the intercept from the model.

You can use the "predict.occMod" function to get predicted estimates based on a model and new set of covariate data. Type "help('predict.occMod')".

Thank you for the prompt response. I understood this part but what if I want to add covariates such that it has a different influence on psiBA and psiBa so that my matrix looks like this:

a1 a2 a3 a4 a5
psiA "1" "0" "0" "0" "0"
psiBA "0" "0" "1" "X" "0"
psiBa "0" "1" "0" "0" "X"

where X is a covariate. Most of the example I saw has only shown the effect of covariate where psiBA(X)=psiBa(X). For different parameters, how to specify covariates for each one of them. There is very little information available online for such analysis in RPresence. I tried different combinations but didn't get the kind of matrix I wanted. Can you please help with how can I generate the above kind of matrix with covariates?

Also, I tried "predict.occMod" but in help file also it is showing "predict" only and I am using the code:

predict(mod2, data$unitcov, param = "psiBA")

But I am not able to define the parameter for which I want to predict real estimates. I tried many combinations but its not working.

So, this covariate affects species B differently when species A is present versus when A is absent, AND it has no effect on species A? Seems strange to me, but here goes...

First, you'll need to create the X covariate a little differently. Normally, RPresence will replicate the site covariates 3 times, once for each parameter (psiA, psiBA, psiBa). In this case, we need a site covariate which contains "0" for psiA, and the X-values for psiBA and psiBa.

Code: Select all

#  get length of detection data (nsites)
 nu=nrow(x1$det.data)
  pao=createPao(data=x1$det.data, frq=x1$frq,
        #  create site covariate where first 'nu' rows contain zero,
        #  and rest of rows contain X (I used a random value in this example)
        unitcov=data.frame(
             X=c(rep(0,nu),round(runif(2*nu)*5,0))
     ))

Now, we can call occMod using the built-in covariates, SP and INT
as well as the covariate, X.

Code: Select all

 EV2SPMOD<-occMod(model=list(psi~-1+SP+INT+X+INT:X,p~SP),
                   data=pao,type="so.2sp.1",outfile='tmp1.out')
  print(EV2SPMOD$dmat$psi)


To get predicted values for the 2-species model, we'll have to create a data-frame
containing all of the covariates, SP, INT and X.

Code: Select all

  newdata=data.frame(
    SP=factor(c(rep('A',nu),rep('B',2*nu))),
    INT=c(rep(0,2*nu),rep(1,nu)),
    X=pao$unitcov$X)

  z = predict(object=EV2SPMOD, newdata=newdata, param='psi')
print(unique(z))

Yes, for example, both the species use a certain patch or water resource, I am particularly interested in how the subordinate species respond to that resource in the presence or absence of the dominant species. The covariate might or might not have an effect on the dominant species but I am interested in how species B, in particular, responds to covariates in the absence or presence of other dominant species.

I am able to run the models based on the codes provided and got the expected results, however, I am having little trouble in making the plots for the predicted models. I want to plot the species effect on covariate X when species A is present and when it is absent.

Is there any way to plot that kind of conditional plots?

Thank You

OK, even if you're not interested in whether the covariate affects species A, you should include a model which has that effect if it's reasonably possible. If there is no effect on species A, the beta estimate will simply be near zero.

Here is code to plot the psi values:

Code: Select all

 # z[1:256, ] = site-specific estimates of psiA
  # z[257:512, ] = site-specific estimates of psiBA
  # z[512:792, ] = site-specific estimates of psiBa
 
  #  add parameter name to z data-frame
  z$param = c(rep('psiA',256),rep('psiBA',256),rep('psiBa',256))
 
  #  plot psi parameters by site
  library(ggplot2)
  ggplot(z,aes(x=rep(1:256,3),y=est)) + xlab('site') +
    geom_point(aes(color=param))

Okay, I will include the covariate effect on species A and will check the models. I want to plot predicted psi values against the covariate and not for each site. The prediction for earlier code also estimates values for just n sites, I want to predict psi estimate for each value of covariate separated by some difference and want to plot that predicted value against the covariate (as done in unmarked package). The code helped to plot psi values against each site.

I added the covariate in the data frame however, for prediction it's only showing est values with SE. I want to include both predicted values and covariate values in a single data frame so as to plot those predicted values against covariates for different parameters.

Thank You

If I understand what you're saying, you want to plot the psi parameters over a range of new covariate values, right? If so, just create a data-frame specifying the new covariate values and call the predict function using the new data-frame. Here is some sample code:

Code: Select all

  EV2SPMOD<-occMod(model=list(psi~-1+SP+INT+X+INT:X,p~SP),
                   data=pao,type="so.2sp.1",outfile='tmp1.out')
  print(EV2SPMOD$dmat$psi)

  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)),
    X=newcov)

  z = predict(object=EV2SPMOD, newdata=newdata, param='psi')

  #  add parameter name to z data-frame
  z$param = c(rep('psiA',nu),rep('psiBA',nu),rep('psiBa',nu))

  #  plot psi parameters by site
  library(ggplot2)
  ggplot(z,aes(x=rep(1:nu,3),y=est)) + xlab('site') +
    geom_line(aes(color=param))

I realized I plotted site vs psi estimates in last code posted. Here is code to plot the covariate versus the psi's:

Code: Select all

  EV2SPMOD<-occMod(model=list(psi~-1+SP+INT+X+INT:X,p~SP),
                   data=pao,type="so.2sp.1",outfile='tmp1.out')
  print(EV2SPMOD$dmat$psi)

  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)),
    X=newcov)

  z = predict(object=EV2SPMOD, newdata=newdata, param='psi')

  #  add parameter name to z data-frame
  z$param = c(rep('psiA',nu),rep('psiBA',nu),rep('psiBa',nu))
  z$X = newdata$X

  #  plot psi parameters by site
  library(ggplot2)
  ggplot(z,aes(x=X,y=est)) + xlab('covariate X') +
    geom_line(aes(color=param))


Thank you so much for the help. I was able to make the expected plot with the effect of the covariate on different parameters. Running the models in RPresence and specifying the parameter conditions is a little different than other packages.

While specifying the covariate values in the data-frame, the length of the data-frame was differing from the earlier data frame that was made to run models where psiA was kept "0", so it was showing error while predicting the model. Then I combined the predicted model (having newdata) with the new data frame and I was able to get the expected plot.

Thank you again.