[R-sig-ME] fixed effects clarification/confirmation (ZIP model)
Ned Dochtermann
ned.dochtermann at gmail.com
Mon Jul 11 21:50:06 CEST 2011
Hi Jarrod,
Thanks for your reply, your description of the fixed effects is essentially how I was treating them, which is reassuring. Regarding transforming the coefficients, I had wanted to plot predicted values along with observed, for which I originally just used exp() but then realized that this was not necessarily appropriate given the random effects. However, since many of models don't have the slope (intercept models only), I hadn't properly considered that complication. For rather my benign plotting purposes I might therefore just stick with using exp() so as to also not have to deal with the additive overdispersion.
Your comment regarding the random term highlights a bigger issue though: my forgetting to specify the level for the random effects is a definite screw up. With an intercept only model the random term should be (I think):
random=~us(at.level(trait,1)):subject
with a slope-intercept model would it then be:
~us(at.level(trait,1) + at.level(trait,1):tl.diff.f):subject ?
Looking at the terms this doesn't seem like it should change the structure of G.
My confusion with the fixed effects design matrix stemmed from looking at the matrix for both the poisson and inflation portion, taking just the first half of the matrix like you describe gave me what I was expecting. Thanks for that, it was a bit disconcerting not seeing what I expected to see; though I should have remembered the inflation portion particularly since in retrospect that portion of the design matrix should have clued me in.
Some days I really wish I could run a simple t-test to address the questions I'm interested in. Even just normally distributed data would be nice!
Thank you very much for your help, it is greatly appreciated.
Ned
-----Original Message-----
From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk]
Sent: Monday, July 11, 2011 8:54 AM
To: Ned Dochtermann
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] fixed effects clarification/confirmation (ZIP model)
Hi Ned,
The estimated Poisson count for each of the four trt/phase combinations
can be obtained as:
trtA and fphaseA = traity + at.level(trait, 1):trtA
trtA and fphaseB = traity + at.level(trait, 1):trtA +
at.level(trait,1):fphaseB
trtB and fphaseA = traity
trtB and fphaseB = traity + at.level(trait, 1):trtB:fphaseB
at.level(trait,1):trtB:fphaseB is therefore the difference between trtB
in fphaseA and trtB in fphaseB.
MCMCglmm uses the standard call to model.matrix so standard R text-books
should give you an understanding of how the contrasts are formed.
The ZIP model is essentially multi-trait with the first trait being the
Poisson part. The fixed effect design matrix can be obtained as:
X<-maleF.zip at X[1:(dim(maleF.zip at X)[1]/2),]
which is often helpful.
The "standard" transformation of the fixed effects is complicated in
this analysis because you have a random intercept-slope model and so the
variance in the random effects is a function of the covariate. Also, is
it your intention to hold these effects constant over the Poisson part
and the zero-inflation part? Perhaps, using at.level(trait,1) within the
us() structure might be more appropriate?
Jarrod
On Thu, 2011-07-07 at 17:58 -0700, Ned Dochtermann wrote:
> Hi all,
>
> I constructed a mixed-effect model about which I'm primarily interested in
> the fixed effects and I wanted to make sure that I was interpreting the
> output appropriately. Based on the MCMCglmm course notes I specified the
> following model:
>
> maleF.zip<-MCMCglmm(attack.male~trait-1+at.level(trait,1):trt*at.level(trait
> ,1):fphase*at.level(trait,1):day.1+
> at.level(trait,1):temp.bar+at.level(trait,1):temp.dev,
> rcov = ~idh(trait):units,
> random=~us(1+tl.diff.f):subject,family="zipoisson",
> data=Compiled,prior=mf.prior,
> nitt=1080000,thin=480,burnin=120000,verbose=FALSE)
>
> trt and fphase each have two levels (A&B for both)
>
> The posterior mode estimates for the fixed effects from a shorter run are:
> posterior.mode(maleF.zip$Sol)
> traitattack.male traitzi_attack.male at.level(trait, 1):trtA
> 11.74390983 -2.52399057 -0.23913921
> at.level(trait, 1):fphaseB at.level(trait, 1):day.1
> at.level(trait, 1):temp.bar
> -0.14206666 0.02169066 -0.43889554
> at.level(trait, 1):temp.dev at.level(trait, 1):trtB:fphaseB
> at.level(trait, 1):trtB:day.1
> -0.05774492 -0.24961100 -0.03413810
> at.level(trait, 1):fphaseB:day.1 at.level(trait,
> 1):trtB:fphaseB:day.1
> -0.02103302 0.11741711
>
> My interpretation of the labeling of these effects is that
> traitzi_attack.male is the inflation latent factor. Of more interest to the
> questions I'm attempting to address, traitattack.male is the group mean for
> trtB with at.level(trait,1):trtA ( henceforth a.l():trtA ) being the
> difference between the group mean for trtA versus trtB. Picking an effect
> for example/clarification, I'm interpreting a.l():trtB:fphaseB as the
> difference from the global mean (traitattack.male) for trtB at fphaseB. Is
> this the correct interpretation? Also, since besides the inflation factor
> the effects are modeled as Poisson processes I've been assuming that the
> standard transformation of the fixed effects (averaged over random effects).
>
> The layout is a bit different from other mixed model outputs (i.e. trtA for
> the first contrast rather than trtB being carried throughout) so I wanted to
> double check. I tried extracting the sparse fixed effects design matrix for
> clarification but it wasn't what I was expecting and didn't end up helping
> me.
>
> Thanks for any help,
> Ned
>
> --
> Ned Dochtermann
> Department of Biology
> University of Nevada, Reno
>
> ned.dochtermann at gmail.com
> http://wolfweb.unr.edu/homepage/mpeacock/Ned.Dochtermann/
> http://www.researcherid.com/rid/A-7146-2010
> --
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list