¿Es posible especificar un modelo lmer sin efectos fijos?

9

En R, ¿cómo especifico el modelo lmer sin efecto fijo global? Por ejemplo, si digo algo como

lmer(y ~ (1 | group) + (0 + x | group), data = my_df)

el modelo ajustado será

yij=a+αi+βixij

¿Cómo encajo modelo?

yij=αi+βixij ?

León
fuente
Mi respuesta fue que lmer(y~0+(1|group)+(0+x|group))funcionaría, pero esto arroja un error.
Mike Lawrence
44
No creo que sea posible especificar un modelo sin un efecto fijo con lmer porque el paquete lme4 está dedicado solo a modelos mixtos (con al menos un efecto fijo y un efecto aleatorio). No recuerdo haber visto modelos de efectos aleatorios en la documentación, en cualquier caso. Podría ser útil solicitarlo en la lista de modelos mixtos R-sig ( stat.ethz.ch/mailman/listinfo/r-sig-mixed-models )
maxTC
1
Solo por curiosidad, ¿por qué querrías hacerlo? Quizás haya otra forma de acercarse a su objetivo subyacente.
jbowman
2
centrar yprimero :)
Macro

Respuestas:

8

Como @Mike Lawrence mencionó que lo más obvio al definir un modelo sin efectos fijos es algo en forma de:

lmer(y ~ -1 + (1|GroupIndicator))

que en realidad es bastante sencillo; uno no define ninguna intersección o una matriz X. La razón básica por la que esto no funciona es que, como señaló @maxTC, "el paquete lme4 está dedicado solo a modelos mixtos ".

En particular, lo que hace el ajuste lmer () es calcular la desviación perfilada resolviendo la regresión de mínimos cuadrados penalizada entre e , así como los efectos aleatorios esféricos y (Ec. (11), Ref. (2)). Computacionalmente, este procedimiento de optimización calcula la descomposición de Cholesky del sistema correspondiente explotando la estructura de bloques del sistema (Ec. (5), Ref. (1)). Al no establecer efectos fijos globales, prácticamente se distorsiona esa estructura de bloques de una manera que el código de lmer () no puede hacer frente. Entre otras cosas, el valor condicional esperado de se basa en pero resuelve paray^ytu0 0tuβ^β^pregunta la solución de un sistema matricial que nunca existió (la matriz en la Ref. (1), o en la Ref. (2)). Entonces obtienes un error como:RXXLX

Error in mer_finalize(ans) : 
  Cholmod error 'invalid xtype' at file:../Cholesky/cholmod_solve.c, line 970

porque, después de todo, no había nada que resolver en primer lugar.

Suponiendo que no desea volver a escribir la función de costo de desviación perfilada lmer (), la solución más fácil se basa en el axioma CS-101: basura adentro, basura afuera .

 N = length(y); Garbage <- rnorm(N); 
 lmer(y ~ -1 + Garbage + (1|GroupIndicator));

Entonces, lo que hacemos es definir una variable que sea solo ruido; como antes, a lmer () se le indica que no use una intersección fija, sino que solo la matriz X nos definió (en este caso, la matriz de columna única Garbage). Esta variable de ruido gaussiana adicional no estará correlacionada con nuestros errores de medición de muestra, así como con su variación de efectos aleatorios. No es necesario decir que cuanto más estructura tenga su modelo, menor será la probabilidad de obtener correlaciones aleatorias no deseadas pero estadísticamente significativas.solunarsiunasolmi

Así LMER () tiene un placebo variable (matriz) para jugar, a la espera que se obtiene el correspondiente sea cero y que no tenías para normalizar sus datos en cualquier forma (centrado ellos, para blanquear los etc.) . Probablemente intentar un par de inicializaciones aleatorias de la matriz placebo tampoco le hará daño. Una nota final para la "Basura": el uso del ruido gaussiano no fue "accidental"; tiene la mayor entropía entre todas las variables aleatorias de igual varianza, por lo que tiene menos posibilidades de proporcionar una ganancia de información.XβX

Claramente, esto es más un truco computacional que una solución, pero permite al usuario especificar efectivamente un modelo más antiguo sin un efecto fijo global. Disculpas por esperar alrededor de las dos referencias. En general, creo que la Ref. (1) es la mejor apuesta para que alguien se dé cuenta de lo que está haciendo lmer (), pero la Ref. (2) está más cerca del espíritu del código real.

Aquí hay un poco de código que muestra la idea anterior:

library(lme4)
N= 500;                 #Number of Samples
nlevA = 25;             #Number of levels in the random effect 
set.seed(0)             #Set the seed
e = rnorm(N); e = 1*(e - mean(e) )/sd(e); #Some errors

GroupIndicator =  sample(nlevA, N, replace=T)   #Random Nvel Classes 

Q = lmer( rnorm(N) ~ (1| GroupIndicator ));      #Dummy regression to get the matrix Zt easily
Z = t(Q@Zt);            #Z matrix

RA <-  rnorm(nlevA )                        #Random Normal Matrix
gammas =c(3*RA/sd(RA))                      #Colour this a bit

y = as.vector(  Z %*% gammas +  e )         #Our measurements are the sum of measurement error (e) and Group specific variance

lmer_native <- lmer(y ~ -1 +(1| GroupIndicator)) #No luck here.
Garbage <- rnorm(N)                              #Prepare the garbage

lmer_fooled <- lmer(y ~ -1 + Garbage+(1| GroupIndicator)) #OK...
summary(lmer_fooled)                             #Hey, it sort of works!

Referencias

  1. Modelos lineales mixtos y mínimos cuadrados penalizados por DM Bates y S. DebRoy, Journal of Multivariate Analysis, Volumen 91 Número 1, octubre de 2004 ( Enlace a la preimpresión )
  2. Métodos computacionales para modelos mixtos por Douglas Bates, junio de 2012. ( Enlace a la fuente )
usεr11852
fuente
¿Hay alguna razón para preferir el enfoque que ha presentado sobre el mencionado por Marco en los comentarios?
russellpierce
2
Si; no altera los datos, por lo que no se necesitan transformaciones posteriores. Como resultado, todos los diagnósticos y productos estándar de lmer (), p. Ej. Las variables ajustadas, los residuos, los niveles de efectos aleatorios, etc., etc., son directamente interpretables, ya que corresponden a su conjunto de datos "verdadero" y no a uno "alterado".
usεr11852