[R-sig-ME] non-nested effects in lme
Andrea Onofri
andrea.onofri at unipg.it
Wed Nov 20 09:08:52 CET 2013
Dear Andrzej,
the use of the auxiliary variables makes the trick! Brilliant
solution! Very many thanks for this.
Regards
Andrea
----------
Andrea Onofri
Department of Agroenvironmental and Crop Sciences
University of Perugia
Italy
On 19 November 2013 22:50, Andrzej Galecki <agalecki at umich.edu> wrote:
> Hi Andrea,
>
> In Chapter 19 of Galecki, Burzykowski recent book titled "Linear
> Mixed-Effects Models Using R: Step-by-Step Approach" you will find an
> example of
> using lme() function for amode with crossed random effects. Related syntax
> (see below) is distributed with nlmeU package available from cran. Note
> that execution time can be fairly long.
>
> Best,
>
> Andrzej
>
> PS. Note the use of auxiliary variables named one1 and one2.
>
> library(lattice)
> data(fcat, package = "nlmeU")
>
> ## ---->>>> NOTE: Code used in Panels R19.1 - R19.7 is stored in Ch19mer.R
> file
>
>
> ###################################################
> ### code chunk: R19.8
> ###################################################
> library(nlme)
>
> fcat1 <- within(fcat, one1 <- one2 <- 1L)
> system.time(
> fm19.2 <-
> lme(scorec ~ 1,
> random = list(
> one1 = pdIdent(~target - 1),
> one2 = pdIdent(~id - 1)),
> data = fcat1))
> fm19.2 # M19.2: (19.2)
>
>
>
>
> On Tue, Nov 19, 2013 at 8:46 AM, Friso Muijsers
> <Friso.muijsers at uni-oldenburg.de> wrote:
>>
>> Ok, I'm sorry for directing you into the wrong direction. Probably someone
>> else has an idea?
>>
>> Am 11/19/2013 2:28 PM, schrieb Andrea Onofri:
>>
>>> Hello Friso,
>>>
>>> very many thanks for your answer. I think that the coding in
>>> StatsExchange is equivalent to:
>>>
>>> lme(Y~ 1, random=~1|A/B, data=X, weights=varIdent(form=~1|A))
>>>
>>> which is actually different from what I am looking for.
>>> Indeed:
>>>
>>> Y <- c(1.6, 2.3, 2.25, 3, 1.6, 2.35, 1.5, 2.85, 1.45, 2.65, 1.95,
>>> 2.65, 1.1, 2.1, 0.7, 2.25, 1.15, 1.65, 0.8, 1.7, 0.95, 1.65,
>>> 0.75, 1.35, 1, 2.05, 0.8, 2, 0.75, 1.9, 0.65, 1.9, 1.4, 2.1,
>>> 1.6, 1.95, 1.05, 1.75, 0.85, 1.75, 1.3, 1.95, 0.95, 1.55, 1,
>>> 1.1, 0.65, 1.05, 1.3, 1.45, 1.05, 0.9, 0.8, 0.9, 0.65, 0.7)
>>> A <- structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L,
>>> 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L,
>>> 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 11L,
>>> 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L
>>> ), .Label = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J",
>>> "K", "L", "M", "N"), class = "factor")
>>> B <- structure(c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L,
>>> 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L,
>>> 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L,
>>> 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("A",
>>> "B", "C", "D"), class = "factor")
>>>
>>> mod <- lmer(Y ~ 1 + (1|A) + (1|B))
>>>
>>> Results in:
>>> ......
>>> Linear mixed model fit by REML
>>> Formula: Y ~ 1 + (1 | A) + (1 | B)
>>> ....
>>> Random effects:
>>> Groups Name Variance Std.Dev.
>>> A (Intercept) 0.180203 0.42450
>>> B (Intercept) 0.163587 0.40446
>>> Residual 0.088169 0.29693
>>> Number of obs: 56, groups: A, 14; B, 4
>>>
>>> Fixed effects:
>>> Estimate Std. Error t value
>>> (Intercept) 1.4839 0.2352 6.308
>>>
>>> While:
>>>
>>> mod2 <- lme(Y ~ 1, random=list(A=~1, B=~1))
>>>
>>> Results in:
>>>
>>> Linear mixed-effects model fit by REML
>>> .......
>>> Random effects:
>>> Formula: ~1 | A
>>> (Intercept)
>>> StdDev: 0.3732374
>>>
>>> Formula: ~1 | B %in% A
>>> (Intercept) Residual
>>> StdDev: 0.4641765 0.1905162
>>>
>>> Fixed effects: Y ~ 1
>>> Value Std.Error DF t-value p-value
>>> (Intercept) 1.483929 0.1201919 42 12.34633 0
>>>
>>> It looks quite different, but perhaps I am missing something? Thank
>>> you again very much
>>>
>>> Andrea Onofri
>>> Department of Agroenvironmental and Crop Sciences
>>> University of Perugia
>>> Italy
>>>
>>>
>>> On 19 November 2013 10:48, Friso Muijsers
>>> <Friso.muijsers at uni-oldenburg.de> wrote:
>>>>
>>>> Am 11/19/2013 9:39 AM, schrieb andrea.onofri at unipg.it:
>>>>>
>>>>> Dear all,
>>>>>
>>>>> I am trying to fit a simple model, relating to a randomised block
>>>>> design where both blocks (A) and treatments (B) are random effects. Coding
>>>>> in lmer, this model would be:
>>>>>
>>>>> model <- lmer(Y ~ 1 + (1|A) + (1|B))
>>>>>
>>>>> However, I would also like to be able to 'manipulate' the correlation
>>>>> structure and thus I assume I have to revert to the lme function in the nlme
>>>>> package. In other cases I have been able to fit non-nested effects in lme by
>>>>> appropriately using the pdMat construct, but, after several efforts, I do
>>>>> not seem to succeed in this simple case. I would greatly appreciate any
>>>>> hints that puts me in the right direction. I thank you very much in advance.
>>>>>
>>>>> Regards
>>>>>
>>>>> Andrea Onofri
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models at r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>
>>>> Hello,
>>>>
>>>> if I understand correctly, you want to specifiy multiple but non-nested
>>>> random effects? If so, I recently found a post on statsExchange where
>>>> someone has found a solution. I didn't check it though, but perhaps it
>>>> fits your needs:
>>>>
>>>> |fit<- lme(Y~ time, random=list(year=~1, date=~time), data=X,
>>>> weights=varIdent(form=~1|year))
>>>>
>>>> or adapted to your general lmer example
>>>>
>>>> ||fit<- lme(Y~ 1, random=list(A=~1, B=~1), data=X,
>>>> weights=varIdent(form=~1|A))||
>>>>
>>>> You can read the details here:
>>>>
>>>>
>>>> http://stats.stackexchange.com/questions/58669/specifying-multiple-separate-random-effects-in-lme
>>>>
>>>> Perhaps that works (didn't check it myself, yet)
>>>>
>>>> Greetings
>>>> |
>>>>
>>>> --
>>>> Friso Muijsers
>>>>
>>>> Institute for Chemistry and Biology of the Marine Environment (ICBM)
>>>> Carl-von-Ossietzky University Oldenburg
>>>> Schleusenstrasse 1
>>>> 26382 Wilhemshaven
>>>>
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>
>> --
>> Friso Muijsers
>>
>> Institute for Chemistry and Biology of the Marine Environment (ICBM)
>> Carl-von-Ossietzky University Oldenburg
>> Schleusenstrasse 1
>> 26382 Wilhemshaven
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
More information about the R-sig-mixed-models
mailing list