From: Daniel Ezra Johnson <johnson4_at_babel.ling.upenn.edu>

Date: Sun, 31 Dec 2006 04:28:40 -0500 (EST)

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)

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 <- data.frame(a)

B <- data.frame(b)

ranef(A.fit)$Item

ranef(B.fit)$Item

6 1.2199e-09

7 7.1989e-10

6 0.293893

7 0.120293

R-help_at_stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon 01 Jan 2007 - 17:07:26 GMT

Date: Sun, 31 Dec 2006 04:28:40 -0500 (EST)

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. in A (%) avg. 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, so sensitively? Surely there is some more reasonable estimate of the Item effect here 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 ***

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

> 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

R-help_at_stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon 01 Jan 2007 - 17:07:26 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Fri 13 Apr 2007 - 22:32:55 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*