1
$\begingroup$

I am trying to investigate the influence of Sex on the response variable Prop_new_trees (proportion data, No_new_used/Total_trees), with joey_number (a factor) as a covariate. Because I have proportion data I have been trying to fit a binomial GLM in R.

> head(exp)
Possum Sex Prop_new_trees joey_number No_prev_used No_new_used Total_trees
1  Cadence   F     0.08333333           1           11           1          12
2 Chechnya   F     0.41666667           1            7           5          12
3     Coco   F     1.00000000           1            0          13          13
4    Comet   F     0.50000000           1            8           8          16
5  Cupcake   F     0.60000000           1           12          18          30
6  Delilah   F     0.25000000           2            3           1           4

If I run a basic binomial GLM

bb1 <- glm(Prop_new_trees ~ Sex + joey_number, family = binomial, data = exp

I get an error message saying non-integer successes. I did a bit of reading and I think I need to add in the weights for each individual - ie. the number of observations per individual (Total_trees) which has been used to calculate the proportions. If I run this:

bb <- glm(No_new_used ~ Sex + joey_number, 
        weights = Total_trees, data = exp)

> summary(bb)

Call:
glm(formula = No_new_used ~ Sex + joey_number, data = exp, weights = Total_trees)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-31.295   -8.189    0.471    5.927   43.631  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   13.200      2.489   5.304 4.05e-05 ***
SexM          -3.767      2.229  -1.690   0.1073    
joey_number   -3.166      1.223  -2.587   0.0181 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 307.1152)

    Null deviance: 9051.9  on 21  degrees of freedom
Residual deviance: 5835.2  on 19  degrees of freedom
AIC: 143.31

Number of Fisher Scoring iterations: 2

> Anova(bb)
Analysis of Deviance Table (Type II tests)

Response: No_new_used
            LR Chisq Df Pr(>Chisq)   
Sex           2.8573  1   0.090960 . 
joey_number   6.6946  1   0.009671 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

My question is, have I interpreted the weights function correctly?

Secondly, my data is positively skewed (there are several zero values, which are all biologically meaningful). If I try using a negative binomial my results become more significant.

enter image description here

enter image description here

bb <- glm.nb(No_new_used ~ Sex + joey_number, 
        weights = Total_trees, data = exp)

> summary(bb)

Call:
glm.nb(formula = No_new_used ~ Sex + joey_number, data = exp, 
    weights = Total_trees, init.theta = 2.170196952, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-7.2052  -2.9716  -0.9308   1.5414   5.2810  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.22237    0.12152  26.517  < 2e-16 ***
SexM        -0.60528    0.10636  -5.691 1.26e-08 ***
joey_number -0.84772    0.07177 -11.811  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(2.1702) family taken to be 1)

    Null deviance: 454.88  on 21  degrees of freedom
Residual deviance: 256.19  on 19  degrees of freedom
AIC: 1295.9

Number of Fisher Scoring iterations: 1


              Theta:  2.170 
          Std. Err.:  0.289 

 2 x log-likelihood:  -1287.950 
> Anova(bb)
Analysis of Deviance Table (Type II tests)

Response: No_new_used
            LR Chisq Df Pr(>Chisq)    
Sex           32.176  1  1.408e-08 ***
joey_number  145.946  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Is it valid to use a negative binomial without specifying theta, and is negative binomial the better option for my skewed data? (I tried an arcsine transformation but that didn't really help). I didn't have any luck with a poisson distribution because I have different total number of observations for each individual (hence using a proportion), although I don't doubt that there is a way to factor this in somehow.

Thanks in advance.

$\endgroup$
  • 1
    $\begingroup$ Try glm(cbind(No_prev_used, No_new_used) ~ Sex + joey_number, family = binomial, data = exp) instead. You won't need weights if you model as cbind(success, failure). $\endgroup$ – Frans Rodenburg Oct 5 '17 at 12:22
  • 1
    $\begingroup$ If you want to use Poisson or negative binomial best to make the response the number of events and add log(total) as an offset. As you surmise for Poisson this is no advance on the binomial.. $\endgroup$ – mdewey Oct 5 '17 at 12:25
  • $\begingroup$ Thanks - the cbind has worked perfectly and spits out a significant p-value for both response variables. I ran a test for dispersion using the P__disp function in msme and it shows the model is overdispersed (dispersion value ~4). I have changed to quasibinomial and now the result is non-significant for Sex (binomial p=0.0019, quasibinomial p=0.12). I'm still a bit concerned that the data is not normally distributed - I'm not sure how strict this is for a binomial GLM. I would be keen to hear anyone's suggestions based on the histogram I've included. Thanks! $\endgroup$ – ecoH Oct 5 '17 at 12:43

Your Answer

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

Browse other questions tagged or ask your own question.