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.
type = "response"
to thelsmeans
call, they will be in range. $\endgroup$