# Weight function for binomial or negative binomial GLM

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.

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.

• 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). – Frans Rodenburg Oct 5 '17 at 12:22