Cette question est en rapport avec https://stats.stackexchange.com/questions/3143/linear-model-with-constraints, mais un scénario légèrement différent. J'ai un modèle linéaire simple à 2 facteurs avec un résultat continu Y
Modèle linéaire avec contraintes d'égalité entre deux variables de facteur
factor1
a ~ 350 valeurs catégoriques, et factor2
a les mêmes ~ 350 catégories. Je veux contraindre le coefficient sur chaque niveau à la somme de zéro à travers les deux facteurs.
(La raison est que chaque niveau de factor1
et factor2
entre soit dans un exemple de formation positivement ou négativement, mais jamais apparaît deux fois dans le même exemple.)
Voici un ensemble de données exemple illustrant la situation, où il y a quatre niveaux de chaque facteur:
Y factor1 factor2
1 -1.2470416 A B
2 4.3368592 C D
3 1.0005147 D A
4 -2.8309146 A C
5 1.7501315 B D
6 -0.8372193 B A
7 3.3542627 C A
8 4.3319422 D C
9 1.4937895 D B
10 2.0951559 A D
11 -2.6610207 C D
12 -4.9917367 D B
13 2.2424169 D A
14 1.0205409 C A
15 -3.4584576 C B
Le modèle statistique que je veux estimer est: $$ y _ {(i, j)} = \ alpha_i- \ beta_j + \ varepsilon _ {(i, j) } $$ où $ (i, j) $ est un résultat qui dépend de la paire. factor1
marques $ i $ et factor2
marques $ j $. Si le groupe A
apparaît dans factor2
, le paramètre sur A
doit être égal au négatif de s'il est apparu dans factor1
. Ainsi, je voudrais mettre $ \ alpha $ égal à $ \ beta $ pour tout $ i $ et $ j $.
je peux estimer une version (absurde) de ce modèle en lm()
assez facilement comme suit:
Y <- c(-1.2470416, 4.3368592 , 1.0005147 , -2.8309146 , 1.7501315 , -0.8372193 , 3.3542627 , 4.3319422 , 1.4937895 , 2.0951559 , -2.6610 207 , -4.9917367 , 2.2424169 , 1.0205409 , -3.4584576)
factor1 <- c("A" , "C" , "D" , "A" , "B" , "B" , "C" , "D" , "D" , "A" , "C" , "D" , "D" , "C" , "C")
factor2 <- c("B", "D", "A", "C", "D", "A", "A", "C", "B", "D", "D", "B", "A", "A", "B")
DF <- data.frame(Y,factor1,factor2)
lm(Y~factor1+factor2,data=DF)
et je reçois la sortie suivante:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5363 2.5856 0.207 0.841
factor1B -0.4579 3.1121 -0.147 0.887
factor1C 0.4047 2.4925 0.162 0.875
factor1D 1.8737 2.4098 0.778 0.459
factor2B -3.6252 2.2050 -1.644 0.139
factor2C -0.7226 2.8903 -0.250 0.809
factor2D 0.7561 2.2094 0.342 0.741
Notez que, théoriquement, factor1C
devrait -factor2C
comme dicté par mon modèle. Ce n'est pas le cas dans la sortie simple lm()
car je n'ai imposé aucune contrainte.
Donc ce que je voudrais faire est d'estimer
Y ~ factor1 + factor2 [subject to factor1+factor2=0 for each level of factor1, factor2]
En clair, ce serait quelque chose comme
model2 <- lm(Y~factor1-factor2, data=DF)
Mais cela est bien sûr pas comment R interprète cette expression (parce que mettre un signe moins dans une instruction model
indique à R d'exclure cette variable du modèle). J'ai lu sur les contrastes, mais je ne pense pas qu'il existe un moyen de le faire. J'ai également lu sur glmc
, mais je n'ai pas vu un moyen simple de l'incorporer pour les facteurs qui ont autant de niveaux. En outre, il n'est pas clair pour moi que la génération d'un nouveau factor3 = factor1-factor2
est une opération bien définie pour ce scénario spécifique. Enfin, j'ai essayé d'exécuter model3 <- lm(Y+factor2 ~ factor1, data=DF)
mais j'ai reçu une erreur. J'ai le sentiment que je devrais créer une matrice de contraintes en faisant une boucle sur les niveaux de chaque variable. Je suis suffisamment nouveau à R pour ne pas savoir exactement comment cela se fait. Toute aide serait appréciée.
Notez qu'il est assez facile de le faire dans Stata, comme suit:
input ID y factor1 factor2
1 -1.2470416 1 2
2 4.3368592 3 4
3 1.0005147 4 1
4 -2.8309146 1 3
5 1.7501315 2 4
6 -0.8372193 2 1
7 3.3542627 3 1
8 4.3319422 4 3
9 1.4937895 4 2
10 2.0951559 1 4
11 -2.6610207 3 4
12 -4.9917367 4 2
13 2.2424169 4 1
14 1.0205409 3 1
15 -3.4584576 3 2
end
constraint 1 2.factor1 = -2.factor2
constraint 2 3.factor1 = -3.factor2
constraint 3 4.factor1 = -4.factor2
cnsreg y i.factor1 i.factor2, constraints(1/3)
qui donne le résultat suivant:
Constrained linear regression Number of obs = 15
F( 3, 11) = 0.73
Prob > F = 0.5554
Root MSE = 2.9875
(1) 2.factor1 + 2.factor2 = 0
(2) 3.factor1 + 3.factor2 = 0
(3) 4.factor1 + 4.factor2 = 0
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
factor1 |
B | 2.104393 1.439085 1.46 0.172 -1.063011 5.271798
C | .5222649 1.377463 0.38 0.712 -2.509511 3.55404
D | .6589209 1.266188 0.52 0.613 -2.127941 3.445783
|
factor2 |
B | -2.104393 1.439085 -1.46 0.172 -5.271798 1.063011
C | -.5222649 1.377463 -0.38 0.712 -3.55404 2.509511
D | -.6589209 1.266188 -0.52 0.613 -3.445783 2.127941
|
_cons | .5054862 .829675 0.61 0.555 -1.320616 2.331589
------------------------------------------------------------------------------
Comment peut-on faire ce qui précède en R?
Je ne sais pas si vous voulez contraindre _coefficients_ sur factor1 et factor2 pour faire une somme nulle, ou si les valeurs sont contraintes à somme à 0 ... – MichaelChirico
Ma compréhension de cette question est que 'factor1' et' factor2' sont parfaitement multicolinéaires. Donc, vous ne pouvez inclure l'un ou l'autre ... – MichaelChirico
Je ne pense pas que ce soit vraiment une question sur le code R/R, et en tant que tel, je ne pense pas que cette question est hors sujet à cet égard. OTOH, je ne suis pas vraiment votre situation, ou comment elle motive votre solution suspectée. D'ailleurs, je ne suis pas certain d'être clair sur ce que votre solution * est * (par exemple, je partage les confusions de @ MichaelChirico). Il pourrait être utile de développer un exemple plus simple avec juste quelques niveaux et un exemple de jeu de données pour y aller, et ensuite ajouter quelques explications supplémentaires. – gung