1 データ

n <- 100
x1 <- 1:n
x2 <- sqrt(x1)
y  <- 0.5 * x1 + 1.0 * x2 + 0.1 * x1 * x2 + rnorm(n)
d <- data.frame(y, x1, x2)

2 フィッティング

Rのlm関数のモデル式で\(説明変数1\times 説明変数2\)と記入すると， 説明変数1，説明変数2，説明変数1と説明変数2の交互作用項の３変数に展開される。

fit <- lm(y ~ x1 * x2, data = d)

library(sjPlot)
tab_model(fit)
  y
Predictors Estimates CI p
(Intercept) 0.34 -1.95 – 2.64 0.766
x1 0.60 0.34 – 0.87 <0.001
x2 0.59 -0.83 – 2.01 0.410
x1 * x2 0.09 0.08 – 0.11 <0.001
Observations 100
R2 / R2 adjusted 1.000 / 1.000

3 計画行列

head(model.matrix(y ~ x1 + x2, data = d))
##   (Intercept) x1       x2
## 1           1  1 1.000000
## 2           1  2 1.414214
## 3           1  3 1.732051
## 4           1  4 2.000000
## 5           1  5 2.236068
## 6           1  6 2.449490
head(model.matrix(y ~ x1 * x2, data = d))
##   (Intercept) x1       x2     x1:x2
## 1           1  1 1.000000  1.000000
## 2           1  2 1.414214  2.828427
## 3           1  3 1.732051  5.196152
## 4           1  4 2.000000  8.000000
## 5           1  5 2.236068 11.180340
## 6           1  6 2.449490 14.696938

4 変数自体を変換するときの注意

I(大文字のアイ)関数で，説明変数の変換式を括ること。

head(model.matrix(y ~ x1^2, data = d)) # 無効
##   (Intercept) x1
## 1           1  1
## 2           1  2
## 3           1  3
## 4           1  4
## 5           1  5
## 6           1  6
head(model.matrix(y ~ I(x1^2), data = d))
##   (Intercept) I(x1^2)
## 1           1       1
## 2           1       4
## 3           1       9
## 4           1      16
## 5           1      25
## 6           1      36

モデル式の書式詳細は次のコマンドで確認できる。

?formula