第4章のRコード

第4章回帰分析の基礎

サンプルデータ

temperature_aug.csv：2014年8月の東京都の気温データ．

wage.csv：1976年米国における男性労働者の賃金に関するデータ．

パッケージの読み込み

library(readr)
library(ggplot2)
library(dplyr)
library(lubridate)

4.1 回帰分析の考え方

4.1.4 数値例：気温と電力使用量

tempdata <- read_csv("temperature_aug.csv")

ggplot(tempdata, aes(x = time, y = temp)) +
  geom_point() +
  xlab("時刻") +
  ylab("気温") +
  theme_gray(base_family = "HiraKakuPro-W3")

4.1.5 データフレームとティブル

tempdata

# A tibble: 744 × 5
   date      time  elec  prec  temp
   <chr>    <dbl> <dbl> <dbl> <dbl>
 1 2014/8/1     0  3193     0  27.9
 2 2014/8/1     1  2960     0  27.9
 3 2014/8/1     2  2807     0  27.1
 4 2014/8/1     3  2748     0  26.8
 5 2014/8/1     4  2735     0  26.9
 6 2014/8/1     5  2736     0  27.3
 7 2014/8/1     6  2950     0  28.3
 8 2014/8/1     7  3336     0  29.4
 9 2014/8/1     8  3863     0  30.2
10 2014/8/1     9  4328     0  32.2
# ℹ 734 more rows

4.1.6 ノンパラメトリック回帰の実行

tempdata %>% 
  summarise(mean_temp = mean(temp))

# A tibble: 1 × 1
  mean_temp
      <dbl>
1      27.7

tempdata %>% 
  summarise(mean_temp = mean(temp),
            sd_temp   = sd(temp),
            max_temp  = max(temp),
            min_temp  = min(temp))

# A tibble: 1 × 4
  mean_temp sd_temp max_temp min_temp
      <dbl>   <dbl>    <dbl>    <dbl>
1      27.7    3.55     35.5     19.8

tempdata %>%
  group_by(time) %>% 
  summarise(mean_temp = mean(temp))

# A tibble: 24 × 2
    time mean_temp
   <dbl>     <dbl>
 1     0      26.4
 2     1      26.1
 3     2      25.9
 4     3      25.8
 5     4      25.7
 6     5      25.8
 7     6      26.4
 8     7      27.3
 9     8      28.1
10     9      29.0
# ℹ 14 more rows

tempdata %>%
  group_by(time) %>%
  summarise(mean_temp = mean(temp)) %>%
  ggplot(aes(x = time, y = mean_temp)) +
  geom_line() +
  xlab("時間") +
  ylab("気温") + 
  theme_gray(base_family = "HiraKakuPro-W3")

tempdata %>% 
  ggplot(aes(x = time, y = temp)) +
  stat_summary(geom = "line", fun = "mean") +
  xlab("時間") + 
  ylab("気温") +
  theme_gray(base_family = "HiraKakuPro-W3")

4.1.7 グラフを重ねる

tempdata %>% 
  ggplot(aes(x = time, y = temp)) +
  geom_point() +
  stat_summary(geom = "line", fun = "mean") +
  xlab("時間") + 
  ylab("気温") +
  coord_cartesian(ylim = c(20, 30)) +
  theme_gray(base_family = "HiraKakuPro-W3")

4.2 単回帰分析

4.2.1 ノンパラメトリック回帰の限界

with(tempdata, cor(temp, elec))

[1] 0.7198181

tempdata %>% 
  summarise(cor(temp, elec))

# A tibble: 1 × 1
  `cor(temp, elec)`
              <dbl>
1             0.720

tempdata %>% 
  with(cor(temp, elec))

[1] 0.7198181

図4.6

tempdata %>% 
  ggplot(aes(x = temp, y = elec)) +
  stat_summary(geom = "line", fun = "mean") + 
  xlab("電気使用量の回帰曲線") +
  ylab("電気使用量（万kw）") +
  theme_gray() +
  theme_gray(base_family = "HiraKakuPro-W3")

