[R-sig-eco] Inference, logistic regression

[Simon Blomberg]

On Mon, 2008-06-02 at 21:38 -0700, Andrew Rominger wrote:
> Dear list,
> 
> Please pardon this beginner's-level question, I feel it's not quite up  
> to the same caliber as recent discussions.
> 
> I'm working with a simple logistic regression model comparing the  
> presence/absence of an insect species against an index of plant  
> species turnover:
> 
> > foo<-glm(bout.psol$pres.de~bout.psol$index,family=binomial)
> 
> The term bout.psol$pres.de is binary 0,1; and bout.psol$index is continuous.
> 
> I'd like to use a likelihood ratio statistic to test the significance  
> of this regression, but I'm a little uncertain as how to proceed.   
> When I call summary(foo), I get...
> 
> Call:
> glm(formula = bout.psol$pres.de ~ bout.psol$index,
>      family = binomial)
> 
> Deviance Residuals:
>      Min       1Q   Median       3Q      Max
> -1.7180  -1.1289   0.6314   1.0323   1.7499
> 
> Coefficients:
>                       Estimate Std. Error z value Pr(>|z|)
> (Intercept)           0.30584    0.23095   1.324  0.18542
> bout.psol.edit$index  0.04552    0.01439   3.163  0.00156 **
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> 
> (Dispersion parameter for binomial family taken to be 1)
> 
>      Null deviance: 130.14  on 93  degrees of freedom
> Residual deviance: 118.17  on 92  degrees of freedom
> AIC: 122.17
> 
> Taking (Null dev) - (Redid dev), I get 11.97, which I assume to be  
> equal to -2log(L,full/L,reduced).  That's the desired test statistic,  
> so is it as simple as calling:
> 
> > pchisq(11.97,df=92)
> [1] 2.911346e-25          ?
> 

The degrees of freedom for the chisq test is 93 - 92 = 1.
 pchisq(11.97, 1, lower.tail=FALSE)
[1] 0.0005406394

Which is pretty close to the Wald test. Wald tests can sometimes be
misleading, since the estimate of the standard error in the denominator
can blow out and cause the test to be not significant even when there is
a big effect size, which is one reason to prefer the LR test. If you
want to get R to do the LR test, fit one model with and without the
covariate, and use anova.

fit.with <- glm(bout.psol$pres.de~bout.psol$index,family=binomial)

fit.without <- glm(bout.psol$pres.de~ 1 ,family=binomial)

anova(fit.without, fit.with)

Cheers,

Simon.

> That's an awfully small p-value, I think I'm interpreting something  
> wrong.  Any advice would be very welcomed.
> 
> Thanks very much in advance
> Andy Rominger
> 
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
-- 
Simon Blomberg, BSc (Hons), PhD, MAppStat. 
Lecturer and Consultant Statistician 
Faculty of Biological and Chemical Sciences 
The University of Queensland 
St. Lucia Queensland 4072 
Australia
Room 320 Goddard Building (8)
T: +61 7 3365 2506
http://www.uq.edu.au/~uqsblomb
email: S.Blomberg1_at_uq.edu.au

Policies:
1.  I will NOT analyse your data for you.
2.  Your deadline is your problem.

The combination of some data and an aching desire for 
an answer does not ensure that a reasonable answer can 
be extracted from a given body of data. - John Tukey.