[R-sig-ME] How to vary residual variance with covariate and per factor in mcmcGLMM

[LOMMEN Suzanne]

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]]