2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ The LS means will be on the same scale as the model's predictions -- which will be the logit scale, not the percentage scale. (Look at the annotations below the output which say just that.) I you add type = "response" to the lsmeans call, they will be in range. $\endgroup$ – rvl Feb 28 '17 at 22:15
  • $\begingroup$ This works great, thanks for pointing out that simple omission. Is it safe to take this as a tacit endorsement of the binomial approach for this sort of data? $\endgroup$ – Megachile Mar 2 '17 at 18:26
  • $\begingroup$ I posted a corrected version of my second comment (since deleted) as an answer $\endgroup$ – rvl Mar 2 '17 at 19:11
3
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ This is great, thanks--exactly the simple solution I was looking for. One follow up question: is it valid to use applecount as a nuisance variable in a model like the last one you suggested? When the response variable is a percent it's distanced from the absolute count but it seems sketchy for count to appear in both the predictor and response variables. There are a couple of cases where including it as a predictor changes the AICc rankings and I want to get it right. $\endgroup$ – Megachile Mar 6 '17 at 18:27
  • $\begingroup$ You're just conditioning on applecount and in principle I don't see any issue with also using it as a predictor. I'm more concerned with the solution I gave from the standpoint of it assumes no variations other than those inherent in the binomial model -- and it seems reasonable to believe that trees vary somewhat from one another. I think people fit something called an overdispersion model in such cases but I'm no expert. I hope somebody else will comment. $\endgroup$ – rvl Mar 6 '17 at 22:48

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.