non-nested effects in lme

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

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

--
Friso Muijsers

Institute for Chemistry and Biology of the Marine Environment (ICBM)
Carl-von-Ossietzky University Oldenburg
Schleusenstrasse 1
26382 Wilhemshaven