4.2.4 Rによる線形回帰分析

lm(elec ~ temp,
   data = tempdata)


Call:
lm(formula = elec ~ temp, data = tempdata)

Coefficients:
(Intercept)         temp  
     -614.3        145.0

4.2.6 単回帰の図示

result <- lm(elec ~ temp,
             data = tempdata)

result$coefficients

(Intercept)        temp 
  -614.2722    145.0303

tempdata %>% 
  ggplot(aes(x = temp, y = elec)) +
  geom_point() +
  geom_abline(intercept = result$coefficients[1],
              slope     = result$coefficients[2]) +
  xlab("気温") +
  ylab("電気使用量（万kw）") +
  theme_gray(base_family = "HiraKakuPro-W3")

tempdata %>% 
  ggplot(aes(x = temp, y = elec)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  xlab("気温") +
  ylab("電気使用量（万kw）") +
  theme_gray(base_family = "HiraKakuPro-W3")

4.3 重回帰分析

4.3.1 説明変数を追加する

tempdata %>% 
  ggplot(aes(x = time, y = elec)) +
  geom_point() +
  xlab("時刻") +
  ylab("電気使用量（万kw）") +
  theme_gray(base_family = "HiraKakuPro-W3")

tempdata$daytime <- 
  (tempdata$time >= 9) & (tempdata$time <= 18)

tempdata <- tempdata %>% 
  mutate(daytime = 1 * (9 <= time & time <= 18),
         elec100 = elec / 100,
         time12  = time %% 12,
         ampm    = ifelse(time < 12, "a.m.", "p.m."))

4.3.2 重回帰分析

lm(elec ~ temp + daytime,
   data = tempdata)


Call:
lm(formula = elec ~ temp + daytime, data = tempdata)

Coefficients:
(Intercept)         temp      daytime  
     -69.55       116.98       555.44

lm(elec ~ temp + daytime + prec,
    data = tempdata)


Call:
lm(formula = elec ~ temp + daytime + prec, data = tempdata)

Coefficients:
(Intercept)         temp      daytime         prec  
    -64.503      116.802      556.148       -3.366

tempdata %>%
  mutate(sunday = 1 * 
           (date == "2014/8/3")  |
           (date == "2014/8/10") |
           (date == "2014/8/17") |
           (date == "2014/8/24") |
           (date == "2014/8/31"))

# A tibble: 744 × 10
   date      time  elec  prec  temp daytime elec100 time12 ampm  sunday
   <chr>    <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>  <dbl> <chr> <lgl> 
 1 2014/8/1     0  3193     0  27.9       0    31.9      0 a.m.  FALSE 
 2 2014/8/1     1  2960     0  27.9       0    29.6      1 a.m.  FALSE 
 3 2014/8/1     2  2807     0  27.1       0    28.1      2 a.m.  FALSE 
 4 2014/8/1     3  2748     0  26.8       0    27.5      3 a.m.  FALSE 
 5 2014/8/1     4  2735     0  26.9       0    27.4      4 a.m.  FALSE 
 6 2014/8/1     5  2736     0  27.3       0    27.4      5 a.m.  FALSE 
 7 2014/8/1     6  2950     0  28.3       0    29.5      6 a.m.  FALSE 
 8 2014/8/1     7  3336     0  29.4       0    33.4      7 a.m.  FALSE 
 9 2014/8/1     8  3863     0  30.2       0    38.6      8 a.m.  FALSE 
10 2014/8/1     9  4328     0  32.2       1    43.3      9 a.m.  FALSE 
# ℹ 734 more rows

tempdata <- tempdata %>% 
  mutate(date = ymd(date))
tempdata

# A tibble: 744 × 9
   date        time  elec  prec  temp daytime elec100 time12 ampm 
   <date>     <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>  <dbl> <chr>
 1 2014-08-01     0  3193     0  27.9       0    31.9      0 a.m. 
 2 2014-08-01     1  2960     0  27.9       0    29.6      1 a.m. 
 3 2014-08-01     2  2807     0  27.1       0    28.1      2 a.m. 
 4 2014-08-01     3  2748     0  26.8       0    27.5      3 a.m. 
 5 2014-08-01     4  2735     0  26.9       0    27.4      4 a.m. 
 6 2014-08-01     5  2736     0  27.3       0    27.4      5 a.m. 
 7 2014-08-01     6  2950     0  28.3       0    29.5      6 a.m. 
 8 2014-08-01     7  3336     0  29.4       0    33.4      7 a.m. 
 9 2014-08-01     8  3863     0  30.2       0    38.6      8 a.m. 
10 2014-08-01     9  4328     0  32.2       1    43.3      9 a.m. 
# ℹ 734 more rows

wday(tempdata$date, label = TRUE, locale = "ja_JP")

  [1] 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 土
 [26] 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 日 日
 [51] 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 月 月 月
 [76] 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 火 火 火 火
[101] 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 水 水 水 水 水
[126] 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 木 木 木 木 木 木
[151] 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 金 金 金 金 金 金 金
[176] 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 土 土 土 土 土 土 土 土
[201] 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 日 日 日 日 日 日 日 日 日
[226] 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 月 月 月 月 月 月 月 月 月 月
[251] 月 月 月 月 月 月 月 月 月 月 月 月 月 月 火 火 火 火 火 火 火 火 火 火 火
[276] 火 火 火 火 火 火 火 火 火 火 火 火 火 水 水 水 水 水 水 水 水 水 水 水 水
[301] 水 水 水 水 水 水 水 水 水 水 水 水 木 木 木 木 木 木 木 木 木 木 木 木 木
[326] 木 木 木 木 木 木 木 木 木 木 木 金 金 金 金 金 金 金 金 金 金 金 金 金 金
[351] 金 金 金 金 金 金 金 金 金 金 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土
[376] 土 土 土 土 土 土 土 土 土 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日
[401] 日 日 日 日 日 日 日 日 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月
[426] 月 月 月 月 月 月 月 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火
[451] 火 火 火 火 火 火 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水
[476] 水 水 水 水 水 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木
[501] 木 木 木 木 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金
[526] 金 金 金 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土
[551] 土 土 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日
[576] 日 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月 月
[601] 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 火 水
[626] 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 水 木 木
[651] 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 木 金 金 金
[676] 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 金 土 土 土 土
[701] 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 土 日 日 日 日 日
[726] 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日 日
Levels: 日 < 月 < 火 < 水 < 木 < 金 < 土

tempdata <- tempdata %>% 
  mutate(dow    = wday(date, label = TRUE, locale = "ja_JP"),
         sunday = 1 * (dow == "日"))

tempdata <- tempdata %>% 
  mutate(recess = 1 * ("2014-08-11" <= date & 
                         date <= "2014-08-16"))

lm(elec ~ temp + daytime + prec + sunday + recess,
   data = tempdata)


Call:
lm(formula = elec ~ temp + daytime + prec + sunday + recess, 
    data = tempdata)

Coefficients:
(Intercept)         temp      daytime         prec       sunday       recess  
     179.07       113.48       563.53        14.27      -448.39      -438.23

179.07 + 113.48 * 28 + 563.53 * 1 + 14.27 * 0 - 448.39 * 0 - 438.23 * 0

[1] 3920.04

result <- lm(elec ~ temp + daytime + prec + sunday + recess,
             data = tempdata)
sum(result$coefficients * c(1, 28, 1, 0, 0, 0))

[1] 3919.964

4.4 決定係数と回帰分析

4.4.2 決定係数の出力

lm(elec ~ temp + daytime + prec + sunday + recess,
   data = tempdata) %>% 
  summary()


Call:
lm(formula = elec ~ temp + daytime + prec + sunday + recess, 
    data = tempdata)

Residuals:
    Min      1Q  Median      3Q     Max 
-825.22 -289.33    1.57  270.66  999.57 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  179.071    114.088   1.570    0.117    
temp         113.477      4.193  27.066   <2e-16 ***
daytime      563.528     29.716  18.964   <2e-16 ***
prec          14.275     14.701   0.971    0.332    
sunday      -448.392     38.125 -11.761   <2e-16 ***
recess      -438.230     35.199 -12.450   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 370.3 on 738 degrees of freedom
Multiple R-squared:  0.7331,    Adjusted R-squared:  0.7313 
F-statistic: 405.5 on 5 and 738 DF,  p-value: < 2.2e-16

lm(elec ~ temp,
   data = tempdata) %>% 
  summary()


Call:
lm(formula = elec ~ temp, data = tempdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-1105.11  -421.64    11.39   394.19  1065.64 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -614.272    143.223  -4.289 2.03e-05 ***
temp         145.030      5.134  28.246  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 496.2 on 742 degrees of freedom
Multiple R-squared:  0.5181,    Adjusted R-squared:  0.5175 
F-statistic: 797.9 on 1 and 742 DF,  p-value: < 2.2e-16

4.4.3 決定係数は「モデルの正しさ」を保証しない

set.seed(2022)

a <- sample(1:6, 15, replace = TRUE)
b <- sample(1:6, 15, replace = TRUE)
c <- sample(1:6, 15, replace = TRUE)
d <- sample(1:6, 15, replace = TRUE)
e <- sample(1:6, 15, replace = TRUE)
f <- sample(1:6, 15, replace = TRUE)
g <- sample(1:6, 15, replace = TRUE)

summary(lm(a ~ b + c + d + e + f + g))


Call:
lm(formula = a ~ b + c + d + e + f + g)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6875 -0.8476  0.2100  0.6657  1.9867 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  5.97169    1.97662   3.021   0.0165 *
b            0.01211    0.34577   0.035   0.9729  
c           -0.27789    0.31640  -0.878   0.4054  
d            0.37124    0.31947   1.162   0.2787  
e           -0.30285    0.22822  -1.327   0.2211  
f           -0.52931    0.27004  -1.960   0.0856 .
g            0.30288    0.26822   1.129   0.2915  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.328 on 8 degrees of freedom
Multiple R-squared:  0.618, Adjusted R-squared:  0.3316 
F-statistic: 2.157 on 6 and 8 DF,  p-value: 0.1552

h <- sample(1:6, 15, replace = TRUE)

lm(a ~ b + c + d + e + f + g + h) %>% summary()


Call:
lm(formula = a ~ b + c + d + e + f + g + h)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.24977 -0.23138 -0.09301  0.53529  0.78575 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  2.82566    1.69156   1.670   0.1388  
b            0.38667    0.26654   1.451   0.1901  
c            0.01289    0.23694   0.054   0.9581  
d            0.05239    0.24250   0.216   0.8351  
e           -0.21587    0.15971  -1.352   0.2185  
f           -0.54584    0.18619  -2.932   0.0220 *
g           -0.11759    0.22833  -0.515   0.6224  
h            0.63108    0.20116   3.137   0.0164 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9152 on 7 degrees of freedom
Multiple R-squared:  0.8413,    Adjusted R-squared:  0.6825 
F-statistic: 5.299 on 7 and 7 DF,  p-value: 0.0214

補足

tempdata <- readr::read_csv("temperature_aug.csv")

result <- lm(elec ~ temp,
             data = tempdata)

ehat <- result$residuals

sum(ehat)

[1] 2.615153e-11

sum(tempdata$temp * ehat)

[1] 1.093596e-09

第4章 回帰分析の基礎

パッケージの読み込み

4.1 回帰分析の考え方

4.1.4 数値例：気温と電力使用量

4.1.5 データフレームとティブル

4.1.6 ノンパラメトリック回帰の実行

4.1.7 グラフを重ねる

4.2 単回帰分析

4.2.1 ノンパラメトリック回帰の限界

図4.6

4.2.4 Rによる線形回帰分析

4.2.6 単回帰の図示

4.3 重回帰分析

4.3.1 説明変数を追加する

4.3.2 重回帰分析

4.4 決定係数と回帰分析

4.4.2 決定係数の出力

4.4.3 決定係数は「モデルの正しさ」を保証しない

補足

第4章回帰分析の基礎