StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

重み付き回帰で2種類の重みを使いたい時

Rで重み付きの回帰を実行する場合、lmとかlmer(ランダム効果を使う場合)のweights引数に重みを指定すればいいだけですが、2種類の重みを同時に使いたいような場合があります。たとえばパネル調査データで、母集団からのサンプリングのバイアスを補正するための重みと、初回調査から脱落しやすい属性の人たちが過小評価されてしまうのを補正するための重みとか。
考え方はいろいろあると思いますが、とりあえずは2つの重みの積をweightsに指定しておけばいいんじゃないかと思います。
というのも、下記のコードで確認しているように、重み付き回帰を実行するとはすなわち、重みの平方根を元の変数にかけてからふつうに回帰をやるのと同じ(切片項や、ランダム効果の分散にも重みを反映する必要がある点に注意)なので。いったん変数変換で重みを反映したデータにべつの重みを上からかけ直すのと、重み同士をあらかじめ乗じてからかけ直すのは同じ意味になります。
まぁでも、重みの合成の仕方も、データの性質によって考え方がいろいろありそうな気はしますが(重みの平均をとるとか)、勉強してないのでわかりません。積を取る方法だと、重みが極端な値になって扱いづらいことがあるようです。


以下のコードは、weights引数を使う場合と、あらかじめ変数に重みの平方根をかけてから回帰する場合の比較です。(2種類の重みの反映の例ではないです。)

> library(dplyr)
> library(lme4)
> 
> ### lm関数での重み付き回帰の原理の確認
> 
> # 変数を作成する
> set.seed(1)
> x1 <- rnorm(mean=0, sd=3, n=100)
> x2 <- rnorm(mean=3, sd=1, n=100)
> y <- 1 + 2*x1 + 3*x2 + rnorm(mean=0, sd=3, n=100)
> 
> # 重みをつくる
> w <- rexp(n=100)  # 重み変数
> 
> # 重みの平方根
> root_w <- sqrt(w)
> 
> # 重みなしの回帰
> lm(y ~ x1 + x2) %>%
+   coef()
(Intercept)          x1          x2 
   1.557262    2.021110    2.839600 
> 
> # 重み付き回帰(重み付き最小二乗法)
> lm(y ~ x1 + x2, weights=w) %>%
+   coef()
(Intercept)          x1          x2 
   2.839395    1.938573    2.345262 
> 
> # 重みの平方根をかけてから回帰
> yw <- y * root_w
> x1w <- x1 * root_w
> x2w <- x2 * root_w
> 
> # 切片を0で消して、重みの平方根を切片とする
> lm(yw ~ 0 + root_w + x1w + x2w) %>%
+   coef()
  root_w      x1w      x2w 
2.839395 1.938573 2.345262 
> 
> ### lmer関数での重み付き回帰の原理の確認
> 
> # 変数を作成する
> set.seed(1)
> group_id <- factor(rep(1:100, 10))
> group_effects <- rep(rnorm(n=100, mean=0, sd=3),10)
> x <- rnorm(mean=1, sd=5, n=1000)
> y <- 2 + 3*x + group_effects + rnorm(mean=0, sd=3, n=1000)
> 
> # 重みと重みの平方根
> w <- rexp(n=1000)  # 重み変数
> root_w <- sqrt(w)
> 
> # 重みなしの回帰
> lmer(y ~ x + (1 | group_id)) %>%
+   fixef()
(Intercept)           x 
   2.230189    3.009346 
> 
> # 変数に重みの平方根をかける
> yw <- y * root_w
> xw <- x * root_w
> 
> # モデル1:weights引数を使用
> lmer(y ~ x + (1 | group_id), weights=w) %>%
+   fixef()
(Intercept)           x 
   2.033027    2.998279 
> 
> # モデル2:重みの平方根をかけた変数でふつうに回帰
> # まず、切片を0でいたん消去して、root_wを新たな切片として取り入れる
> # ランダム効果については、分散を重みによってスケーリングする。
> # ランダム効果の記法のうち|の前の部分はもともと、分散を表しているので、
> # そこに重みの平方根を入れればよい。
> # グループidにかけるわけではない。
> lmer(yw ~ 0 + root_w + xw + (root_w | group_id)) %>%
+   fixef()
  root_w       xw 
2.033027 2.998279 


ちなみに、lmにweightsを指定した場合は重み付き最小二乗法をやっていて、重み付き二乗誤差を最小化するようになってるはずですが、lmerはデフォルトでは制限付き最尤法(REML)なのでちょっとややこしいです。