My Notes

統計学とかR(R言語)とかPython3の覚え書きとか走り書きとか。 座右の銘にしたい: All work and no play makes Jack a dull boy.

R(R言語)で、重回帰分析における予備的解析(1変数ごとの予備的解析、要約統計量、度数、(不偏)標準偏差、データのグラフ化、ヒストグラム、箱ヒゲ図、ドットプロット、幹葉図)

Rコード

#
# 使用するデータと解説は『SPSSによる回帰分析』p28~33
#


#
# 重回帰分析における予備的解析
# 1変数ごとの予備的解析
#


#
# 使用するデータについて
#
# ある製品の重量y(単位 : mg)と、
# その製品を製造するときの条件x1(熱処理時間 : 秒)、x2(乾燥時間 : 秒)、
# x3(熱処理時間 : 秒)、x4(幅寸法 : mm)を30個の製品を無作為抽出して調べた結果。
#
# このデータを使用して、製造するときの条件x1, x2, x3, x4で製品の重量yを予測する式を作る。
# 回帰式 y = b0 + b1x1 + b2x2 + b3x3 + b4x4
# (説明変数が2つ以上あるので重回帰分析を適用する問題となる)。
#


番号_vec <- c(1:30)

x1_vec <- c(15, 10, 9, 10, 11, 9, 10, 15, 12, 9,
            12, 14, 13, 11, 13, 12, 15, 16, 17, 18,
            16, 15, 15, 18, 19, 20, 8, 11, 18, 12)

x2_vec <- c(105, 99, 75, 103, 102, 87, 95, 111, 110, 105,
            132, 135, 130, 122, 125, 96, 105, 115, 98, 109,
            103, 100, 101, 114, 120, 121, 115, 124, 118, 104)

x3_vec <- c(35, 37, 31, 42, 52, 55, 35, 40, 38, 42,
            46, 48, 47, 45, 49, 41, 50, 44, 45, 46,
            44, 43, 43, 50, 55, 56, 50, 40, 39, 35)

x4_vec <- c(216, 220, 185, 215, 182, 212, 180, 203, 198, 183,
            214, 223, 226, 186, 215, 201, 203, 236, 184, 237,
            224, 220, 197, 195, 224, 221, 195, 231, 186, 197)

y_vec <- c(55, 51, 47, 52, 55, 48, 53, 67, 57, 50,
           59, 75, 62, 56, 73, 60, 88, 74, 80, 80,
           79, 68, 72, 89, 90, 95, 45, 66, 65, 62)

df <- data.frame(番号 = 番号_vec,
                 X1 = x1_vec,
                 X2 = x2_vec,
                 X3 = x3_vec,
                 X4 = x4_vec,
                 Y = y_vec)

df


# 1変数の解析
#
# 要約統計量
summary(df)
# 度数
my_freq <- c(length(df$X1), length(df$X2), length(df$X3), length(df$X4), length(df$Y))
dim(my_freq) <- c(1, 5)
colnames(my_freq) <- c("X1", "X2", "X3", "X4", "Y")
my_freq
# (不偏)標準偏差
my_sd <- c(sd(df$X1), sd(df$X2), sd(df$X3), sd(df$X4), sd(df$Y))
dim(my_sd) <- c(1, 5)
colnames(my_sd) <- c("X1", "X2", "X3", "X4", "Y")
my_sd


# データのグラフ化
par(family = "Osaka") # macOSやOS Xで文字化けするなら。
par(mfrow = c(3, 2))
# ヒストグラム
x1 <- hist(df$X1, right = FALSE)
x1

x2 <- hist(df$X2, right = FALSE)
x2

x3 <- hist(df$X3, right = FALSE)
x3

x4 <- hist(df$X4, right = FALSE)
x4

y <- hist(df$Y, right = FALSE)
y


# 箱ヒゲ図
par(mfrow = c(3, 2))
boxplot(df$X1, xlab = "X1")
boxplot(df$X2, xlab = "X2")
boxplot(df$X3, xlab = "X3")
boxplot(df$X4, xlab = "X4")
boxplot(df$Y, xlab = "Y")


# ドットプロット
#
# データ数が少ない場合は、ヒストグラムや箱ヒゲ図より、ドットプロットがいい。
#
par(mfrow = c(3, 2))
dotchart(df$X1, xlab = "X1")
dotchart(df$X2, xlab = "X2")
dotchart(df$X3, xlab = "X3")
dotchart(df$X4, xlab = "X4")
dotchart(df$Y, xlab = "Y")
# (ドットプロットはggplot2のがいいかな)。


# 幹葉図
#
# データ数が多い場合には、不向き。
#
# (df$X1だけよくわからない数値が出る。試しに、stem(df$番号)でやってみたところ、
# ちゃんとした数値が出るんだが...)。
stem(df$X1)
stem(df$X2)
stem(df$X3)
stem(df$X4)
stem(df$Y)

R Console

