My Notes

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

アイスクリーム統計学 第5章 偏相関をR(R言語)で 試しに偏相関係数を計算する関数を作成

Rコード

# アイスクリーム統計学 ch05




#
# 偏相関
#




# 最低気温のデータも含めたすべてのデータ
データ番号_vec <- c(1:20)

最高気温_vec <- c(33, 33, 34, 34, 35, 35, 34, 32, 28, 35,
                33, 28, 32, 25, 28, 30, 29, 32, 34, 35)

最低気温_vec <- c(22, 26, 27, 28, 28, 27, 28, 25, 24, 24,
                26, 25, 23, 22, 21, 23, 23, 25, 26, 27)

客数_vec <- c(382, 324, 338, 317, 341, 360, 339, 329, 218, 402,
            342, 205, 368, 196, 304, 294, 275, 336, 384, 368)

df <- data.frame(データ番号 = データ番号_vec,
                最高気温 = 最高気温_vec,
                最低気温 = 最低気温_vec,
                客数 = 客数_vec)

df




summary(df)




# 最低気温と客数の相関
#
# 散布図
par(family = "Osaka")
plot(df$最低気温, df$客数,
    pch = 18, col = "blue", cex = 3,
    xlim = c(20, 30), ylim = c(0, 500),
    xlab = "最低気温", ylab = "客数", main = "最低気温と客数の関係")




# 最低気温と客数の相関係数
cor(df$最低気温, df$客数)
round(cor(df$最低気温, df$客数) ,3)
cor.test(df$最低気温, df$客数) # (デフォルトのピアソンでいいのか不安はある。
                             # ピアソンを使用できる前提となる要件に合致してるのかどうかなど)
# p-value = 0.1958なので、相関の有意性の検定の帰無仮説である、"変数は無関係である。r = 0"は、
# 棄却できない。




# 偏相関
#
# 各相関係数
最高気温と客数の相関係数 <- round(cor(df$最高気温, df$客数), 3)
最高気温と客数の相関係数

最低気温と客数の相関係数 <- round(cor(df$最低気温, df$客数) ,3)
最低気温と客数の相関係数

最高気温と最低気温の相関係数 <- round(cor(df$最高気温, df$最低気温), 3)
最高気温と最低気温の相関係数

# aの影響を除いたbとyの偏相関係数
# aは最高気温、bは最低気温、yは客数
# 数式 (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
rab <- 最高気温と最低気温の相関係数
rab

ray <- 最高気温と客数の相関係数
ray

rby <- 最低気温と客数の相関係数
rby

