[R] (no subject)

Daniel Ezra Johnson johnson4 at babel.ling.upenn.edu
Sun Dec 31 18:27:41 CET 2006


> > If one compares the random effect estimates, in fact, one sees that
> > they are in the correct proportion, with the expected signs. They are
> > just approximately eight orders of magnitude too small. Is this a bug?
>
> BLUPs are essentially shrinkage estimates, where shrinkage is
> determined with magnitude of variance. Lower variance more
> shrinkage towards the mean - zero in this case. So this is not a bug.
>
> Gregor

I doubled each data set by duplicating each Subject. There are now
46 subjects in each group instead of 23. So now A and B differ in 2
observations out of 322. The lmer resuls are sensible.

I know BLUPs are not like regular parameters, but doubling (or cutting
in half) the data should not, in my opinion, cause this behavior. There
is a lot of room for the original A variance estimate to be lower than B,
maybe it should be 0.05, 0.01, 0.05, or whatever, but not < .0000000005.

"DOUBLED" (Laplace)
A
Groups  Name        Variance Std.Dev.
Subject (Intercept) 1.759    1.326
Item    (Intercept) 0.178    0.422
number of obs: 322, groups: Subject, 46; Item, 7
B
Groups  Name        Variance Std.Dev.
Subject (Intercept) 1.892    1.376
Item    (Intercept) 0.319    0.564
number of obs: 322, groups: Subject, 46; Item, 7

"ORIGINAL" (Laplace)
A
Groups Name Variance Std.Dev.
Subject (Intercept) 1.63    1.28
Item (Intercept) 5.00e-10 2.24e-05
number of obs: 161, groups: Subject, 23; Item, 7
B
Groups  Name        Variance Std.Dev.
Subject (Intercept) 1.712    1.308
Item    (Intercept) 0.109    0.330
number of obs: 161, groups: Subject, 23; Item, 7

By the way, using the PQL method instead of Laplace "fails"
even more badly with the original data (and gives very
different estimates with the doubled data).

"DOUBLED" (PQL)
A
Groups  Name        Variance Std.Dev.
Subject (Intercept) 2.997    1.731
Item    (Intercept) 0.509    0.713
number of obs: 322, groups: Subject, 46; Item, 7
B
Subject (Intercept) 3.317    1.821
Item    (Intercept) 0.725    0.852
number of obs: 322, groups: Subject, 46; Item, 7
"ORIGINAL" (PQL)
A
1: Estimated variance for factor Item
  is effectively zero
in: LMEopt(x = mer, value = cv)
2: Estimated variance for factors Subject, Item
  is effectively zero
in: LMEopt(x = mer, value = cv)
B
Estimated variance for factors Subject, Item
  is effectively zero
in: LMEopt(x = mer, value = cv)

I understand that the between-Item variance is low, and probably
it is no greater than what you would expect to occur by chance,
but isn't that what hypothesis testing is for (anova, etc.)?

Is my best way around the algorithm returning zero to do
what I have done above, with my real data? That is, duplicate
(or triplicate) Subjects to increase my data size, and thereby
get a reasonably comparable (if wrong) estimate of the Item variance?
Zero is not a reasonable estimate in any of these data sets.

Thanks,
Daniel

> A.fit # "DOUBLED" DATA A
Generalized linear mixed model fit using Laplace
Formula: Response ~ (1 | Subject) + (1 | Item)
    Data: A
Family: binomial(logit link)
AIC BIC logLik deviance
232 243   -113      226
Random effects:
Groups  Name        Variance Std.Dev.
Subject (Intercept) 1.759    1.326
Item    (Intercept) 0.178    0.422
number of obs: 322, groups: Subject, 46; Item, 7

Estimated scale (compare to  1 )  0.81576

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.601      0.327   -7.95  1.9e-15 ***
---

> B.fit  # "DOUBLED" DATA B
Generalized linear mixed model fit using Laplace
Formula: Response ~ (1 | Subject) + (1 | Item)
    Data: B
Family: binomial(logit link)
AIC BIC logLik deviance
237 249   -116      231
Random effects:
Groups  Name        Variance Std.Dev.
Subject (Intercept) 1.892    1.376
Item    (Intercept) 0.319    0.564
number of obs: 322, groups: Subject, 46; Item, 7

Estimated scale (compare to  1 )  0.8052

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    -2.61       0.36   -7.25  4.1e-13 ***

> ranef(A.fit)$Item  # "DOUBLED" DATA A
   (Intercept)
1   -0.084969
2   -0.084969
3   -0.084969
4    0.325780
5   -0.308044
6    0.515590        # 10 1's, 36 0's
7    0.325780
> ranef(B.fit)$Item  # "DOUBLED" DATA B
   (Intercept)
1    -0.11962
2    -0.11962
3    -0.11962
4     0.42389
5    -0.43555
6     0.88322        # 12 1's, 34 0's
7     0.42389

On Dec 31, 2006, at 1:19 AM, Daniel Ezra Johnson wrote:

I am fitting models to the responses to a questionnaire that has seven 
yes/no questions (Item). For each combination of Subject and Item, the 
variable Response is coded as 0 or 1.

I want to include random effects for both Subject and Item. While I 
understand that the datasets are fairly small, and there are a lot of 
invariant subjects, I do not understand something that is happening here, 
and in comparing other subsets of the data.

In the data below, which has been adjusted to show this phenomenon 
clearly, the Subject random effect variance is comparable for A (1.63) and 
B (1.712), but the Item random effect variance comes out as 0.109 for B 
and essentially zero for A (5.00e-10).

Note that the only difference between data set A and data set B occurs on 
row 19, where a single instance of Response is changed.

