Pcount function - Missing estimates in global model

[Whitney Wiest]

Hi Unmarked Users,

I'm pretty new to the package. I'm using the pcount and occu
functions to determine which of my survey covariates most strongly
influence detection probablility for salt marsh bird data. Each time
I run the global model for either function I am missing estimates for
the same factor levels of different variables.

Here I thought noise0 was the intercept, however, I'm wondering if
sky0, wind0, and date1 are all rolled up in the intercept estimate?
All covariates are factors, except for temp (numeric).

> (clra.abun.8global <- pcount(~noise+sky+temp+wind+date ~1, clra.abun))

Call:
pcount(formula = ~noise + sky + temp + wind + date ~ 1, data =
clra.abun)

Abundance:
Estimate SE z P(>|z|)
1.23 0.214 5.76 8.64e-09

Detection:
Estimate SE z P(>|z|)
(Intercept) 0.1153 0.8264 0.139 8.89e-01
noise1 0.4087 0.2331 1.753 7.96e-02
noise2 0.4752 0.2941 1.616 1.06e-01
noise3 0.3594 1.1556 0.311 7.56e-01
sky1 0.2026 0.2300 0.881 3.78e-01
sky2 0.4068 0.2019 2.014 4.40e-02
sky5 -0.4254 0.8318 -0.511 6.09e-01
temp -0.0656 0.0298 -2.204 2.75e-02
wind1 -0.3391 0.3168 -1.070 2.84e-01
wind2 -0.6424 0.3099 -2.073 3.82e-02
wind3 -0.5942 0.3435 -1.730 8.36e-02
wind4 -0.4546 0.4458 -1.020 3.08e-01
date2 -0.9081 0.1833 -4.955 7.25e-07
date3 -1.8589 0.4584 -4.055 5.01e-05


A separate issue: For this model: occu(~temp~1, clra.occu)), I
continue to receive the warning message
In sqrt(diag(vcov(obj))) : NaNs produced. There are no missing
values in the data.

In the results:
Occupancy:
Estimate SE z P(>|z|)
16 NaN NaN NaN

Any advice/guidance is greatly appreciated. Thanks for your time!

Whitney

[Matthew Giovanni]

Hi Whitney.  I think you are correct in assuming, for your categorical covariates, that noise 0, sky0, wind0, and date1 are all represented by the model's intercept.  So if you calculated estimates of pdet for sky 1, 2, and 5, they would all be relative to sky0, but not to one another.  If you wanted to change the reference level for the intercept to sky 5, for example, so you could compere the influence of sky 5 with the other levels, then you could use something like:

data=read.csv("your data",header=TRUE)
data$sky=relevel(data$sky,ref="sky5")

This is one way to generate something akin to pairwise comparisons from your model.

Richard, Andy, and Ian could diagnose your warning far more effectively than me, but I usually get that when I have missing values or not enough data to converge on an mle parameter estimate....I think.  Sometimes centering or z-standardizing the data for your continuous covariates before fitting models will help w convergence issues.  I assume "temp" represents temperature or some other continuous variable, so you might try standardizing it before fitting the occu model with something like:

umf=csvToUMF("data.csv",long=FALSE,type="unmarkedFramePCount")
sc=siteCovs(umf)
sc[,5:10]=scale(sc[,5:10])
siteCovs(umf)=sc
Null=pcount(~1 ~1+offset(log(AREA)),umf)

where "5" represents the first column with some continuous covariate (like "temp") following your "site" column (1) and columns 2-4 with your repeated-visit count data (in this case, 3 replicate visits), while "10" would be your last covariate data column preceding your pdet covariates.  This is all assuming your data is in "wide" rather than "long" format.  I hope this helps!

Matt

[Richard Chandler]

Hi Whitney,

I would just add that these questions are more related to R and linear models than to unmarked. The way that categorical variables are converted to dummy variables is controlled by R's "contrasts" option. The default is "contr.treatment" -- you could set it to "contr.sum" using:

