カテゴリカルデータで制約付の回帰。
結局結論をまとめると
y=a1*x1+a2*x2+a3*x3+Cで、
x=[1,2,3], x1,x2,x3ダミー変数のとき(0,1)、
a<-NULL
a[1]<-100
a[2]<-300
a[3]<-400
aa<-rbinom(size=2,1000,prob=0.5)+1
y<-a[aa]+rnorm(1000)
x<-factor(aa)
lm(y~factor(x),list(x = contr.sum))
結果
Call:
lm(formula = y ~ x, contrasts = list(x = contr.sum))Coefficients:
(Intercept) x1 x2
266.69 -166.68 33.33
とすれば、回帰係数a1:x1とa2:x2が求められる。a3=(0-a1-a2)となる。
一方、
lm(y~factor(x),list(x = contr.treatment))
Call:
lm(formula = y ~ x, contrasts = list(x = contr.treatment))Coefficients:
(Intercept) x2 x3
100 200 300
とすると、a1がベースラインモデルがもとまる、係数u2:x2とu3:x3とすると。
y=0*x1+u2*x1+u3*x3+C
のイメージ。
理由は、contr.sumのコーディングは、
y=a1*x1+a2*x2+a3*x3+C
をa3=-a1-a2としたときの計算、
y=a1*(x1-x3)+a2*(x2-x3)をおこなっているため。
具体的には、コーディングが
u1 u2
x1 1 0
x2 0 1
x3 -1 -1
となっている。
一方、デフォルトのcontr.treatmentの場合は、
u1 u2
x1 0 0
x2 1 0
x3 0 1
なので、x1が0でデフォルト扱いになっている。
ちなみに、以下のようにすれば、コーディングが変数ごとにわけられる。
lm(y~x1+x2,contrasts=list(x1=contr.sum,x2=contr.treatment))
まとめて設定するときは
def_con<-options("contrasts")
#変更、１変数名が順序なし、２変数名が順序ありの場合
options(contrasts=c("contr.sum","contr.sum"))
#ひとしきりの処理
#元に戻す
options(contrasts=def_con$contrasts)
参考（デフォルト値）
$contrasts
unordered ordered
"contr.treatment" "contr.poly"
以下は、例
a<-NULL a[1]<-100 a[2]<-300 a[3]<-400 b<-NULL b[1]<-30 b[2]<-20 aa<-rbinom(size=2,1000,prob=0.5)+1 bb<-rbinom(size=1,1000,prob=0.5)+1 y<-a[aa]+b[bb] a1<-rep(0,length(aa)) a1[which(aa==1)]<-1 a2<-rep(0,length(aa)) a2[which(aa==2)]<-1 a3<-rep(0,length(aa)) a3[which(aa==3)]<-1 b1<-rep(0,length(bb)) b1[which(bb==1)]<-1 b2<-rep(0,length(bb)) b2[which(bb==2)]<-1 lm(y~a1+a2+a3+b1+b2+0) y<-b[bb]+rnorm(1000) lm(y~factor(bb),list(b = contr.sum)) y<-a[aa]+rnorm(1000) x<-factor(aa) lm(y~x,contrasts=list(x = contr.sum)) y2<--166.6*a1+33.3*a2+(0+166.6-33.3)*a3+266.7 y<-a[aa]+b[bb]+rnorm(1000,sd=0.001) x1<-factor(aa) x2<-factor(bb) lm(y~x1+x2,contrasts=list(x1 = contr.sum,x2=contr.sum)) y2<--166.6*a1+33.3*a2+(0+166.6-33.3)*a3+5.00*b1-5.00*b2+291.49 plot(y,y2) curve(1.0*x,add=T) y<-a[aa]+b[bb]+rnorm(1000,sd=0.001) x1<-factor(aa) x2<-factor(bb) lm(y~x1+x2,contrasts=list(x1=contr.treatment,x2=contr.treatment)) y2<--0*a1+200*a2+300*a3+0*b1-10*b2+130 plot(y,y2) curve(1.0*x,add=T) y<-a[aa]+b[bb]+rnorm(1000,sd=0.001) x1<-factor(aa) x2<-factor(bb) lm(y~x1+x2,contrasts=list(x1=contr.sum,x2=contr.treatment)) y2<- -166.67*a1+33.33*a2+(1+166.7-33.33)*a3+0*b1-10*b2+296 plot(y,y2) curve(1.0*x,add=T)
例：
y=x+e
x1:水準１
x2:水準２
x3:水準3
の場合、通常は、
y=a1*x1+a2*x2+e
としてx3の分を消して計算するが、
a1とa2とa3の効果をそれぞれ明確にしたいので、制約付の回帰
y=a1x1+a2x2+a3x3+e
かつ制約条件 a1+a2+a3=0
としたい。
その場合の方法はcontrastsのパラメータをいじれば解決する。
例：通常
#３水準のデータを発生 x<-rbinom(1000,prob=0.5,size=2) #説明変数を発生 y<-2*x+rnorm(length(x)) #回帰 m<-glm(y~x) summary(m)
出力はfactor1とfactor2の係数のみ表示している。
Call:
glm(formula = y ~ factor(x))Deviance Residuals:
Min 1Q Median 3Q Max
- 3.13863 -0.63221 -0.00327 0.62930 3.00819
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.006192 0.066530 -0.093 0.926
factor(x)1 2.059249 0.080387 25.617 <2e-16 ***
factor(x)2 3.978692 0.088639 44.887 <2e-16 ***
- -
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.9914657)
Null deviance: 2994.38 on 999 degrees of freedom
Residual deviance: 988.49 on 997 degrees of freedom
AIC: 2834.3Number of Fisher Scoring iterations: 2
制約版
#３水準のデータを発生 x<-rbinom(1000,prob=0.5,size=2) #説明変数を発生 y<-2*x+rnorm(length(x)) #回帰 m<-glm(y~xb+0, contrasts = list(xb = contr.sum)) summary(m)
出力
summary(m)
Call:
glm(formula = y ~ xb + 0, contrasts = list(xb = contr.sum))Deviance Residuals:
Min 1Q Median 3Q Max
- 3.1360 -0.6501 0.0303 0.6536 2.8382
Coefficients:
Estimate Std. Error t value Pr(>|t|)
xb0 -0.07778 0.06520 -1.193 0.233
xb1 1.93244 0.04422 43.703 <2e-16 ***
xb2 3.90256 0.05740 67.990 <2e-16 ***
- -
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.9521565)
Null deviance: 7170.7 on 1000 degrees of freedom
Residual deviance: 949.3 on 997 degrees of freedom
AIC: 2793.8Number of Fisher Scoring iterations: 2
定数問題。
例２）
参考：
How to fit a glm with sum to zero constraints in R (no reference level)
http://stats.stackexchange.com/questions/162381/how-to-fit-a-glm-with-sum-to-zero-constraints-in-r-no-refere
nce-level
http://stats.stackexchange.com/questions/3143/linear-model-with-constraints
ANOVA君/平方和のタイプ/タイプⅢ平方和の計算法 - 井関龍太のページ http://riseki.php.xdomain.jp/index.php?ANOVA%E5%90%9B%2F%E5%B9%B3%E6%96%B9%E5%92%8C%E3%81%AE%E3%82%B
http://faculty.nps.edu/sebuttre/home/R/contrasts.html
Setting and Keeping Contrasts
http://faculty.nps.edu/sebuttre/home/R/contrasts.html
R aov関数を分解する
http://nakhirot.hatenablog.com/entry/20130717/1374065833
optionsはグローバルな(リセットしない限り常に有効な)オプションを設定する関数。Rのヘルプによると、contrastsはmodel fittingにおいて使われるオプションのことで、前者がunordered factorsに対して使われるもの、後者がordered factorsに対して使われるもの。"contr.helmert"は第2水準を第1水準と比較し、第3水準を第1水準と第2水準の平均と比較…するオプション。"contr.poly"は"直交多公式系"で比較を行うオプション。(詳細要確認)。on.exitはこの関数(aov関数)の実行後にオプションを保存しておく関数。
