# Finding stability in a time series of success/failure data — Mixed Model Random Effects

I have a very large data set of success/failure data for hundreds of subjects. The chance of success is rather small but does/should stabilize over time. Subjects have varying levels of observations/trials.

The aim is to find out how many observations we need for a subject until their performance 'stabilizes' and becomes more reliable.

I've previously used split-half reliability to determine how many observations we would need to see stability. The issue with this data set is that there are a lot of 0's so when taking random samples for split half reliability it is entirely possible that the correlation that gets returned is rather high given that people have a lot of zeros in both random samples.

As such, I was thinking to use a mixed model to try and solve this issue, however I'm running into a similar problem with all of the 0's (failures). The approach I was going to use comes from Gelman & Hill's book Data Analysis Using Regression and Multilevel/Hierarchical Models (pg. 268).

Here is a small example of the data:

# Load libraries

library(tidyverse)
library(zoo)
library(lme4)

# Create data (Success = 1 | Failure = 0)

set.seed(4)
subject <- rep(c("A", "B", "C"), times = c(6, 20, 40))
variable <- sample(c(0, 1), size = length(subject), replace = T, prob = c(0.9, 0.1))

df <- data.frame(subject, variable)

df <- df %>%
group_by(subject) %>%
mutate(Trial = seq_along(subject),
Success_Prob = cummean(variable)) %>%
as.data.frame

subject variable Trial Success_Prob
1        A        0     1    0.0000000
2        A        0     2    0.0000000
3        A        0     3    0.0000000
4        A        0     4    0.0000000
5        A        0     5    0.0000000
6        A        0     6    0.0000000
7        B        0     1    0.0000000
8        B        1     2    0.5000000
9        B        1     3    0.6666667
10       B        0     4    0.5000000
11       B        0     5    0.4000000
12       B        0     6    0.3333333
13       B        0     7    0.2857143
14       B        1     8    0.3750000
15       B        0     9    0.3333333

# Build mixed model

fit <- lmer(Success_Prob ~ 1 + (1|subject), data = df)
summary(fit)

Linear mixed model fit by REML ['lmerMod']
Formula: Success_Prob ~ 1 + (1 | subject)
Data: df

REML criterion at convergence: -86.8

Scaled residuals:
Min      1Q  Median      3Q     Max
-3.3715 -0.3627  0.0464  0.4162  3.3207

Random effects:
Groups   Name        Variance Std.Dev.
subject  (Intercept) 0.05108  0.2260
Residual             0.01272  0.1128
Number of obs: 66, groups:  subject, 3

Fixed effects:
Estimate Std. Error t value
(Intercept)   0.6110     0.1318   4.636


Gelman recommends dividing the variance of the intercept random effect by the variance of the residual random effect. Then taking that answer and dividing 1 by it. In the book example, the residual random effect (0.76) is larger than the intercept random effect (0.16). Thus, when calculating the number of observations needed for the county in his model to start to move from a pooled model to a no-pooling estimate the answer comes out to 23. My answer, since the residual is smaller than the intercept random effect doesn't make intuitive sense.

# From Gelman's book, using the random effects variance

1 /(0.2260^2 / 0.1128^2) # 0.249


Does anyone have an idea how they might approach this problem?