#偏相関係数の手計算 最高気温の影響を取り除いた偏相関係数
偏相関係数1 <- (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
偏相関係数1
round(偏相関係数1, 3)

# 最低気温の影響を取り除いた偏相関係数
偏相関係数2 <- (ray - (rby * rab)) / (sqrt(1 - rby^2) * sqrt(1 - rab^2))
偏相関係数2
round(偏相関係数2, 3)


# 試しに関数化
偏相関係数_func1 <- function(ray, rby, rab) {
    (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
}
偏相関係数_func1(最高気温と客数の相関係数, 最低気温と客数の相関係数, 最高気温と最低気温の相関係数)
# もっと一般化できるが、とりあえず。
# 上記、関数では具体的すぎて、汎用性がない。
# まともな関数にするなら、cor()の計算も関数内で行えばいい。
# 引数は、例えば、最高気温、最低気温、客数にする。
# 引数をベクトルにすれば、もう少しまともな関数になるだろう。
#
# 例えば、
偏相関係数_func2 <- function(a, b, y) {
    rab <- round(cor(a, b), 3)

    ray <- round(cor(a, y), 3)

    rby <- round(cor(b, y), 3)

    (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
}
偏相関係数_func2(df$最高気温, df$最低気温, df$客数) # 第1引数に影響を除外したいデータ
偏相関係数_func2(df$最低気温, df$最高気温, df$客数)
# だが、これでも一般化は足りない。
# 例外処理もない。
# いまのところ使用頻度が高いとは思えないので、これ以上深入りはしない。

R Console

> # アイスクリーム統計学 ch05
> 
> 
> 
> 
> #
> # 偏相関
> #
> 
> 
> 
> 
> # 最低気温のデータも含めたすべてのデータ
> データ番号_vec <- c(1:20)
> 
> 最高気温_vec <- c(33, 33, 34, 34, 35, 35, 34, 32, 28, 35,
+                 33, 28, 32, 25, 28, 30, 29, 32, 34, 35)
> 
> 最低気温_vec <- c(22, 26, 27, 28, 28, 27, 28, 25, 24, 24,
+                 26, 25, 23, 22, 21, 23, 23, 25, 26, 27)
> 
> 客数_vec <- c(382, 324, 338, 317, 341, 360, 339, 329, 218, 402,
+             342, 205, 368, 196, 304, 294, 275, 336, 384, 368)
> 
> df <- data.frame(データ番号 = データ番号_vec,
+                 最高気温 = 最高気温_vec,
+                 最低気温 = 最低気温_vec,
+                 客数 = 客数_vec)
> 
> df
   データ番号 最高気温 最低気温 客数
1           1       33       22  382
2           2       33       26  324
3           3       34       27  338
4           4       34       28  317
5           5       35       28  341
6           6       35       27  360
7           7       34       28  339
8           8       32       25  329
9           9       28       24  218
10         10       35       24  402
11         11       33       26  342
12         12       28       25  205
13         13       32       23  368
14         14       25       22  196
15         15       28       21  304
16         16       30       23  294
17         17       29       23  275
18         18       32       25  336
19         19       34       26  384
20         20       35       27  368
> 
> 
> 
> 
> summary(df)
   データ番号       最高気温        最低気温       客数      
 Min.   : 1.00   Min.   :25.00   Min.   :21   Min.   :196.0  
 1st Qu.: 5.75   1st Qu.:29.75   1st Qu.:23   1st Qu.:301.5  
 Median :10.50   Median :33.00   Median :25   Median :337.0  
 Mean   :10.50   Mean   :31.95   Mean   :25   Mean   :321.1  
 3rd Qu.:15.25   3rd Qu.:34.00   3rd Qu.:27   3rd Qu.:362.0  
 Max.   :20.00   Max.   :35.00   Max.   :28   Max.   :402.0  
> 
> 
> 
> 
> # 最低気温と客数の相関
> #
> # 散布図
> par(family = "Osaka")
> plot(df$最低気温, df$客数,
+     pch = 18, col = "blue", cex = 3,
+     xlim = c(20, 30), ylim = c(0, 500),
+     xlab = "最低気温", ylab = "客数", main = "最低気温と客数の関係")
> 
> 
> 
> 
> # 最低気温と客数の相関係数
> cor(df$最低気温, df$客数)
[1] 0.3019116
> round(cor(df$最低気温, df$客数) ,3)
[1] 0.302
> cor.test(df$最低気温, df$客数) # (デフォルトのピアソンでいいのか不安はある。

    Pearson's product-moment correlation

data:  df$最低気温 and df$客数
t = 1.3436, df = 18, p-value = 0.1958
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1622917  0.6566963
sample estimates:
      cor 
0.3019116 

>                              # ピアソンを使用できる前提となる要件に合致してるのかどうかなど)
> # p-value = 0.1958なので、相関の有意性の検定の帰無仮説である、"変数は無関係である。r = 0"は、
> # 棄却できない。
> 
> 
> 
> 
> # 偏相関
> #
> # 各相関係数
> 最高気温と客数の相関係数 <- round(cor(df$最高気温, df$客数), 3)
> 最高気温と客数の相関係数
[1] 0.87
> 
> 最低気温と客数の相関係数 <- round(cor(df$最低気温, df$客数) ,3)
> 最低気温と客数の相関係数
[1] 0.302
> 
> 最高気温と最低気温の相関係数 <- round(cor(df$最高気温, df$最低気温), 3)
> 最高気温と最低気温の相関係数
[1] 0.706
> 
> # aの影響を除いたbとyの偏相関係数
> # aは最高気温、bは最低気温、yは客数
> # 数式 (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
> rab <- 最高気温と最低気温の相関係数
> rab
[1] 0.706
> 
> ray <- 最高気温と客数の相関係数
> ray
[1] 0.87
> 
> rby <- 最低気温と客数の相関係数
> rby
[1] 0.302
> 
> #偏相関係数の手計算 最高気温の影響を取り除いた偏相関係数
> 偏相関係数1 <- (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
> 偏相関係数1
[1] -0.894139
> round(偏相関係数1, 3)
[1] -0.894
> 
> # 最低気温の影響を取り除いた偏相関係数
> 偏相関係数2 <- (ray - (rby * rab)) / (sqrt(1 - rby^2) * sqrt(1 - rab^2))
> 偏相関係数2
[1] 0.9728118
> round(偏相関係数2, 3)
[1] 0.973
> 
> 
> # 試しに関数化
> 偏相関係数_func1 <- function(ray, rby, rab) {
+     (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
+ }
> 偏相関係数_func1(最高気温と客数の相関係数, 最低気温と客数の相関係数, 最高気温と最低気温の相関係数)
[1] -0.894139
> # もっと一般化できるが、とりあえず。
> # 上記、関数では具体的すぎて、汎用性がない。
> # まともな関数にするなら、cor()の計算も関数内で行えばいい。
> # 引数は、例えば、最高気温、最低気温、客数にする。
> # 引数をベクトルにすれば、もう少しまともな関数になるだろう。
> #
> # 例えば、
> 偏相関係数_func2 <- function(a, b, y) {
+     rab <- round(cor(a, b), 3)
+ 
+     ray <- round(cor(a, y), 3)
+ 
+     rby <- round(cor(b, y), 3)
+ 
+     (rby - (ray * rab)) / (sqrt(1 - ray^2) * sqrt(1 - rab^2))
+ }
> 偏相関係数_func2(df$最高気温, df$最低気温, df$客数) # 第1引数に影響を除外したいデータ
[1] -0.894139
> 偏相関係数_func2(df$最低気温, df$最高気温, df$客数)
[1] 0.9728118
> # だが、これでも一般化は足りない。
> # 例外処理もない。
> # いまのところ使用頻度が高いとは思えないので、これ以上深入りはしない。

散布図のスクリーンショット

f:id:my_notes:20170702035803p:plain

参考Webサイト

アイスクリーム統計学にようこそ!

書籍はこち

統計学がわかる 【回帰分析・因子分析編】 (ファーストブック)

統計学がわかる 【回帰分析・因子分析編】 (ファーストブック)