R: 順序尺度の説明変数 [統計]
順序尺度変数を説明変数とした線形モデルというのをやってみる。
まずはデータを用意。説明変数はA, B, C, Dの4水準。x0は名義尺度、x1は順序尺度にしてみる。
set.seed(1234) lv <- c("A", "B", "C", "D") rp <- 10 x <- rep(c("A", "B", "C", "D"), each = rp) x0 <- factor(x, ordered = FALSE) x1 <- factor(x, ordered = TRUE) y <- rep(c(2, 2.5, 5, 5.5), each = rp) + rnorm(length(x), 0, 0.5)
データを表示する。
> x0 [1] A A A A A A A A A A B B B B B B B B B B C C C C C C C C C C D D D D D D D D D D Levels: A B C D > x1 [1] A A A A A A A A A A B B B B B B B B B B C C C C C C C C C C D D D D D D D D D D Levels: A < B < C < D > y [1] 1.966021 1.607139 2.076757 2.018520 2.463188 1.958114 2.047564 2.528151 1.939917 2.181881 [11] 1.806544 1.995214 2.772020 2.280906 2.203969 1.730568 2.134171 2.891899 1.452519 2.591360 [21] 4.860152 5.279696 5.646491 5.583011 5.412890 4.996073 4.362742 5.219878 5.123312 4.667861 [31] 5.203847 6.077237 4.557721 6.033058 4.724702 5.129178 5.087102 5.466280 4.164443 5.935812
説明変数それぞれについてのyの平均値を表示する。
> sapply(lv, function(i) mean(y[x == i])) A B C D 2.078725 2.185917 5.115211 5.237938
グラフを表示する。
> plot(y ~ x0, las = 1, xlab = "x", ylab = "y")
線形モデルにあてはめてみる。
> summary(lm(y ~ x0)) Call: lm(formula = y ~ x0) Residuals: Min 1Q Median 3Q Max -1.07349 -0.20679 -0.03263 0.31938 0.83930 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.0787 0.1476 14.085 3.23e-16 *** x0B 0.1072 0.2087 0.514 0.611 x0C 3.0365 0.2087 14.548 < 2e-16 *** x0D 3.1592 0.2087 15.136 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4667 on 36 degrees of freedom Multiple R-squared: 0.9221, Adjusted R-squared: 0.9156 F-statistic: 142 on 3 and 36 DF, p-value: < 2.2e-16 > summary(lm(y ~ x1)) Call: lm(formula = y ~ x1) Residuals: Min 1Q Median 3Q Max -1.07349 -0.20679 -0.03263 0.31938 0.83930 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.654448 0.073792 49.524 < 2e-16 *** x1.L 2.774274 0.147584 18.798 < 2e-16 *** x1.Q 0.007768 0.147584 0.053 0.958 x1.C -1.258609 0.147584 -8.528 3.65e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4667 on 36 degrees of freedom Multiple R-squared: 0.9221, Adjusted R-squared: 0.9156 F-statistic: 142 on 3 and 36 DF, p-value: < 2.2e-16
順序尺度の時の説明変数が何かというと、粕谷さんのところに解説があるが、対比行列になるとのこと。4水準の場合は以下のようなもの。
> contr.poly(4) .L .Q .C [1,] -0.6708204 0.5 -0.2236068 [2,] -0.2236068 -0.5 0.6708204 [3,] 0.2236068 -0.5 -0.6708204 [4,] 0.6708204 0.5 0.2236068
[追記]
箱ひげ図はわかりにくいというご意見がありましたので、XYプロットも。
plot(y ~ as.numeric(x0), las = 1, xlab = "x", ylab = "y", type = "p", xaxt = "n", col = rgb(0, 0, 0, 0.5), cex = 2, lwd =2) axis(1, 1:4, lv)
タグ:R
箱ヒゲ図ってどうもイメージし難い。
by 春分 (2013-02-02 21:21)
最近だと、半透明重ねうちプロットにした方がわかりやすいかもしれませんね。
by hiroki (2013-02-03 16:29)