1
$\begingroup$

The goal of this regression is to determine whether the amount of leaf disk that an insect consumed varied by what tree the leaf material came from. I'll acknowledge upfront that my coding is rarely pretty/efficient, but hopefully it works (usually).

  • Variables:
    • Response: pctrans; the percent of a 7 mm diameter leaf disk that was consumed. Values have been transformed to fit (0,1).
    • Explanatory: tree; a categorical (factor) variable of six tree types.

When I use betareg(), which as I understand it, is best suited to data of this sort, I get no significance:

model.beta <- betareg(pctrans ~ tree, data=BT.data, link="logit")
modelnull.beta <- betareg(pctrans ~ tree, data=BT.data, link="logit")
lrtest(model1.beta, modelnull.beta)

Results:

Call:
betareg(formula = pctrans ~ tree, data = BT.data, link = "logit")

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
-2.7716 -0.5800  0.0472  0.5351  3.5109 

Coefficients (mean model with logit link):
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.111504   0.069191 -16.064  < 2e-16 ***
treeBC3F3   -0.050940   0.095889  -0.531  0.59525    
treeD54     -0.279927   0.096470  -2.902  0.00371 ** 
treeD58     -0.034000   0.095716  -0.355  0.72242    
treeEllis1  -0.006764   0.095175  -0.071  0.94334    
treeQing     0.785992   0.094003   8.361  < 2e-16 ***

Phi coefficients (precision model with identity link):
      Estimate Std. Error z value Pr(>|z|)    
(phi)   3.5549     0.1352   26.29   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Type of estimator: ML (maximum likelihood)
Log-likelihood: 529.8 on 7 Df
Pseudo R-squared: 0.1105
Number of iterations: 20 (BFGS) + 2 (Fisher scoring) 

Likelihood ratio test

Model 1: pctrans ~ tree
Model 2: pctrans ~ 1
  #Df LogLik Df  Chisq Pr(>Chisq)    
1   7 529.82                         
2   2 460.70 -5 138.25  < 2.2e-16 ***

As I've been told, since the model is significantly worse than the null, no comparisons can be made between treatment means.

HOWEVER... If I run the same model using glm the model is significantly better than the null.

beta.glm <- glm(pctrans ~ tree, data=BT.data, family=quasibinomial)

Results:

Call:
glm(formula = pctrans ~ tree, family = quasibinomial, data = BT.data)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.94474  -0.38492  -0.08785   0.22725   1.80291  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.22601    0.07643 -16.042  < 2e-16 ***
treeBC3F3    0.06826    0.10660   0.640  0.52205    
treeD54     -0.33864    0.11312  -2.994  0.00281 ** 
treeD58     -0.19878    0.11062  -1.797  0.07260 .  
treeEllis1  -0.07763    0.10808  -0.718  0.47276    
treeQing     0.88596    0.09978   8.879  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.2069603)

    Null deviance: 307.54  on 1240  degrees of freedom
Residual deviance: 267.59  on 1235  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

Analysis of Deviance Table
Model: quasibinomial, link: logit
Response: pctrans
Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                  1240     307.54              
tree  5   39.951      1235     267.59 < 2.2e-16 ***

Where do I go from here?

$\endgroup$
3
  • $\begingroup$ This is not a programming, but a stats question. Accordingly, I have voted for migration to the appropriate stackexchange site. $\endgroup$
    – Roland
    Jun 27 '16 at 14:35
  • $\begingroup$ My apologies! I'll be happy to make the move myself! $\endgroup$ Jun 27 '16 at 14:43
  • $\begingroup$ Just wait a few more minutes. The move got 4 out of 5 votes already. $\endgroup$
    – Roland
    Jun 27 '16 at 14:45
3
$\begingroup$

The quasi-binomial and the beta regression model lead to rather similar results: The tree type has a highly significant effect on the response. The two tree types that clearly differ from the others are D54 (with a somewhat lower response) and Qing (with a much higher response). Even the coefficients on the logit scale are roughly comparable: -0.3 for D54 and +0.8 for Qing. Hence, I'm not sure why you think that there is a problem or inconsistency.

$\endgroup$
5
  • $\begingroup$ I was told that if a model has no significance (compared to the null model), then no conclusions could be made regarding the model. In the case of the beta regression, the model is significantly worse than the null (using a likelihood ratio test), so I assume that the significance of the D54 and Qing are meaningless in the context of that model. Edit: I'd be extremely happy to find out I was misinformed! $\endgroup$ Jun 27 '16 at 15:24
  • 1
    $\begingroup$ @RegalPlatypus Well, you would be well advised to listen to Achim rather than whoever your earlier advisor might have been. I only saw this answer after finishing mine and am relieved that we agree. I don't usually rely on "argument by authority" but if AZ is the authority, I'm almost certain. $\endgroup$
    – DWin
    Jun 27 '16 at 15:30
  • $\begingroup$ Thank you, @Achim Zeileis, I was misreading a crucial part of the results. $\endgroup$ Jun 27 '16 at 15:40
  • 2
    $\begingroup$ @RegalPlatypus: The likelihood ratio test always compares the null hypothesis that a simpler/restricted model is sufficient against a more complex model under the alternative. In your case the null model with just an intercept is clearly rejected. [And for DWin: So we can add "argument by majority" to "argument by authority" now. But I hope that the case is clear enough that RegalPlatypus will come to the right insights himself/herself as well.] $\endgroup$ Jun 27 '16 at 15:53
  • $\begingroup$ I had been used to using drop1() or anova() for liklihood ratio tests which gives residual deviance in the output. I'd been reading the Log Likelihood column in the lrtest() as residual deviance. Like a said, a silly mistake on my part. $\endgroup$ Jun 27 '16 at 16:03
2
$\begingroup$

Looking at these two sections of the output, I am scratching my head about what you "have been told":

#-------
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.22601    0.07643 -16.042  < 2e-16 ***
treeBC3F3    0.06826    0.10660   0.640  0.52205    
treeD54     -0.33864    0.11312  -2.994  0.00281 ** 
treeD58     -0.19878    0.11062  -1.797  0.07260 .  
treeEllis1  -0.07763    0.10808  -0.718  0.47276    
treeQing     0.88596    0.09978   8.879  < 2e-16 ***

#---------------
Model 1: pctrans ~ tree
Model 2: pctrans ~ 1
  #Df LogLik Df  Chisq Pr(>Chisq)    
1   7 529.82                         
2   2 460.70 -5 138.25  < 2.2e-16 ***

The global test for significance is "highly significant", since it is a "good thing" to have a p-value near zero (or equivalently to have the difference in loglikelihood be large). There are 2 particular trees whose beta coefficients are sufficiently different from the (Intercept) to also be "highly significant".

As far as I can see the two models are telling you the same story as far as whether the full and null models are different and also identifying the same two factors (trees) as the ones to focus on with further discussion.

$\endgroup$
1
  • 1
    $\begingroup$ Thanks, gentlemen. I withdraw my question in a cloud of shame. $\endgroup$ Jun 27 '16 at 15:38

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .