[R-sig-ME] random intercept and random slope

[Douglas Bates]

2012/2/9 André Barbosa <andre.frainer at emg.umu.se>:
> Dear list members,
> I have read several threads on this list about the use of random variables and its interpretation. I seem to have learned a lot about model fitting, from plotting the raw data and checking the slopes and intercepts, to getting rid of the p-value mindset in which I had had my basic statistics courses at university.
> However, since my statistical courses never covered Mixed Effect Models, I am still unsure if I am doing the right thing or not. Having said that, would you please take a look at the following data and see if my rational is correct? The model is quite simple, I believe.
> I have 8 different plant species, which were mixed two-by-two in all possible combinations. So, each species has 7 pairs + it being alone (monoculture). My response variable is a ratio (observed /expected productivity values, where observed is the productivity of species “a” achieved when mixed with another species, and expected is its value when in monoculture).
> Each species had its own nutrient content analyzed. As I had three nutrient variables measured from each plant, I calculated indices of dissimilarity for each of those pairs.
> My starting model (without specifying random or fixed effects) would be:
> ratio ~ dissimilarity | species
> I expected, based on previous studies, that the relationship between ratio and dissimilarity would yield different slopes for each species, from negative to positive – expectation confirmed by potting a xyplot function of my data. Thus, species should be random. Looking at the xyplot of my data, I also see that the intercepts are somehow variable, ranging between 0.5 and 1.5 (response data points do not extend much further from this range, either). For this reason, I thought on including intercepts as random, as well, which leave me without fixed variables.
> So, I decided to test:
> lmer(ratio ~ 1 + (dissimilarity|species))
> Here follows a subset of my data:
> species            pair            ratio            dissimilarity
> a            a+b            1.090935            1.870297012
> a            a+c            1.182509            0.691033781
> a            a+d            1.505538            1.441237522
> a            a+e            1.547295            0.953060747
> a            a+f            1.463782            1.306913498
> a            a+g            1.197587            1.331087471
> a            a+h            1.113263            1.097840225
> b            b+a            0.899969            1.870297012
> b            b+c            1.102478            1.548604304
> b            b+d            1.218110            1.669409077
> b            b+e            1.095748            1.536191709
> b            b+f            1.306822            1.579788658
> b            b+g            1.299480            1.084382658
> b            b+h            1.219945            1.137927922
> c            c+a            1.092199            1.441237522
> c            c+b            1.486702            1.669409077
> c            c+d            0.847517            1.688612802
> c            c+e            0.210150            0.651183878
> c            c+f            1.459219            1.064428069
> c            c+g            0.87810            0.590191888
> c            c+h            0.91223            1.37455314
> d            d+a            1.32486            0.953060747
> d            d+b            1.37737            1.536191709
> d            d+c            1.23869            1.310607287
> d            d+e            1.15714            0.651183878
> d            d+f            0.97390            1.22540051
> d            d+g            0.92355            0.640371057
> d            d+h            0.79097            1.224534999
>
> The output from my model, using the whole data set is:
>> summary(random.model)
> Linear mixed model fit by REML
> Formula: ratio~ 1 + (dissimilarity | species)
>   Data: k
>   AIC   BIC  logLik deviance REMLdev
>  11.29 21.42 -0.6465   -3.273   1.293
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr
>  species  (Intercept) 0.092206 0.30365
>          dissimilarity        0.030348 0.17421  -1.000
>  Residual             0.048348 0.21988
> Number of obs: 56, groups: species, 8
>
> Fixed effects:
>            Estimate Std. Error t value
> (Intercept)  1.18024 0.04087   28.88
>
> ----
>
> My questions:
>
> 1. Is my approach correct in plotting the independent variable “dissimilarity” as random intercept, as in ~ 1 + (dissimilarity|species)?

Probably not.  The random effects are defined to have mean 0 so any
non-zero mean must occur in the fixed effects.  Generally the model to
be fit is

ratio ~ 1 + dissumilarity + (1 + dissimilarlity|species)

which is equivalent to

ratio ~ dissimilarity + (dissimilarity|specties)

It is a matter of taste whether to include the 1+ or not in a model
formula like this.  I find that doing so emphasizes to me where each
of the parameters come from.

> 2. On a publication, can I report the variance component of the random terms (in percentage) as the main result of the test?
> 3. Would it be possible to run McMC on a model that does not have Fixed Effects?
> 4. Would there be any other metrics that I should report as well?
>
> I am sure that my questions are pretty basic, but I would strongly appreciate any input from you. Thank you!
> Andre
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models