How to analyse count data that include 0 as a category?

Question

I'm testing the yield of corn subject to two treatments: Temperature (Cold and Warm) and Light Color (Red, Blue, Natural.) The number of plants per plot that produced no cobs, one cob, two cobs and three cobs were counted.

#Several replicates were made, but only one plot per combination of treatments 
#is shown in this example.
temp = c('cold','cold','cold','warm','warm','warm') 
light = c('red','blue','natural','red','blue','natural')
zer_cob = c(8,7,3,0,2,0)
one_cob = c(5,2,7,1,0,2)
two_cob = c(0,2,4,8,7,0)
thr_cob = c(1,0,0,9,8,6)
total_plants = zer_cob + one_cob + two_cob + thr_cob
data = data.frame(temp, light,zer_cob,one_cob,two_cob,thr_cob,total_plants) 

Twenty plants were planted in each plot, but only the number in total_plants survived at the end - so the number of plants per plot varies. I want to determine if the treatments produce different yield. My initial plan was to obtain a weighted average of produced cobs for each row and perform an ANOVA, but this will be incomplete as the number of plants per plot differ and I don't want to lose that information (and cold has a higher plant count in the zero-one bracket, while warm does in the two-three bracket). Someone has a recommendation of how to better analyse this dataset?

Answer 1

I am not sure whether this should be counted as a comment, as I still ask for detail, or as an answer, as I propose a method and code.

In a first step you will have to define what you consider to be yield. Counting cobs I assume you want to yield cobs. As you have counted cobs, not weighed them, the count of cobs is the natural $y$. If you also consider plants as yield you either need to find a good way to weigh cobs against plants or do two different calculations.

Cobs come in whole numbers (0, 1, 2, ...) and there is no theoretical upper bound so a Poisson process is a good first candidate. To do Poisson regression in R you'll first need to compute the number of cobs as in

cobs <- one_cob + 2*two_cob+ 3*thr_cob

From which follows the following code:

temp = c('cold','cold','cold','warm','warm','warm') 
light = c('red','blue','natural','red','blue','natural')
zer_cob = c(8,7,3,0,2,0)
one_cob = c(5,2,7,1,0,2)
two_cob = c(0,2,4,8,7,0)
thr_cob = c(1,0,0,9,8,6)
total_plants = zer_cob + one_cob + two_cob + thr_cob
cobs = one_cob + 2*two_cob+ 3*thr_cob
data = data.frame(temp, light,zer_cob,one_cob,two_cob,thr_cob,total_plants, cobs) 

mod <- glm(cobs ~ temp * light, data = data, family = poisson)
summary(mod)

plot(data$cobs, predict(mod, type = "response"), xlab = "observed cobs", ylab="predicted cobs")

with Residual deviance: 1.5987e-14 can we assume, these are not really observed data but made up for the purpose of asking the question?