Est-il possible de spécifier un modèle lmer sans aucun effet fixe?

9

Dans R, comment spécifier le modèle lmer sans effet fixe global? Par exemple, si je dis quelque chose comme

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

le modèle équipé sera

yij=a+αi+βixij

Comment puis-je adapter le modèle

yij=αi+βixij ?

Leo
la source
Ma réponse a été que lmer(y~0+(1|group)+(0+x|group))cela fonctionnerait, mais cela donne une erreur.
Mike Lawrence
4
Je ne pense pas qu'il soit possible de spécifier un modèle sans effet fixe avec lmer car le package lme4 est dédié aux modèles mixtes uniquement (avec au moins un effet fixe et un effet aléatoire). Je ne me souviens pas avoir vu des modèles d'effets aléatoires dans la documentation, en tout cas. Il pourrait être utile de le demander dans la liste des modèles mixtes R-sig ( stat.ethz.ch/mailman/listinfo/r-sig-mixed-models )
maxTC
1
Par simple curiosité, pourquoi voudriez-vous le faire? Il y a peut-être une autre façon d'approcher votre objectif sous-jacent.
jbowman
2
center yfirst :)
Macro

Réponses:

8

Comme @Mike Lawrence l'a mentionné, la chose évidente à faire lors de la définition d'un modèle sans effets fixes est quelque chose sous la forme de:

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

ce qui est en fait assez simple; on ne définit pas d'interception ou une matrice X. La raison fondamentale pour laquelle cela ne fonctionne pas est que, comme l'a souligné @maxTC, "le package lme4 est dédié aux modèles mixtes uniquement ".

En particulier, l'ajustement lmer () calcule la déviance profilée en résolvant la régression du moindre carré pénalisé entre les et ainsi que les effets aléatoires sphériques et (Eq. (11), réf. (2)). Calculée, cette procédure d'optimisation calcule la décomposition de Cholesky du système correspondant en exploitant la structure de blocs du système (Eq. (5), Réf. (1)). La définition d'aucun effet fixe global ne déforme pratiquement cette structure de bloc d'une manière que le code de lmer () ne peut pas supporter. Entre autres choses, la valeur conditionnelle attendue de est basée sur celle de mais la résolution dey^yu0uβ^β^demande la solution d'un système matriciel qui n'a jamais existé (la matrice dans la réf. (1), ou dans la réf. (2)). Vous obtenez donc une erreur comme:RXXLX

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

car après tout, il n'y avait rien à résoudre en premier lieu.

En supposant que vous ne voulez pas réécrire la fonction de coût de déviance profilée lmer (), la solution la plus simple est basée sur l'axiome CS-101: garbage in, garbage out .

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

Donc, ce que nous faisons, c'est définir une variable qui est juste du bruit; comme auparavant, lmer () est chargé de ne pas utiliser d'interception fixe, mais seule la matrice X nous a définis (dans ce cas, la matrice à colonne unique Garbage). Cette variable de bruit gaussien supplémentaire sera en attente non corrélée à nos erreurs de mesure d'échantillon ainsi qu'à votre variance d'effets aléatoires. Inutile de dire que plus votre modèle est structuré, plus la probabilité d'obtenir des corrélations aléatoires indésirables mais statistiquement significatives est faible.Garbage

Donc, lmer () a une variable placebo (matrice) pour jouer avec, dans l'attente, vous obtiendrez la associée à zéro et vous n'avez pas eu à normaliser vos données de quelque façon (les centrer, les blanchir, etc.) . Essayer probablement une initialisation aléatoire de la matrice placebo ne fera pas de mal non plus. Une dernière note pour le "Garbage": l'utilisation du bruit gaussien n'était pas "accidentelle"; il a la plus grande entropie parmi toutes les variables aléatoires de variance égale, donc le moins de chance de fournir un gain d'information.XβX

Clairement, c'est plus une astuce de calcul qu'une solution, mais cela permet à l'utilisateur de spécifier efficacement un modèle lmer sans effet fixe global. Toutes mes excuses pour avoir espéré les deux références. En général, je pense que la référence (1) est le meilleur pari pour quiconque de réaliser ce que fait lmer (), mais la référence (2) est plus proche de l'esprit du code actuel.

Voici un peu de code illustrant l'idée ci-dessus:

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!

Références:

  1. Modèles mixtes linéaires et moindres carrés pénalisés par DM Bates et S. DebRoy, Journal of Multivariate Analysis, Volume 91 Numéro 1, octobre 2004 ( Lien vers la préimpression )
  2. Méthodes de calcul pour les modèles mixtes par Douglas Bates, juin 2012. ( Lien vers la source )
usεr11852
la source
Y a-t-il une raison de préférer l'approche que vous avez proposée à celle mentionnée par Marco dans les commentaires?
russellpierce
2
Oui; vous ne modifiez pas les données, donc aucune transformation arrière n'est nécessaire. En conséquence, tous les diagnostics et goodies standard par lmer (), par exemple. les variables ajustées, les résidus, les niveaux d'effets aléatoires, etc., etc. sont directement interprétables car ils correspondent à votre "vrai" ensemble de données et non à un "altéré".
usεr11852