I would recommend phrasing the question in terms of the analytic model and then finding the right software package--R or otherwise. The individuals are nested in rounds which are nested in sessions nested in treatments nested in fisheries. The individuals have attrition. I am not sure how you handle the implicit missing data created by this "bidding out" process. This could be a big problem for any model. This creates an intent to treat (the initial universe) vs. intent to treat for compilers problem (the ones who survive bidding out). Know which estimand you are interested in. This is a research and not an analytic question. Rosenbaum (2002), Observational Studies, has a good summary of this issue.

It sounds like the main dependency you wish to control for is the individual over time because the choice in time 2 is influenced by time 1 for that individual. So really, you are not measuring real time, but time to event or sequence. If not, as.factor(session) might solve your problem. However, why not just use a gee approach using an auto-regressive model for correlation and a fixed effect for the session? See Hubbard et al. (2010). Have you done a Hausman test to see if a random effects model explains the data better than a fixed effects model?

Using a factor would be an unconditional fixed effect approach. This would be fine if the sessions were unique. I am assuming that no one can switch sessions out of sequence. But if they cannot switch sessions, why can't you just recode the two variables? I am not sure from the example why you cannot create a session-round variable where each time period is expressed by a session round in sequence. So Session 1b Round 1 would Session.Round 1b.1 or 1, Session.Round 1b.2 would be 2 up to Session.Round 3.12 = 36. I am assuming that the sessions occur simultaneously or have mid-session crossovers so that would be not possible?

Thanks Richard, an update on this. We have tried to use a GEE approach as you advocate to deal with the autocorrelation. Our model is as follows

geeglm(Nash_decision~fishery+treatment+state_of_resource+period+Experience+Cumulative, family=binomial(link="logit"), data=data1, id=playerno., corstr = "ar1")

We recoded player and session (playerno.) to indicate a different player playing in each session (i.e. 1a - 8a for players 1-8 in session 1 and then 1b-8b for players 1-8 in session 2) and included period (i.e. round) in the model to cater for the time effect. We assume that playerno is the individual sample units in our dataset.

What we have found through testing however is that the QIC score is the same no matter what correlation matrix we specify (ar1 or exchangeable or independence).

geeglm(Nash_decision~fishery+treat+period+state_of_resource+Experience+Cumulative, family=binomial(link="logit"),

+ data=data1, id=playerno, corstr = "ar1")

> gee2 = update(gee1, corstr = "independence")

> gee3 = update(gee1, corstr = "exchangeable")

> sapply(list(gee1, gee2, gee3), QIC)

[,1] [,2] [,3]

QIC 871.02 871.02 871.02

Quasi Lik -434.02 -434.02 -434.02

Trace 1.49 1.49 1.49

px 1296.00 1296.00 1296.00

Also there is no Auto-correlation (checked with ACF) in each case, but if we remove period auto-correlation arises in some lags. Is this suggesting that the correlation between rounds is not so strong, so including period as factor is enough to deal with that? Is this an issue that the QIC is the same across the same models with different correlation structures?

Also do we need to consider the nesting factors (i.e. rounds nested within session)? We have read that you can ignore nested design through specifying the correlation structure correctly in GEE? Is this correct?

Thanks!