What is the appropriate way to test for a shift in probability using multiple logistic regression?

Question

We have data on the day in which a butterfly pupates (forms a cocoon) in the summer/fall of 2 different pairs of years (2005-2006 vs 2009-2010). At the time that the pupa forms it can either be in diapause (suspended animation) or not.

Individual butterfly caterpillars were randomly selected from a wild population and observed in outdoor enclosures until they pupated. At this time the pupa was determined to be either in diapause or not.

We would like to know if there is a shift in the probability that a butterfly pupa will be in diapause or not on a given day between years. In other words, does a butterfly pupa have a significantly different probability of being in diapause later (or earlier) in 2005 vs 2006?

We have analyzed each pair of years using separate multiple logistic regressions with Day of pupation and Year as factors and are interpreting a significant Year effect as a significant shift in diapause probability on a given day between years.

Is this the appropriate way to approach this problem?

Thank you.

Data Format:

Individual   Diapause Pupation_Day Year
1            1          200        2005
2            0          201        2005
.            .           .          .
.            .           .          .
.            .           .          . 
n            1          300        2006

Answer 1

Your initial logistic regression approach seems reasonable, assuming good diagnostics. However, if I understand the scientific question right, you'll need an interaction between year and day to see differences between years in the tendency to be in diapause as day increases. The main effect of year only says that one year has a higher per day probability of being in diapause, which is not, I think, what you care about.

Alternatively

The other way to think about the problem is to have days as a dependent variable. This means you would reverse the conditional probability formulation that you implicitly started with by not asking 'are later pupae more likely to be in diapause?' i.e. P(diapause | day, year), but rather 'do pupae in dipause tend to occur later in the year' i.e. P(day | diapause, year).

In such a model the observations are Pupation_Day which is classified by Diapause status which is nested within Year. The effect of interest is a difference of differences, specifically: the difference between the difference in mean Pupation_Day for Diapause==1 and Diapause==0 in Year==2005 versus Year=2006. (I'm not still sure what role the 'pairs of years' play in your design though.)

The only potentially tricky bits I see are remembering that Diapause is a random effect, and making a sensible assumption about the conditional distribution of Pupation_Day. Frankly, I'd try Normal with shared variance first and then see if the model diagnostics object. If they do, I'd crack open JAGS/BUGS, and just write out the full model, which would make inference about difference in differences easier too.

This alternative approach may be too much machinery for the scientific question though. If we could tell what the ideal inferential endpoint was it would be easier to recommend an approach.