How to vary residual variance with covariate and per factor in mcmcGLMM

Fri, Mar 31, 2017 2:15 AM

Many thanks!
Suzanne

From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk]
Sent: 30 March 2017 07:44
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

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><mailto:j.hadfield at ed.ac.uk>
Sent: 28 March 2017 17:59:43
To: LOMMEN Suzanne; r-sig-mixed-models at r-project.org<mailto: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><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/><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/><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<mailto:R-sig-mixed-models at r-project.org> mailing list

https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

How to vary residual variance with covariate and per factor in mcmcGLMM

Thread (2 messages)