[R-sig-ME] How to vary residual variance with covariate and per factor in mcmcGLMM
Jarrod Hadfield
j.hadfield at ed.ac.uk
Thu Mar 30 07:43:33 CEST 2017
cc-ed back to the list...
DIC isn't meaningful in this context - I would ignore it.
My suggestion assumes that the variance increases to the square of the
covariate. Sometimes people square-root the covariate (if it is
positive) if they want a positive linear relationship. A negative
relationship is more tricky. You could have 1/logsize1 but this is going
to behave badly at values close to zero, or perhaps 1/size1 or
1/sqrt(size1).
Cheers,
Jarrod
On 29/03/2017 21:13, LOMMEN Suzanne wrote:
>
>
> Dear Jarrod,
>
>
> Many thanks for the swift and helpful answer!
>
>
> Implementing your suggestions works but indeed gives very awkward
> values and even a negative DIC (which I never experienced before).
>
>
> Defining instead only the random factors (random = ~ site + subplot +
> idh(logsize1:treatment:year):units) and no rcov gives much better
> results, but in both cases the estimates of the variances increase
> with the covariate logsize1, while the data suggests it should
> decrease at high values. I assume this increase with higher covariate
> values is not a build-in assumption, but rather based on fitting my
> data? (in which case the fit is not good)?
>
>
> The DIC values of the more reasonable alternative described above is
> 2980, while the DIC of my model previously defined (variance invarient
> for the covariate, but differing in value for each treatment*year
> combination defined in rcov) gives a DIC of ca 4900. I know I cannot
> compare these at all since the R- and G-structure are different, but
> is it worrying to get values that differ so extremely from each other
> while fixed effects are the same?
>
>
> Cheers,
>
> Suzanne
>
>
>
>
> ------------------------------------------------------------------------
> *From:* Jarrod Hadfield <j.hadfield at ed.ac.uk>
> *Sent:* 28 March 2017 17:59:43
> *To:* LOMMEN Suzanne; r-sig-mixed-models at r-project.org
> *Subject:* Re: [R-sig-ME] How to vary residual variance with covariate
> and per factor in mcmcGLMM
>
> Hi,
>
> 1/ This is correct
>
> 2/ add idh(logsize1:treatment:year):units to the random effects. This
> gives 6 variances which we will call Vs1 ... Vs6 (s stands for slope)
> and these complement your 6 variances specified in rcov. We'll call
> these Vi1 ... Vi6 (s stands for intercept). The residual variance for
> the first treatment/year category (for example) at a specific value of
> logsize1 is then Vi1+Vs1*logsize1^2.
>
> 3/ You will find it hard to get precise estimates of these variance
> functions, and so if you can safely simplify it I would.
>
> Cheers,
>
> Jarrod
>
>
>
> On 28/03/2017 15:39, LOMMEN Suzanne via R-sig-mixed-models wrote:
>> Dear all,
>>
>>
>> I am grateful to the active community on this mailing list! Reading the posts already helped me solving many pRoblems.
>>
>>
>> I am quite new to mcmcGLMM and haven't managed to find out how to specify that the residual variance varies with explanatory factors AND a continuous covariate. I hope someone could advice me on the case below. My apologies for the long explanation.
>>
>>
>> I investigate how a herbivore treatment throughout the growing season ("treatment") affects individual plant size at the end of the season (log size2). The treatment was applied at 4 sites, and data were collected from ca 10 plants in each of ca 12 permanent subplots per site, in each of 3 years. The plant is an annual, so plants are only measured in 1 single year.
>>
>> Because not all plants survived , the total number of datapoints is 843.
>>
>>
>>
>> Other explanatory factors I included in my model are "year" (3 levels), the continuous covariate plant size at the start of the season ("log size1"), the continuous covariate "baresoil" (as the fraction of soil not covered by vegetation in the subplot), and their interactions.
>>
>> I also include the random factors "subplot" nested within "site" (although in mcmcGLMM it is not really nested, as I read in other posts).
>>
>> (Note: I included "year" as a fixed effect and not as a random effect because plant growth and the treatment effect varied dramatically from year to year, and this is what I want to capture in my analysis.)
>>
>>
>> I am interested in the effect on mean plant size AND in the variance, because I will use these two parameters to model a probability density function of logsize2 given logsize1 for each of all treatment*year combinations (2 treatments * 3 years = 6 combinations).
>>
>>
>> I noticed that variance varies per treatment*year combination, and in some of these combinations it also decreases with an increase in the covariate logsize1. Therefore, I would like to specify the residual variance per treatment*year combination as a function of logsize1.
>>
>>
>> So far, I have specified the variance as a constant for each of the treamnet*year combinations:
>>
>> rcov= ~idh(treatment*year):units,
>>
>> and adjusted the prior accordingly:
>>
>> R = list(V = diag(6), n = 0.002).
>>
>>
>> The complete model looks like this:
>> priorY1q1 = list(R = list(V = diag(6), n = 0.002),
>> G = list(G1 = list(V = 1, n = 0.002), G2 = list(V = 1, n = 0.002)))
>> modelY0bq1 <- MCMCglmm(logsize2 ~ logsize1 + treatment + baresoil + year +
>> treatment:logsize1 + baresoil:logsize1+
>> year:treatment + year:logsize1 + year:baresoil +
>> year:treatment:logsize1 + year:baresoil:logsize1 ,
>> random = ~ site + subplot, family = "gaussian",
>> rcov= ~idh(treatment*year):units,
>> data =pd_grow, prior = priorY1q1, burnin = 5000, nitt = 805000, thin=200)
>>
>>
>> I would be very glad to get comments on
>>
>> 1. if this is correctly specified so far (the results seem to make sense)
>>
>> 2. how to specify that the variance ALSO varies with logsize1, but possibly in a different way for each treatment*year combination?
>>
>> 3. If there is a risk to over-specify the variance this way?
>>
>>
>> Many thanks,
>>
>>
>> Kind regards,
>>
>> Suzanne Lommen
>>
>> PostDoc<http://www.unifr.ch/ecology/groupmueller/group/post-doctoral-associates/suzanne-lommen> in the Plant Population Biology group of Prof. Dr. Heinz M�ller-Sch�rer
>>
>> Co-management of the EU-COST Action 'SMARTER<http://www.ragweed.eu/>': Sustainable management of Ambrosia artemisiifolia in Europe
>> Coordination of the SMARTER Task Force Population Dynamics<http://ragweed.eu/task-force-population-dynamics/>
>>
>> Department of Biology/Ecology & Evolution
>> University of Fribourg/P�rolles
>> Chemin du Mus�e 10
>> CH-1700 Fribourg, SWITZERLAND
>> Tel: + (41) (0) 26-300 88 48 direct
>> Tel: + (41) (0) 26-300 88 50 secr.
>> Fax: + (41) (0) 26-300 9741
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20170330/a8ae77e3/attachment.pl>
More information about the R-sig-mixed-models
mailing list