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)
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)
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 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 > B.fit
Generalized linear mixed model fit using Laplace
Formula: Response ~ (1 | Subject) + (1 | Item)
    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 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 > 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
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
7    0.120293
	[[alternative HTML version deleted]]

