2017-04-05 3 views
1

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 YModè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?

+1

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

+1

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

+1

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

Répondre

0

Comme il est indiqué dans le plus populaire (mais non acceptée) réponse https://stats.stackexchange.com/questions/3143/linear-model-with-constraints, ce problème est facilement résolu en créant une nouvelle variable qui est la différence dans le « one-hot » facteurs codés.

En Stata, on peut le faire comme suit:

* one-hot encode each of the factors 
qui tab factor1, gen(f1dum) 
qui tab factor2, gen(f2dum) 

* generate difference in one-hot vectors 
forv x=1/4{ 
    gen fdiffdum`x' = f1dum`x'-f2dum`x' 
} 

* regress y on differenced one-hot vectors 
reg y fdiffdum2 fdiffdum3 fdiffdum4 

Ce qui donne le résultat suivant:

 Source |  SS   df  MS  Number of obs =  15 
-------------+---------------------------------- F(3, 11)  =  0.73 
     Model | 19.5429062   3 6.51430205 Prob > F  = 0.5554 
    Residual | 98.1766922  11 8.92515383 R-squared  = 0.1660 
-------------+---------------------------------- Adj R-squared = -0.0614 
     Total | 117.719598  14 8.40854274 Root MSE  = 2.9875 

------------------------------------------------------------------------------ 
     y |  Coef. Std. Err.  t P>|t|  [95% Conf. Interval] 
-------------+---------------------------------------------------------------- 
    fdiffdum2 | 2.104393 1.439085  1.46 0.172 -1.063011 5.271798 
    fdiffdum3 | .5222648 1.377463  0.38 0.712 -2.509511  3.55404 
    fdiffdum4 | .6589209 1.266188  0.52 0.613 -2.127941 3.445783 
     _cons | .5054862 .829675  0.61 0.555 -1.320616 2.331589 
------------------------------------------------------------------------------ 

Dans R, on peut le faire comme suit:

factor1mat <- model.matrix(~factor1, DF) 
factor2mat <- model.matrix(~factor2, DF) 

factordiffmat <- factor1mat - factor2mat 

summary(lm(Y~factordiffmat, data=DF)) 

Coefficients: (1 not defined because of singularities) 
         Estimate Std. Error t value Pr(>|t|) 
(Intercept)    0.5055  0.8297 0.609 0.555 
factordiffmat(Intercept)  NA   NA  NA  NA 
factordiffmatfactor1B  2.1044  1.4391 1.462 0.172 
factordiffmatfactor1C  0.5223  1.3775 0.379 0.712 
factordiffmatfactor1D  0.6589  1.2662 0.520 0.613