Appropriate GLM for proportion response variable

Question

I have a dataset that includes the number of apples harvested from 31 trees and the number of apples from each tree affected by different forms of pest damage. I also have insect abundance data (averaged from multiple samples throughout the season) associated with each of those trees, and I'd like to see if I can correlate abundance of certain insect taxa with pest damage types. Here's a simplified example of that data:

apple.example <- data.frame(treeID = paste("A", c(1:5)), 
     applecount = c(1,1,13,14,5), chew = c(0,0,2,2,1), 
     codling = c(0,0,8,6,4), spiders = c(1.875,2,1.286,2,1.67))

Each apple can have more than one kind of damage, but I don't have per-apple level data. Lots of trees had no apples, and of those that did there are plenty of 0% and 100% frequency of damage types.

apple.example[3:4] <- apple.example[3:4]/apple.example$applecount

Since I don't care about the counts of pest damage (confounded by yield) but rather their relative frequency, my initial thought was to convert these values into percentages and analyze them using a binomial GLM. When I do that, R complains that the response values are not all 0 or 1. A few posts I've seen suggest that isn't necessarily an issue, but what gives me pause is that the lsmeans predicted values are not bounded between 0 and 1 as a binomial distribution ought to be. Also, AICc model selection results are not finding relationships among the data even where I might expect one, so I'm not 100% confident this is a good way to model this data.

m1 <- glm(chew~spiders, data=apple.example, family="binomial")

Discussion here suggests that a beta regression is better for this kind of data, so I switched to that, incorporating the transformation here to circumvent the 0/1 error there.

y.transf.betareg <- function(y){
  n.obs <- sum(!is.na(y))
  (y * (n.obs - 1) + 0.5) / n.obs
}

m2 <- betareg(y.transf.betareg(chew)~spiders, data=apple.example)

Both of these feel like bad workarounds for something that ought to be straightforward, though? The betareg method provides more resolution between models in AICc output, but neither are particularly good. I guess I'm just looking for some reassurance that one of these is a valid approach with this data, or, ideally, a more appropriate method.

Answer 1

First, you should incorporate the number of trials in the model. I illustrate below using what I feel is a more elegant way to obtain the fractions, and creating new variables pchew and pcodling:

> apple.example <- data.frame(treeID = paste("A", c(1:5)), 
+                             applecount = c(1,1,13,14,5), 
+                             chew = c(0,0,2,2,1), 
+                             codling = c(0,0,8,6,4), 
+                             spiders = c(1.875,2,1.286,2,1.67))
> apple.example <- transform(apple.example, 
+                            pchew = chew/applecount, 
+                            pcodling = codling/applecount)
> m2 <- glm(pchew ~ spiders, weights = applecount, 
+           data = apple.example, family = "binomial")

> library(lsmeans)
> lsmeans(m2, "spiders", at = list(spiders = 1:3), type = "response")
 spiders      prob        SE df   asymp.LCL asymp.UCL
       1 0.1730220 0.1549970 NA 0.024423399 0.6361675
       2 0.1351918 0.0826546 NA 0.037635848 0.3845711
       3 0.1045869 0.1925788 NA 0.002070816 0.8679789

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale

Alternatively, you specify the same model via the success and failure counts (making the transform step unnecessary) as follows:

> m2 <- glm(cbind(chew, applecount - chew) ~ spiders, 
+           family = "binomial", data = apple.example)