Error messages: 'Covariate' not found and more

[Kalatzin]

Hi,
I have been running occupancy models for my masters thesis, originally I was working in Program PRESENCE applying a 2-step approach to 2-species occupancy modeling. After finishing my analysis this week I realized I had to make a major revision to my occurrence data after I detected an error (introduced by the pivot table I used to compile the occurrences). This error meant I had to redo my whole analysis (that I've been working on for days if not weeks) after fixing the error in the data. However, the thought of manually re-inputting the models for all 6 analyses (local-scale & landscape-scale) sounded way too daunting in Program PRESENCE as I had (collectively) ran ~200+ models between all 6 analyses (single-species mods for both species and 2-species models times 2 for local-scale covs. and landscape scale covs.) so I thought about trying to re-run my analyses in RPresence to increase efficiency in inputting the models.

After finding the package documentation and writing some script, I got the constant models to run fine as well as the models with a covariate effect on detection. I also had success running 2 models with site covariate effects on occupancy. However, when it came to trying to run any of the models with sample covariates (survcovs), same ones I used in models with cov. effects on detection, I got the error message:

'Error in eval(predvars, data, env) : object 'temp' not found'

This happened for each covariate, I then tried adding 'cov.list= wtoad$survcovs' into the occmod() code to see if that would help and got the same error. Then I tried:

mod9.toad<-occMod(model=list(psi~survcov$temp,p~1),cov.list = wtoad$survcov, data= wtoad, type="so", conf=0.95, outfile='mod9.toad.out', modfitboot=500)

and instead got this error:
Error in `.rowNamesDF<-`(x, value = value) : invalid 'row.names' length
In addition: Warning messages:
1: In getbetas(v[i], psi.dm) : NAs introduced by coercion
2: In sqrt(var) : NaNs produced

I tried pulling the covariates out as a matrix and data frame and calling the object into the model via cov.list= with no luck as a matrix or data frame, and got the same error messages from before. I tried everything I could think of, searched for possible solutions etc. to no avail.

At this point, I'm out of time to try and troubleshoot this and make this work as my thesis is due in a couple weeks. So I decided to abandon this approach that I hoped would be a solution and am back to rerunning all 6 analyses in Program PRESENCE. I really hope I can get this done on time to submit my draft to my committee ASAP. I am posting to see if anyone had any insight as to why this happened and what could have been done if there's a next time.

Lastly, I saw in another forum post that the latest version of RPresence is 4.13.51, however, no matter how many times I installed the package it downloaded v. 4.13.49, has anyone else had an issue with this? The documentation for RPresence is so sparse, I've had no luck tracking down which version of R is needed for the package to run or download the latest one. I even updated R at the beginning since I had an older version.

If you're reading this and got this far, thank you for taking the time. Wishing you a better week than I've been having. <3

[jhines]

Hi, I'm sorry you're having trouble just before your thesis is due. Without seeing the code, I can't tell whats wrong, but one thing I see in the code you posted is it appears you are trying to model occupancy (psi) as a function of a survey covariate. Since occupancy is constant across surveys, it cannot be modelled as a function of a survey covariate. If you'd like to send me your code and data, I'd be happy to diagnose the error.

Jim (jhines@usgs.gov)

[Kalatzin]

Hi Jim,
Thank you for taking the time to reply, and for the consolation.

Yeah I can share the code. Also, I didn't catch that I wrote the models as constant, so to add the intercept it should be as written as: psi~1+temp,p~1 , correct? I edited the code to reflect that but it still gives me the same error message: 'Error in eval(predvars, data, env) : object 'temp' not found'

Code: Select all

# Bring in the western toad PRESENCE files
wtoad= readPao('wtoad_local_analysis.pao')
## fit some models
mod1.toad<-occMod(model=list(psi~1,p~1),data=wtoad,type="so", conf=0.95, outfile='mod1.toad.out', modfitboot=500)
mod2.toad<-occMod(model=list(psi~1,p~1+SURVEY),data=wtoad,type="so", conf=0.95, outfile='mod2.toad.out', modfitboot=500)
mod3.toad<-occMod(model=list(psi~1,p~1+SampleAge),p.cov=data$survcov, data=wtoad,type="so", conf=0.95, outfile='mod3.toad.out', modfitboot=500)
mod4.toad<-occMod(model=list(psi~1,p~1+temp), p.cov= wtoad$survcov, data=wtoad,type="so" , conf=0.95, outfile='mod4.toad.out', modfitboot=500)
mod5.toad<-occMod(model=list(psi~1,p~1+flow), p.cov= wtoad$survcov,data=wtoad,type="so", conf=0.95, outfile='mod5.toad.out', modfitboot=500)
mod6.toad<-occMod(model=list(psi~1,p~1+depth),p.cov=wtoad$survcov,data=wtoad,type="so", conf=0.95, outfile='mod6.toad.out', modfitboot=500)
mod7.toad<-occMod(model=list(psi~1+em.veg,p~1),psi.cov= wtoad$unitcov,data=wtoad,type="so", conf=0.95, outfile='mod7.toad.out', modfitboot=500)
mod8.toad<-occMod(model=list(psi~1+tree.veg,p~1),psi.cov= wtoad$unitvcov,data=wtoad,type="so", conf=0.95, outfile='mod8.toad.out', modfitboot=500)
  ## the code below is where R gives me the error I mentioned in my post ##
mod9.toad<-occMod(model=list(psi~1+temp,p~1),psi.cov = wtoad$survcov,data= wtoad,type="so", conf=0.95, outfile='mod9.toad.out', modfitboot=500)
mod10.toad<-occMod(model=list(psi~1+flow,p~1),psi.cov= wtoad$survcov, data=wtoad,type="so", conf=0.95, outfile='mod10.toad.out', modfitboot=500)
mod11.toad<-occMod(model=list(psi~depth,p~1),psi.cov=  wtoad$survcov, data=wtoad,type="so", conf=0.95, outfile='mod11.toad.out', modfitboot=500)
mod12.toad<-occMod(model=list(psi~1+tempsqr,p~1),psi.cov= wtoad$survcov, data=wtoad,type="so", conf=0.95, outfile='mod12.toad.out', modfitboot=500)



I also originally tried the code as

Code: Select all

mod9.toad<-occMod(model=list(psi~1+temp,p~1), data= wtoad,type="so", conf=0.95, outfile='mod9.toad.out', modfitboot=500)

and then as:

Code: Select all

mod9.toad<-occMod(model=list(psi~1+temp,p~1),cov.list = wtoad$survcov, data= wtoad,type="so", conf=0.95, outfile='mod9.toad.out', modfitboot=500)

before I used the psi.cov= within the occMod() function, but it gives me the same error when I try to use any of the survey covariates (temp, flow, depth, tempsqr) which are all continuous variables, however the unitcovs ran just fine and those are binary, discrete variables.

If I can get the code to run the models I need, I could theoretically get through the next part of my analysis much faster, the landscape-scale analyses. I've barely managed to get to the second step of the modeling process for the local-scale analyses. The code I provided was for running part of the first-step (single-season, single species occu models),and hope it would work smoothly for the 2-species occupancy modeling too.

[jhines]

There is no need to add the intercept term in the formula when you have covariates as R will add it automatically. So, the formula would simply be: psi~temp. I assume that temp is a "site" covariate (ie., does not change from survey to survey). If not, temp can be used to model detection (p), but not occupancy. Actually, it might actually run, but it would only use the temp data from the first survey. If you really want temperature to be an influencing factor for occupancy, you could compute an average temperature over all surveys for each site and make that a site covariate.

If you have the covariates incorporated in the pao object, there is no need to include the additional arguments in occMod for them (p.cov, psi.cov). If you did have covariates which are not in the pao object, the way to add them to the occMod function is:

Code: Select all

modx.toad <- occMod(model=list(psi~dist, p~temp),
                    cov.list=list(psi.cov=new_site_cov_dataframe,
                                   p.cov=new_survey_cov_dataframe)
)

You don't need to do a goodness-of-fit test for every model, only the most general model in the model-set (usually the one with the most parameters). If that model fits, then the AIC table is OK to use to select the most parsimonious model. When you do a goodness-of-fit test, it is recommended to do many bootstraps (eg., 10,000) since the detection-history expected values are so small.

If you would send me the pao file off-list, I will run those models to make sure they work.

[Kalatzin]

Hi,
Thank you for clarifying that. My original models didn't have psi~1+ covariate, it was written as psi~covariate, etc. from what I had seen in previous posts.

The covariates temp, depth, flow, and tempsqr are not site covariates, they are sample covariates that change across sites and visits.

Thank you for clarifying about the bootstrapping too! I was running individual covariates first so I can narrow it down to 2-3 relevant covariates (out of the 6 covariates I have) for a global model because my sample size is only 30 sites and I didn't want to fit a global model with more than 3 covariates.

Okay I'll send the .pao file.
I really appreciate your help on this! It feels like I was overboard and someone threw me a life-saver.

[jhines]

One thing I forgot to mention is that you can look at the design matrix in the output object.

print(mod4.frog$dmat)
$psi
a1
psi "1"

$p
b1 b2
p1 "1" "p.Temp"
p2 "1" "p.Temp"
p3 "1" "p.Temp"
p4 "1" "p.Temp"
p5 "1" "p.Temp"
p6 "1" "p.Temp"
p7 "1" "p.Temp"
p8 "1" "p.Temp"

You can see that the design matrix for p has an intercept column (all 1's) and a column for the effect of temp on detection (p.Temp). This was produced with the formula, p~Temp. RPresence is smart enough to not add an extra intercept if you use the formula, p~1+Temp. That formula produces the same design matrix. The design matrix is a great way to make sure you know how the model is structured, based on the formula given.