> #
> # 使用するデータと解説は『SPSSによる回帰分析』p28~33
> #
> 
> 
> #
> # 重回帰分析における予備的解析
> # 1変数ごとの予備的解析
> #
> 
> 
> #
> # 使用するデータについて
> #
> # ある製品の重量y(単位 : mg)と、
> # その製品を製造するときの条件x1(熱処理時間 : 秒)、x2(乾燥時間 : 秒)、
> # x3(熱処理時間 : 秒)、x4(幅寸法 : mm)を30個の製品を無作為抽出して調べた結果。
> #
> # このデータを使用して、製造するときの条件x1, x2, x3, x4で製品の重量yを予測する式を作る。
> # 回帰式 y = b0 + b1x1 + b2x2 + b3x3 + b4x4
> # (説明変数が2つ以上あるので重回帰分析を適用する問題となる)。
> #
> 
> 
> 番号_vec <- c(1:30)
> 
> x1_vec <- c(15, 10, 9, 10, 11, 9, 10, 15, 12, 9,
+             12, 14, 13, 11, 13, 12, 15, 16, 17, 18,
+             16, 15, 15, 18, 19, 20, 8, 11, 18, 12)
> 
> x2_vec <- c(105, 99, 75, 103, 102, 87, 95, 111, 110, 105,
+             132, 135, 130, 122, 125, 96, 105, 115, 98, 109,
+             103, 100, 101, 114, 120, 121, 115, 124, 118, 104)
> 
> x3_vec <- c(35, 37, 31, 42, 52, 55, 35, 40, 38, 42,
+             46, 48, 47, 45, 49, 41, 50, 44, 45, 46,
+             44, 43, 43, 50, 55, 56, 50, 40, 39, 35)
> 
> x4_vec <- c(216, 220, 185, 215, 182, 212, 180, 203, 198, 183,
+             214, 223, 226, 186, 215, 201, 203, 236, 184, 237,
+             224, 220, 197, 195, 224, 221, 195, 231, 186, 197)
> 
> y_vec <- c(55, 51, 47, 52, 55, 48, 53, 67, 57, 50,
+            59, 75, 62, 56, 73, 60, 88, 74, 80, 80,
+            79, 68, 72, 89, 90, 95, 45, 66, 65, 62)
> 
> df <- data.frame(番号 = 番号_vec,
+                  X1 = x1_vec,
+                  X2 = x2_vec,
+                  X3 = x3_vec,
+                  X4 = x4_vec,
+                  Y = y_vec)
> 
> df
   番号 X1  X2 X3  X4  Y
1     1 15 105 35 216 55
2     2 10  99 37 220 51
3     3  9  75 31 185 47
4     4 10 103 42 215 52
5     5 11 102 52 182 55
6     6  9  87 55 212 48
7     7 10  95 35 180 53
8     8 15 111 40 203 67
9     9 12 110 38 198 57
10   10  9 105 42 183 50
11   11 12 132 46 214 59
12   12 14 135 48 223 75
13   13 13 130 47 226 62
14   14 11 122 45 186 56
15   15 13 125 49 215 73
16   16 12  96 41 201 60
17   17 15 105 50 203 88
18   18 16 115 44 236 74
19   19 17  98 45 184 80
20   20 18 109 46 237 80
21   21 16 103 44 224 79
22   22 15 100 43 220 68
23   23 15 101 43 197 72
24   24 18 114 50 195 89
25   25 19 120 55 224 90
26   26 20 121 56 221 95
27   27  8 115 50 195 45
28   28 11 124 40 231 66
29   29 18 118 39 186 65
30   30 12 104 35 197 62
> 
> 
> # 1変数の解析
> #
> # 要約統計量
> summary(df)
      番号             X1              X2              X3              X4              Y        
 Min.   : 1.00   Min.   : 8.00   Min.   : 75.0   Min.   :31.00   Min.   :180.0   Min.   :45.00  
 1st Qu.: 8.25   1st Qu.:11.00   1st Qu.:101.2   1st Qu.:40.00   1st Qu.:195.0   1st Qu.:55.00  
 Median :15.50   Median :13.00   Median :107.0   Median :44.00   Median :207.5   Median :63.50  
 Mean   :15.50   Mean   :13.43   Mean   :109.3   Mean   :44.10   Mean   :207.0   Mean   :65.77  
 3rd Qu.:22.75   3rd Qu.:15.75   3rd Qu.:119.5   3rd Qu.:48.75   3rd Qu.:220.8   3rd Qu.:74.75  
 Max.   :30.00   Max.   :20.00   Max.   :135.0   Max.   :56.00   Max.   :237.0   Max.   :95.00  
> # 度数
> my_freq <- c(length(df$X1), length(df$X2), length(df$X3), length(df$X4), length(df$Y))
> dim(my_freq) <- c(1, 5)
> colnames(my_freq) <- c("X1", "X2", "X3", "X4", "Y")
> my_freq
     X1 X2 X3 X4  Y
