SSブログ

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")

Rplot01.png

線形モデルにあてはめてみる。

> 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)

Rplot02.png


タグ:R
nice!(1)  コメント(2)  トラックバック(0) 
共通テーマ:学問

nice! 1

コメント 2

春分

箱ヒゲ図ってどうもイメージし難い。
by 春分 (2013-02-02 21:21) 

hiroki

最近だと、半透明重ねうちプロットにした方がわかりやすいかもしれませんね。
by hiroki (2013-02-03 16:29) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0