1
$\begingroup$

I have trouble with interpreting the output of a general linear model, which I have fitted to some data. Here are my data:

     Gall  Site healthy unhealthy
1  absent Site1    1750        35
2  absent Site2    1642       146
3 present Site1    1146        23
4 present Site2     333        30

I have made a barplot to visualize this better:

My barplot

I want to test whether healthy or unhealthy (saplings) was induced by the galls being present or absent. The saplings were randomly collected in each site. Therefore, I fit a generalized linear model to these data:

abc <- glm(cbind(healthy,unhealthy)~Site*Gall, family="binomial", data=f1)

I summarize my model

summary(abc)

Call:
glm(formula = cbind(healthy, unhealthy) ~ Site * Gall, family = "binomial", 
    data = f1)

Deviance Residuals: 
[1]  0  0  0  0

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)            3.912023   0.170713  22.916  < 2e-16 ***
SiteSite2             -1.491959   0.191314  -7.798 6.27e-15 ***
Gallpresent           -0.003484   0.271097  -0.013    0.990    
SiteSite2:Gallpresent -0.009634   0.342474  -0.028    0.978    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.1069e+02  on 3  degrees of freedom
Residual deviance: 4.5963e-14  on 0  degrees of freedom
AIC: 30.235

Number of Fisher Scoring iterations: 3

Now I try to simplify the model:

drop1(abc, test="Chi")

Single term deletions

Model:
cbind(healthy, unhealthy) ~ Site * Gall
          Df   Deviance    AIC        LRT Pr(>Chi)
<none>       0.00000000 30.235                    
Site:Gall  1 0.00079137 28.236 0.00079137   0.9776

I omit the two-way interaction between Site and Gall:

> f1.rid1 <- update(abc, .~. - Site:Gall)

I test each 1-way term, each at a time:

> drop1(f1.rid1, test="Chi")
Single term deletions

Model:
cbind(healthy, unhealthy) ~ Site + Gall
       Df Deviance     AIC     LRT Pr(>Chi)    
<none>       0.001  28.236                     
Site    1  104.027 130.262 104.027   <2e-16 ***
Gall    1    0.004  26.239   0.003   0.9542    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model tells me to omit Gall, so my final model is

> f1.rid2 <- update(f1.rid1, .~., - Gall)
> summary(f1.rid2)

Call:
glm(formula = cbind(healthy, unhealthy) ~ Site + Gall, family = "binomial", 
    data = f1)

Deviance Residuals: 
        1          2          3          4  
-0.014033   0.007093   0.017279  -0.015671  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.914420   0.148106  26.430   <2e-16 ***
SiteSite2   -1.494968   0.158755  -9.417   <2e-16 ***
Gallpresent -0.009518   0.165722  -0.057    0.954    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1.1069e+02  on 3  degrees of freedom
Residual deviance: 7.9137e-04  on 1  degrees of freedom
AIC: 28.236

Number of Fisher Scoring iterations: 3

My question is:

Why am I dropping Gall in my model? I know the output tells me to do so, but by looking at the initial barplot, I figured the presence or absence of Gall must have a significant effect on healthy or unhealthy.

$\endgroup$
  • $\begingroup$ why are you running generalised linear model for categorical predictors? $\endgroup$ – Adam Quek May 13 '17 at 7:29
  • $\begingroup$ Why not ? I am trying to figure out wether I can drop site or gall for my final model which would explain the counts $\endgroup$ – RLover May 14 '17 at 13:18
  • $\begingroup$ You don't need to use drop1(); note that the p-values are the same as in the summary() output. Also, there is an error in your update(f1.rid1, .~., - Gall), the 2nd comma shouldn't be there. That's why Gallpresent is still in the last model. However, you don't need to drop it just because it isn't significant. $\endgroup$ – gung - Reinstate Monica May 17 '17 at 15:24
3
$\begingroup$

Calculating the proportion unhealthy in each group and visualizing that might help you understand the model better. There is no effect of Gall and a relatively strong effect of Site with no interaction present:

data1 <- data.frame(matrix(c( 0,0,1750,35, 
                              0,1,1642, 146,
                              1,0,1146,23,
                              1,1,333,30  ), 4, 4,byrow=TRUE))

names(data1) <- c("Gall", "Site","healthy","unhealthy")

data1$prop <- data1 $unhealthy/ ( data1 $healthy + data1 $unhealthy )

data1

>      Gall Site healthy unhealthy       prop
>    1    0    0    1750        35 0.01960784
>    2    0    1    1642       146 0.08165548
>    3    1    0    1146        23 0.01967494
>    4    1    1     333        30 0.08264463

library(ggplot2)

ggplot (data1, aes(Gall, Site, y=prop, fill=Gall)) +
        geom_bar(stat="identity") + 
        facet_grid (~Site)

GGplot Image

Indeed, you can get these proportions directly from your fitted model as well.

glmfit1 <- glm( cbind(healthy,unhealthy) ~ Site * Gall, 
                      family = "binomial", data = data1  )

1- predict(glmfit1, type = "response")

>          1          2          3          4 
> 0.01960784 0.08165548 0.01967494 0.08264463

So, your variable selection procedure rightly drops Gall and the Site:Gall interaction term as they contribute nothing to the model. Though, I would question the need to do any variable selection at all here as the model interpretation was sufficiently clear in the initial (full) model.

$\endgroup$
2
$\begingroup$

@Brett has done a good job addressing the issue of why the model suggests Gall should be dropped. I gather your question is motivated by an apparent inconsistency between your original plot and the modeling results. Let me address that.

You state in the comments that, "I am trying to figure out wether I can drop site or gall for my final model which would explain the counts". Your data are counts in some sense, but not in the sense we normally mean by counts (e.g., Poisson data). Instead, you have counts out of a known total. That means your data are binomial. You do correctly specify the binomial distribution in the code for the model fits, but your barplot does not account for the known totals in displaying the counts. That is why it gives the impression that there is an interaction and an effect of Gall. @Brett's plot does display the counts as proportions of a known total, and so does not create that misleading perception.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.