[1,] 30 30 30 30 30
> # (不偏)標準偏差
> my_sd <- c(sd(df$X1), sd(df$X2), sd(df$X3), sd(df$X4), sd(df$Y))
> dim(my_sd) <- c(1, 5)
> colnames(my_sd) <- c("X1", "X2", "X3", "X4", "Y")
> my_sd
           X1       X2     X3       X4        Y
[1,] 3.349558 13.52686 6.4077 17.41082 14.06557
> 
> 
> # データのグラフ化
> par(family = "Osaka") # macOSやOS Xで文字化けするなら。
> par(mfrow = c(3, 2))
> # ヒストグラム
> x1 <- hist(df$X1, right = FALSE)
> x1
$breaks
[1]  8 10 12 14 16 18 20

$counts
[1] 4 6 6 6 3 5

$density
[1] 0.06666667 0.10000000 0.10000000 0.10000000 0.05000000 0.08333333

$mids
[1]  9 11 13 15 17 19

$xname
[1] "df$X1"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
> 
> x2 <- hist(df$X2, right = FALSE)
> x2
$breaks
[1]  70  80  90 100 110 120 130 140

$counts
[1]  1  1  4 10  6  5  3

$density
[1] 0.003333333 0.003333333 0.013333333 0.033333333 0.020000000 0.016666667 0.010000000

$mids
[1]  75  85  95 105 115 125 135

$xname
[1] "df$X2"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
> 
> x3 <- hist(df$X3, right = FALSE)
> x3
$breaks
[1] 30 35 40 45 50 55 60

$counts
[1] 1 6 9 7 4 3

$density
[1] 0.006666667 0.040000000 0.060000000 0.046666667 0.026666667 0.020000000

$mids
[1] 32.5 37.5 42.5 47.5 52.5 57.5

$xname
[1] "df$X3"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
> 
> x4 <- hist(df$X4, right = FALSE)
> x4
$breaks
[1] 180 190 200 210 220 230 240

$counts
[1] 7 5 3 5 7 3

$density
[1] 0.02333333 0.01666667 0.01000000 0.01666667 0.02333333 0.01000000

$mids
[1] 185 195 205 215 225 235

$xname
[1] "df$X4"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
> 
> y <- hist(df$Y, right = FALSE)
> y
$breaks
[1]  40  50  60  70  80  90 100

$counts
[1] 3 9 7 5 4 2

$density
[1] 0.010000000 0.030000000 0.023333333 0.016666667 0.013333333 0.006666667

$mids
[1] 45 55 65 75 85 95

$xname
[1] "df$Y"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
> 
> 
> # 箱ヒゲ図
> par(mfrow = c(3, 2))
> boxplot(df$X1, xlab = "X1")
> boxplot(df$X2, xlab = "X2")
> boxplot(df$X3, xlab = "X3")
> boxplot(df$X4, xlab = "X4")
> boxplot(df$Y, xlab = "Y")
> 
> 
> # ドットプロット
> #
> # データ数が少ない場合は、ヒストグラムや箱ヒゲ図より、ドットプロットがいい。
> #
> par(mfrow = c(3, 2))
> dotchart(df$X1, xlab = "X1")
> dotchart(df$X2, xlab = "X2")
> dotchart(df$X3, xlab = "X3")
> dotchart(df$X4, xlab = "X4")
> dotchart(df$Y, xlab = "Y")
> # (ドットプロットはggplot2のがいいかな)。
> 
> 
> # 幹葉図
> #
> # データ数が多い場合には、不向き。
> #
> # (df$X1だけよくわからない数値が出る。試しに、stem(df$番号)でやってみたところ、
> # ちゃんとした数値が出るんだが...)。
> stem(df$X1)

  The decimal point is at the |

   8 | 0000
  10 | 000000
  12 | 000000
  14 | 000000
  16 | 000
  18 | 0000
  20 | 0

> stem(df$X2)

  The decimal point is 1 digit(s) to the right of the |

   7 | 5
   8 | 7
   9 | 5689
  10 | 0123345559
  11 | 014558
  12 | 01245
  13 | 025

> stem(df$X3)

  The decimal point is 1 digit(s) to the right of the |

  3 | 1
  3 | 555789
  4 | 001223344
  4 | 5566789
  5 | 0002
  5 | 556

> stem(df$X4)

  The decimal point is 1 digit(s) to the right of the |

  18 | 0234566
  19 | 55778
  20 | 133
  21 | 24556
  22 | 0013446
  23 | 167

> stem(df$Y)

  The decimal point is 1 digit(s) to the right of the |

  4 | 578
  5 | 012355679
  6 | 0225678
  7 | 23459
  8 | 0089
  9 | 05

ヒストグラム、箱ヒゲ図、ドットプロットのスクリーンショット

f:id:my_notes:20170819165455p:plain

f:id:my_notes:20170819165510p:plain

f:id:my_notes:20170819165524p:plain

参考文献

SPSSによる回帰分析

SPSSによる回帰分析