Qual seria uma imagem ilustrativa para modelos lineares mistos?

Diga que você está na biblioteca do seu departamento de estatística e se depara com um livro com a figura a seguir na primeira página.

Você provavelmente pensará que este é um livro sobre coisas de regressão linear.

Qual seria a imagem que faria você pensar em modelos lineares mistos?

Para uma palestra, usei a figura a seguir, baseada no sleepstudyconjunto de dados do pacote lme4 . A idéia era ilustrar a diferença entre ajustes de regressão independentes de dados específicos de sujeitos (cinza) versus previsões de modelos de efeitos aleatórios, especialmente que (1) os valores previstos do modelo de efeitos aleatórios são estimadores de encolhimento e que (2) trajetórias de indivíduos compartilham uma inclinação comum com um modelo somente de interceptação aleatória (laranja). As distribuições de interceptações de sujeitos são mostradas como estimativas de densidade de kernel no eixo y ( código R ).

_{(As curvas de densidade se estendem além da faixa de valores observados porque há relativamente poucas observações.)}

Um gráfico mais "convencional" pode ser o próximo, que é de Doug Bates (disponível no site R-forge para lme4 , por exemplo, 4Longitudinal.R ), onde podemos adicionar dados individuais em cada painel.

Então, algo não "extremamente elegante", mas mostrando interceptações aleatórias e inclinações também com R. (acho que seria ainda mais legal se mostrassem as equações reais também)

N =100; set.seed(123);


x1 = runif(N)*3; readings1 <- 2*x1 + 1.0 + rnorm(N)*.99;
x2 = runif(N)*3; readings2 <- 3*x2 + 1.5 + rnorm(N)*.99;
x3 = runif(N)*3; readings3 <- 4*x3 + 2.0 + rnorm(N)*.99;
x4 = runif(N)*3; readings4 <- 5*x4 + 2.5 + rnorm(N)*.99;
x5 = runif(N)*3; readings5 <- 6*x5 + 3.0 + rnorm(N)*.99;

X = c(x1,x2,x3,x4,x5);
Y = c(readings1,readings2,readings3,readings4,readings5)
Grouping  = c(rep(1,N),rep(2,N),rep(3,N),rep(4,N),rep(5,N))

library(lme4);
LMERFIT <- lmer(Y ~ 1+ X+ (X|Grouping))

RIaS <-unlist( ranef(LMERFIT)) #Random Intercepts and Slopes
FixedEff <- fixef(LMERFIT)    # Fixed Intercept and Slope

png('SampleLMERFIT_withRandomSlopes_and_Intercepts.png', width=800,height=450,units="px" )
par(mfrow=c(1,2))
plot(X,Y,xlab="x",ylab="readings")
plot(x1,readings1, xlim=c(0,3), ylim=c(min(Y)-1,max(Y)+1), pch=16,xlab="x",ylab="readings" )
points(x2,readings2, col='red', pch=16)
points(x3,readings3, col='green', pch=16)
points(x4,readings4, col='blue', pch=16)
points(x5,readings5, col='orange', pch=16)
abline(v=(seq(-1,4 ,1)), col="lightgray", lty="dotted");        
abline(h=(seq( -1,25 ,1)), col="lightgray", lty="dotted")   

lines(x1,FixedEff[1]+ (RIaS[6] + FixedEff[2])* x1+ RIaS[1], col='black')
lines(x2,FixedEff[1]+ (RIaS[7] + FixedEff[2])* x2+ RIaS[2], col='red')
lines(x3,FixedEff[1]+ (RIaS[8] + FixedEff[2])* x3+ RIaS[3], col='green')
lines(x4,FixedEff[1]+ (RIaS[9] + FixedEff[2])* x4+ RIaS[4], col='blue')
lines(x5,FixedEff[1]+ (RIaS[10]+ FixedEff[2])* x5+ RIaS[5], col='orange') 
legend(0, 24, c("Group1","Group2","Group3","Group4","Group5" ), lty=c(1,1), col=c('black','red', 'green','blue','orange'))
dev.off()

Este gráfico, extraído da documentação do nlmefit do Matlab, me parece um exemplo real do conceito de interceptações e inclinações aleatórias. Provavelmente algo mostrando grupos de heterocedasticidade nos resíduos de um gráfico OLS também seria bastante padrão, mas eu não daria uma "solução".