Logistic regression on binary response or ANOVA on proportions/percentages

Question

I am analyzing survival of seedlings in Styrofoam blocks (known as styroblocks). These blocks have a certain size and contain cavities in which the seedlings are planted. All cavities within a block have the same dimensions (diameter and depth). Different blocks come with different cavity dimensions. I am analyzing whether there is an influence of cavity size on survival. I am using 5 different styroblocks with 5 different cavity volumes and plant seedlings in all cavities. For the analysis, however, I only consider 15 seedlings from the center of each block. Furthermore, I plant 3 different plant varieties and repeat each variety X styroblock combination 3 times. Lastly, this experiment was replicated at 2 different nurseries to test for nursery effects.

Here’s an example dataset:

xx <- data.frame(Nursery = rep(c("Nursery A", "Nursery B"), each = 675),
Styroblock = rep(c("Block A", "Block B", "Block C", "Block D", "Block E"), 6, each = 45),
Variety = rep(c("Variety A", "Variety B", "Variety C"), 2, each  = 225),
Replicate = rep(c(1,2,3), 30, each = 15), Survival = sample(c(1,0),1350, replace=T, prob=c(.85,.15)))

Here’s my first option:

fit = glm(Survival ~ Nursery * Variety * Styroblock, data =xx, family=binomial(link="logit"))
summary(fit) ## at this point some model simplification should be performed
library(car)
Anova(fit, type="II", test="Wald")

HOWEVER, since cavities within styroblocks are spatially dependent, for my other measurements such as height (not shown here), I took an average across all 15 seedlings to get one measurement per styroblock (which avoids potential pseudo-replication). Given this, I summarized dead/alive (binary) by styroblock and clone, resulting in a proportion (or percentage) instead. This also reduces my sample size from 1350 to 90. See here:

require(plyr)
xx.sum<-ddply(xx, .(Nursery, Styroblock, Variety, Replicate), summarise, Survival = sum(Survival)/15)

Question 1: Is this the correct way to do?

If yes, I am not sure if I can follow up with ANOVA on the proportions/percentages due to unequal variances:

require(ggplot2)
ggplot(xx.sum, aes(x=Styroblock, y=Survival))+geom_boxplot()+facet_grid(Nursery~Variety)

Once option would be to transform, which I tried using arcsine and square root transformation but this does not help.

Question 2: How would you suggest proceeding in this case, i.e. which test is most appropriate to understand whether my predictors or predictor combinations influence survival?

Thanks!

Answer 1

Answering my own questions which I asked a couple of months ago:

No, use a GLMM with a random effect to account for the repeated measures on container, and a binomial distribution for your survival data.

You just need to add a container ID:

xx$ContainerID <- with(xx, paste(Variety, Styroblock, Replicate, sep="-"))

require(lme4)
fit <- glmer(Survival ~ Nursery * Variety * Styroblock + (1|Nursery:ContainerID),
             data = xx, family = binomial(link = "logit"))
summary(fit)

# TO GET P-VALUES FOR ALL FIXED EFFECTS (however this will take some time)
require(afex)
mixed(Survival ~ Nursery * Variety * Styroblock + (1|Nursery:ContainerID),
             data = xx, family = binomial(link = "logit"), method = "LRT")

Answer 2

my opinion is to not overthink it

There are 30 categories defined by 2(nurseries) * 5(blocks) * 3(variety). In the simulated example there is a lot of data. So you can get a good estimate for each category.

Next, what is your question? if you are trying to predict survival based on your 3 variables then just look at the category the plant is in and give this as you best estimate.

if you want to know if nursery does not influence the survival you could simply compare two models one with and one without the nursery variable,

i'd clarify the question to get something more precise