Hi Ben,

Thank you very much for your reply. I  apologize for the cross-post.  I
posted to the regular list serve first and received a suggestion from
someone else that I post here instead. I'm not sure if it's possible to
remove my original post but I would be happy to do so if that is possible.
I should have been clear originally that I had posted on the regular
listserve first.

Your response was very useful in helping me move along in this analysis.  I
have addressed several of your comments in more detail below but to
summarize my email, these are my current questions:

1.  Is it possible for me to still use the nlme package?  If so, do I use
the pdBlocked structure to incorporate my random effects? And if I do use
this structure, could someone point me towards a resource to figure out how
to create my groupedData?  I'm sure my lack of understanding here is a
statisitical, rather than R code issue.
2.  It is unclear to me how to deal with the problem of having my
categorical variables occur in the fixed and random effects.  Any insight
on this would be helpful.
3.  Based on the paper you sent and a look at my data, it does seem like I
need to incorporate the trt interactions but this exacerbates the problem
of an overfitted model.
4.  I tried the model that you kindly provided and received an error.  See
below for details.

I was switching to lme4 because I didn't think that I could use the correct
random terms in nlme.  I would actually prefer to stick with nlme if it is
possible to do so. I have a few simpler models that I have already run in
nlme so it would be nice to not have to re-run all of them.  I understand
the p-values provided are a bit controversial but journals like them anyway
so they are nice to have.

My best guess for how to do this based on what I could figure out and your
suggested lmer model is the following

m1<-lme(y=trt*state*season*site, data=group,
random=pdBlocked(list(pdIdent(~state|ID), pdIdent(~season|ID-1))),
method=("REML"))

However, I have not been able to find any documentation on how to correctly
create the groupedData object.  To make my model run, I used

group=groupedData(y~1|state, data=d, order.groups=FALSE)

If there is a way to run this model correctly in nlme, I would greatly
appreciate your insight on this.

You also raised a couple of other issues in your email that I am not sure
how to address. You stated

>   In general, categorical predictors should not appear in both
> the fixed and random parts of the model -- that constitutes overfitting.
>
> My effect of interest in this model is actually the interaction between
treatment and state because I expect that the treated individuals but not
the control will decline in (scrotal) size after (contraceptive) treatment.
 Thus, assuming I am not missing something, state must be in my fixed
effects.  I expect that there will be seasonal differences in size due to
breeding season.  This seasonal difference is interesting but ancillary to
the main question so perhaps I could take season out of the fixed effects.

The other issue you raised was this one:


> Depending on the structure of your data you should also strongly
> consider including interactions between treatment and your random
> factors (see Schielzeth, Holger, and Wolfgang
> Forstmeier. 2009. “Conclusions Beyond Support: Overconfident Estimates
> in Mixed Models.” Behavioral Ecology 20 (2):416–420.
> doi:10.1093/beheco/arn145.
> http://www.ncbi.nlm.nih.gov/pubmed/19461866.).  However, you can only
> include treatment interactions that make sense.


I asked about this and he said that the interactions with ID were most
important but this paper does make a good argument for it.  It seems like
including these these would result in a very complicated random component
in my model. It seems to me that I would need to include the trt
interaction with season, site, and state (but not ID as you correctly
assumed that each individual is only given 1 treatment).  I admit that I am
already feeling pretty overwhelmed with the mixed modelling approach but
you do have a good point.

If I were to include these, my guess is that my model would then look like

m1<-lmer(y~trt*state*site*season+(state|ID)+(0+season|ID)+(site|trt)+(state|trt)+(season|trt)
)

Based on your point above, it seems like this model would be way overfit.


It sounds like each
> individual is allocated to a single treatment, so (guessing that your
> statistician means that you want the *interactions* of ID with state
> and season) I think you should probably use

y~trt*season*state*site+(state|ID)+(0+season|ID)
>

Yes, the interactions of ID with state and season is exactly what I am
looking for.
Also, your suggestion of a log transformation did take care of the problem
of the unequal variances so thankfully, I should not have to deal with that
in a more complicated way.
Unfortunately, the model as above but on log transformed data resulted in
the following error
In mer_finalize(ans) : singular convergence (7)

It is unclear to me why this is occurring.

Thank you again for any help you can provide me.  I greatly appreciate it!

Sara



>
>
>
> > I am also trying to do this
> > with 2 factors instead of one.
> > The overall structure of my data is
> > Structure of my dataset (d)
> >   ..@ frame   :'data.frame': 173 obs. of  10 variables:
> >   .. ..$ y     : int [1:173] 209 382 448 353 224 112 198 273 495 622 ...
> >   .. ..$ trt   : Factor w/ 2 levels "G","S": 1 1 1 1 1 1 1 1 1 1 ...
> >   .. ..$ season: Factor w/ 2 levels "Breeding","NonBreeding": 2 2 2 2 2
> 2 2
> > 1 1 1 ...
> >   .. ..$ state : Factor w/ 2 levels "Post","Pre": 1 1 1 1 1 1 1 1 1 1 ...
> >   .. ..$ site  : Factor w/ 2 levels "Mrak","Orchard": 1 1 1 1 1 2 2 1 1 1
> > ...
> >   .. ..$ G     : num [1:173] 1 1 1 1 1 1 1 1 1 1 ...
> >   .. ..$ ID    : Factor w/ 59 levels "100","103","105",..: 11 19 22 58 6
> 36
> > 50 11 20 24 ...
> >   .. ..$ S     : num [1:173] 0 0 0 0 0 0 0 0 0 0 ...
> >   .. ..$ Pre   : num [1:173] 0 0 0 0 0 0 0 0 0 0 ...
> >   .. ..$ Post  : num [1:173] 1 1 1 1 1 1 1 1 1 1 ...
> >
>
> > m1<-lmer(y~trt*season*state*site+(0+G|ID)+(0+S|ID) +(0+Pre|ID) +
> > (0+Post|ID)+(1|season), data=d)
>
>   Your basic problem here is that it doesn't make sense to test the
> interaction
> of treatment with ID (as you have effectively done here), because each
> individual is only in one treatment category ...
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]