help(contrasts)
options(contrasts=c("contr.sum", "contr.poly"))

If you want estimates for each level of your factor, you can use the predict method in unmarked. See help(unmarked ) for an example.

Or, you can remove the intercept from the model formula using something like

~sky-1.

Richard

_____________________________________
Richard Chandler, post-doc
USGS Patuxent Wildlife Research Center
301-497-5696

From:	Matthew Giovanni <matthew...@gmail.com>
To:	unma...@googlegroups.com
Date:	09/25/2011 09:45 PM
Subject:	Re: [unmarked] Pcount function - Missing estimates in global model
Sent by:	unma...@googlegroups.com

---

[Whitney Wiest]

Hi Richard and Matt,

Thanks for your help a few weeks ago!

UCLA has a good page on contrasts, which helped clear up the confusion
surrounding the intercept estimates : http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm

Also, my data is in "long" format, but I was able to figure out how to
scale temp from your code, which solved my Nan problem in both the
temp and global models.

I'm free and clear of Nans, except for one null model. Since this is
happening with the null model and not a covariate, I'm not sure how to
go about fixing this. Back when I was using Presence this model gave
me trouble as well and I eventually got it to work by only using data
from 6 surveys vs. 8 surveys (data from 7 surveys wouldn't work
either).

> sesp.occu=unmarkedFrameOccu(y=sesp[,4:11], obsCovs = list(noise = sesp[,12:19], temp = sesp[,20:27], sky = sesp[,28:35], wind = sesp[,36:43], date = sesp[,44:51]))
>
> (sesp.occu.8null <- occu(~1~1, sesp.occu))

Call:
occu(formula = ~1 ~ 1, data = sesp.occu)

Occupancy:
Estimate SE z P(>|z|)

23.8 NaN NaN NaN

Detection:
Estimate SE z P(>|z|)

1.76 0.13 13.5 8.52e-42

AIC: 396.7377
Warning messages:
1: In truncateToBinary(y) :
Some observations were > 1. These were truncated to 1.
2: In sqrt(diag(vcov(obj))) : NaNs produced

Thoughts? Thanks!

Whitney

> On Sun, Sep 25, 2011 at 7:10 PM, Whitney Wiest <whitney.wi...@gmail.com>

[Richard Chandler]

Hi Whitney,

Thanks for the link. I didn't know about user-defined contrasts.

The issue that you face now is that your estimate of psi is near a parameter boundary, ie psi is very close to 1. This makes it hard to find the MLEs of the and obtain SEs from the hessian.

You could try a few different things.

1. Use starting values and play around with some of optim()'s control arguments. See ?optim.

m <- occu(~1 ~1, sesp.occu, starts=c(5, 2), control=list(rel.tol=1e-30))
backTransform(m, type="state")

2. Use the (non-)parametric bootstrap to compute standard errors and confidence intervals.

3. It seems like you have count data, so you could avoid this issue by modeling abundance

Best of luck,

Richard


_____________________________________
Richard Chandler, post-doc
USGS Patuxent Wildlife Research Center
301-497-5696

From:	Whitney Wiest <whitne...@gmail.com>
To:	unmarked <unma...@googlegroups.com>
Date:	10/19/2011 01:19 PM
Subject:	[unmarked] Re: Pcount function - Missing estimates in global model
Sent by:	unma...@googlegroups.com

[Andy Royle]

hi all, sorry but I didn't follow the whole thread here but I think Richard hit on the main issue which is that psi_hat = 1.0.

Is it the case that all of the sites have detections? In that case, there is not much modeling to be done...

regards

andy

 

J. Andy Royle
Research Statistician

USGS Patuxent Wildlife Research Center

12100 Beech Forest Rd.
Laurel, MD 20708
http://profile.usgs.gov/professional/mypage.php?name=aroyle
andy_...@usgs.gov
phone: 301-497-5846
fax: 301-497-5545

Book: "Hierarchical Modeling and Inference in Ecology: The Analysis of Data from Populations, Metapopulations and Communities" by J. A. Royle and R.M. Dorazio.  