Item			avg. Response in A		avg. Response in B
1				9%				9%
2				9%				9%
3				9%				9%
4				17%				17%
5				4%				4%
6				22%		<->		26%
7				17%				17%

Why does the Item random effect sometimes "crash" to zero? Surely there is 
some more reasonable estimate of the Item effect than zero. The items 
still have clearly different Response behavior.

If one compares the random effect estimates, in fact, one sees that they 
are in the correct proportion, with the expected signs. They are just 
approximately eight orders of magnitude too small. Is this a bug?

More broadly, is it hopeless to analyze this data in this manner, or else, 
what should I try doing differently? It would be very useful to be able to 
have reliable estimates of random effect sizes, even when they are rather 
small.

I've included replicable code below, sorry that I did not know how to make 
it more compact!

a1 <- c(0,0,0,0,0,0,0)
a2 <- c(0,0,0,0,0,0,0)
a3 <- c(0,0,0,0,0,0,0)
a4 <- c(0,0,0,0,0,0,0)
a5 <- c(0,0,0,0,0,0,0)
a6 <- c(0,0,0,0,0,0,0)
a7 <- c(0,0,0,0,0,0,0)
a8 <- c(0,0,0,0,0,0,0)
a9 <- c(0,0,0,0,0,0,0)
a10 <- c(0,0,0,0,0,0,0)
a11 <- c(0,0,0,0,0,0,0)
a12 <- c(0,0,0,0,0,0,0)
a13 <- c(0,0,0,0,0,0,1)
a14 <- c(0,0,0,0,0,0,1)
a15 <- c(0,0,0,0,0,1,0)
a16 <- c(0,0,0,0,1,0,0)
a17 <- c(0,0,0,1,0,0,0)
a18 <- c(0,0,1,0,0,0,0)
a19 <- c(0,1,0,0,0,0,0)  # differs
a20 <- c(0,1,0,0,0,1,0)
a21 <- c(0,0,0,1,0,1,1)
a22 <- c(1,0,0,1,0,1,1)
a23 <- c(1,0,1,1,0,1,0)
aa <- 
rbind(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23)

b1 <- c(0,0,0,0,0,0,0)
b2 <- c(0,0,0,0,0,0,0)
b3 <- c(0,0,0,0,0,0,0)
b4 <- c(0,0,0,0,0,0,0)
b5 <- c(0,0,0,0,0,0,0)
b6 <- c(0,0,0,0,0,0,0)
b7 <- c(0,0,0,0,0,0,0)
b8 <- c(0,0,0,0,0,0,0)
b9 <- c(0,0,0,0,0,0,0)
b10 <- c(0,0,0,0,0,0,0)
b11 <- c(0,0,0,0,0,0,0)
b12 <- c(0,0,0,0,0,0,0)
b13 <- c(0,0,0,0,0,0,1)
b14 <- c(0,0,0,0,0,0,1)
b15 <- c(0,0,0,0,0,1,0)
b16 <- c(0,0,0,0,1,0,0)
b17 <- c(0,0,0,1,0,0,0)
b18 <- c(0,0,1,0,0,0,0)
b19 <- c(0,1,0,0,0,1,0)  # differs
b20 <- c(0,1,0,0,0,1,0)
b21 <- c(0,0,0,1,0,1,1)
b22 <- c(1,0,0,1,0,1,1)
b23 <- c(1,0,1,1,0,1,0)
bb <- 
rbind(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23)

a <- array(0, c(161,3), list(NULL, c("Subject","Item","Response")))
 	for (s in c(1:23))
 		for (i in c(1:7))
 			a[7*(s-1)+i,] <- c(s,i,aa[s,i])
A <- data.frame(a)

b <- array(0, c(161,3), list(NULL,c("Subject","Item","Response")))
 	for (s in c(1:23))
 		for (i in c(1:7))
 			b[7*(s-1)+i,] <- c(s,i,bb[s,i])
B <- data.frame(b)

A.fit <- lmer(Response~(1|Subject)+(1|Item),A,binomial)
B.fit <- lmer(Response~(1|Subject)+(1|Item),B,binomial)
A.fit
B.fit
ranef(A.fit)$Item
ranef(B.fit)$Item

Generalized linear mixed model fit using Laplace
Formula: Response ~ (1 | Subject) + (1 | Item)
    Data: A
Family: binomial(logit link)
AIC BIC logLik deviance
120 129  -56.8      114
Random effects:
Groups  Name        Variance Std.Dev.
Subject (Intercept) 1.63e+00 1.28e+00
Item    (Intercept) 5.00e-10 2.24e-05
number of obs: 161, groups: Subject, 23; Item, 7

Estimated scale (compare to  1 )  0.83326

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.517      0.395   -6.38  1.8e-10 ***

    Data: B
Family: binomial(logit link)
AIC BIC logLik deviance
123 133  -58.6      117
Random effects:
Groups  Name        Variance Std.Dev.
Subject (Intercept) 1.712    1.308
Item    (Intercept) 0.109    0.330
number of obs: 161, groups: Subject, 23; Item, 7

Estimated scale (compare to  1 )  0.81445

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.498      0.415   -6.02  1.8e-09 ***
---
> ranef(A.fit)$Item
   (Intercept)
1 -2.8011e-10
2 -2.8011e-10
3 -2.8011e-10
4  7.1989e-10
5 -7.8011e-10
6  1.2199e-09    # 5 1's, 18 0's
7  7.1989e-10
> ranef(B.fit)$Item
   (Intercept)
1   -0.056937
2   -0.056937
3   -0.056937
4    0.120293
5   -0.146925
6    0.293893    # 6 1's, 17 0's
7    0.120293



More information about the R-help mailing list