Factor significance and group level contrasts following analysis w/ pcount

[Steven Brady]

Hi All:

I have a data set comprised of 16 sites where I recorded point count observations five times during each of two consecutive years (160 observations total).  I have been successful in using pcount() from [unmarked] in R for  the initial analysis.  Specifically, I have run a model with two site level covs (1. % land development and 2. pond area) and one obs level cov (date) modeled in a way as to predict abundance at each site.  

 

I would like to run a secondary additional analysis in more of an anova style.  That is, rather than specifying the predictor as % land development (being a continuous variable), I would like to specify a factor composed of three groups (specifically, classes of land cover, e.g. FOREST; URBAN; AGRICULTURAL).  I have gotten the model to run this way, and to report the output in typical R fashion, e.g. 

 

fm2 <- pcount(~date ~ land_class + log10(area), data=duck)

A summary of the model produces the anova table (values omitted here for concision) for Abundance and Detection parameters:

Abundance:                                   

(Intercept)                   

classFOREST       

classURBAN        

log10(area)      

 

Detection: 

(Intercept) 

date           

 

As I understand it, the above output reflects the estimates, etc for FOREST and URBAN land covers, and their significance with regard to the overall intercept.  My two questions are these:

 

1. How would I evaluate the overall significance of the land_class factor?  I tried using anova(full model, reduced model) for a LL ratio test, but that does not seem to be supported.

2. How would I specify a set of contrasts to evaluate the significance of each of the groups within the land_class factor? e.g. how can I evaluate whether the mean abundance in FOREST differs from AGRICULTURAL from URBAN?  

 

Many thanks in advance for any help you can provide on this topic.

Best,

Steve

__________________________________________

Steven P. Brady, Ph.D. Candidate

School of Forestry & Environmental Studies

Yale University

370 Prospect Street

New Haven, CT 06511

Email: steven...@yale.edu

Phone: 203-432-5321 Fax: 203-432-3929 

Web: http://www.cbc.yale.edu/people/skelly/steveb.html 

[Andy Royle]

hi Steve,

 in both cases I think you'll have to use standard methods to obtain the result directly --  e.g., to test the significance of a treatment effect you can obtain the log-likelihood for the full and reduced models and the difference should have a chi-square on 2*(dffull-dfred) df. i.e., fit both models and extract the relevant stuff from the resulting object.

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: Steven Brady <steven...@yale.edu>
Sent by: unma...@googlegroups.com
Date: 11/30/2010 05:07PM
Subject: [unmarked] Factor significance and group level contrasts following analysis w/ pcount

[Matthew Giovanni]

Hi Steve-

I think your model with the land class categorical variable should allow you to compare abundance/density among category types by back-transforming your model coefficients and SE's to abundance predictions and CI's.  I don't think you can estimate a coefficient for land class per se, like you can for % land development, with N-mixture models because it's categorical, but I may be mistaken.

You might consider fitting models with land class and % development interacting (if you haven't already), so that you can assess effects of % development per category of land class.

So might fit some set of competing models such as:

null <- pcount(~date ~ 1 + log10(area), data=duck)
fm1 <- pcount(~date ~ land_class + log10(area), data=duck)

fm2 <- pcount(~date ~ devel + log10(area), data=duck)

fm3 <- pcount(~date ~ land_class*devel + log10(area), data=duck)

then select your most informative model with AIC:

models=fitList(null=null, fm1=fm1, fm2=fm2, fm3=fm3)
modSel(models)

and then make predictions from, for example, your interaction model:

DF=data.frame(land_class=rep(c("FOREST", "URBAN", "AGRICULTURAL"), each=25), devel=seq(0,100,4), area=mean size of your sampling plots)
predict(fm3, newdata=DF, type="state", se=TRUE, backTran=TRUE, appendData=TRUE)

and maybe predictions from your interaction model of abundance averaged across % development per land class type:

DF=data.frame(land_class=c("FOREST", "URBAN", "AGRICULTURAL"), devel=mean(data$devel), area=mean size of your sampling plots)
predict(fm3, newdata=DF, type="state", se=TRUE, backTran=TRUE, appendData=TRUE)

HTH.

Matt

____________________________________
Matt Giovanni
Postdoctoral Visiting Research Fellow
Canadian Wildlife Service
Prairie and Northern Region
Environment Canada
2365 Albert St., Room 300
Regina, SK S4P 4K1
306-780-5308 work
402-617-3764 mobile
Research website:
http://sites.google.com/site/matthewgiovanni/

[Richard Chandler]

Hi Steve,

A likelihood ratio test can be done in unmarked: 

> LRT(full.model, reduced.model)

To test for differences among all levels of a factor, you can change the reference level of your factor using relevel() and run the model again. Currently, it appears as though your reference level is AG and so you have already compared AG vs FOREST and AG vs. URBAN. To test for differences between FOREST and URBAN, try something like this:

> sc <- siteCovs(duck)
> sc$land_class <- relevel(sc$land_class, ref="FOREST")
> siteCovs(duck) <- sc

> pcount(~date ~ land_class + log10(area), data=duck)

You might also want to look at

> help(contrasts)


Richard

[brady]

All:

Thank you for the many useful suggestions. I tried to implement your
advice regarding log-likelihood tests, but encountered some errors.

I tried to extract the log-likelihood for each reduced model, but
received the following error

> logLik(fm2)
Error in UseMethod("logLik") :
no applicable method for 'logLik' applied to an object of class
"unmarkedFitPCount"

I also tried to compare via log-likelihood ratio test the full and
reduced models, and received an error about NaNs produced. (Yet, both
models were fit successfully without error.)

> LRT(fm2, fm2.1)
Chisq DF Pr(>Chisq)
1 -122.7824 -2 NaN
Warning message:
In pchisq(q, df, lower.tail, log.p) : NaNs produced

I would love to hear any additional feedback on this. Many thanks
again!

- Steve

On Nov 30, 5:56 pm, Andy Royle <aro...@usgs.gov> wrote:
> hi Steve,
>
>  in both cases I think you'll have to use standard methods to obtain the result directly --  e.g., to test the significance of a treatment effect you can obtain the log-likelihood for the full and reduced models and the difference should have a chi-square on 2*(dffull-dfred) df. i.e., fit both models and extract the relevant stuff from the resulting object.
>
> 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: Steven Brady <steven...@yale.edu>
> Sent by: unma...@googlegroups.com
> Date: 11/30/2010 05:07PM

> Subject: [unmarked] Factor significance and group level contrasts following analysis w/ pcountHi All:

> Web: <A...
>
> read more »

[Richard Chandler]

Hi Steve,

To extract the *negative* log-likelihood, you could use

> fm@negLogLike

I will add a logLik method in the next release. As for LRT, there seems to be a dumb mistake in my code, but I bet it will work if you change the order of your models:

> LRT(fm2.1, fm2)


Richard

[Steven Brady]

That did the trick!

Many thanks,

Steve