unmarked: A very useful R package for fitting certain hierarchical models using likelihood methods. Available from: http://cran.case.edu/web/packages/unmarked/index.html

-----unma...@googlegroups.com wrote: -----

To: unma...@googlegroups.com
From: Richard Chandler <rcha...@usgs.gov>
Sent by: unma...@googlegroups.com
Date: 10/20/2011 09:08AM
Subject: Re: [unmarked] Re: Pcount function - Missing estimates in global model

[Whitney Wiest]

Thanks to you both. Richard, your first suggestion worked great.

Yes, almost all of the sites have detections. We are comparing the
changes in occupancy vs. abundance estimates depending on number of
surveys using our count data. I think I'm all set now.

Thanks for the help.

Whitney

On Oct 20, 9:27 am, Andy Royle <aro...@usgs.gov> wrote:
> hi all, sorry but I didn't follow the whole thread here but I think Richard hit on the main issue which is that psi_hat = 1.0.
>
> Is it the case that all of the sites have detections? In that case, there is not much modeling to be done...
>
> regards
>
> andy
>
>
>  
>
> J. Andy Royle
> Research Statistician
> USGS Patuxent Wildlife Research Center
> 12100 Beech Forest Rd.

> Laurel, MD 20708http://profile.usgs.gov/professional/mypage.php?name=aroylean...@usgs.gov

> phone: 301-497-5846
> fax: 301-497-5545
> Book: "Hierarchical Modeling and Inference in Ecology: The Analysis of Data from Populations, Metapopulations and Communities" by J. A. Royle and R.M. Dorazio.  
> unmarked: A very useful R package for fitting certain hierarchical models using likelihood methods. Available from:http://cran.case.edu/web/packages/unmarked/index.html
>

> -----unma...@googlegroups.com wrote: -----To: unma...@googlegroups.com

> From: Richard Chandler <rcha...@usgs.gov>
> Sent by: unma...@googlegroups.com
> Date: 10/20/2011 09:08AM

> Subject: Re: [unmarked] Re: Pcount function - Missing estimates in global modelHi Whitney,Thanks for the link. I didn't know about user-defined contrasts.The issue that you face now is that your estimate of psi is near a parameter boundary, ie psi is very close to 1. This makes it hard to find the MLEs of the and obtain SEs from the hessian.You could try a few different things.1. Use starting values and play around with some of optim()'s control arguments. See ?optim.m <- occu(~1 ~1, sesp.occu, starts=c(5, 2), control=list(rel.tol=1e-30))backTransform(m, type="state")2. Use the (non-)parametric bootstrap to compute standard errors and confidence intervals.3. It seems like you have count data, so you could avoid this issue by modeling abundanceBest of luck,Richard_____________________________________

> Richard Chandler, post-doc
> USGS Patuxent Wildlife Research Center

> 301-497-5696From:Whitney Wiest <whitne...@gmail.com>To:unmarked <unma...@googlegroups.com>Date:10/19/2011 01:19 PMSubject:[unmarked] Re: Pcount function - Missing estimates in global modelSent by:unma...@googlegroups.comHi Richard and Matt,

> > This is all assuming your data is...
>
> read more »

[Matt Giovanni]

Related to this earlier thread, can anyone tell me how to use the relevel function with UMFs?  For example:

> umf=csvToUMF("BOBO.csv",long=FALSE,type="unmarkedFramePCount")
> umf$MT=relevel(umf$MT,ref="NH")
Error in umf$MT : $ operator not defined for this S4 class

Thanks everyone.

Matt

[Ian Fiske]

Matt,

You need to extract the proper data.frame. If MT is an observation
level covariate, do

obsCovs(umf)$MT <- relevel(obsCovs(umf)$MT, ref = "NH")

If MT is a site-level covariate, do

siteCovs(umf)$MT <- relevel(siteCovs(umf)$MT, ref = "NH")

Hope that helps,
Ian