重回帰分析

最小二乗法による重回帰パラメータの算出

  • read.csv関数を使って,以下のURLのデータを読み込み,オブジェクト「S16」に代入(付値)するプログラムを書け.
    • https://fnakai.web.nitech.ac.jp/data/14_LP/16Students.csv
  • sum関数を用いて\(\sum x, \sum y, \sum z, \sum x^2, \sum z^2,\sum xy,\sum xz,\sum zy\)を求めよ.
  • 重回帰パラメータを導出する正規方程式に,上記の値を当てはめた連立方程式を書け(解かなくて良い).
  • lm関数を使ってパラメータ\(a,b,c\)を求め,回帰式を書け.

sum関数について

以下を実行してみてほしい.sum関数は,sum( )の括弧の中にベクトルを入れると,その成分の和を返すことがわかる.

a <- c(1,1,1,1,1,1)
sum(a)

以下を実行してみてほしい.「*」はかけ算の演算子である.また,S16$Weight[i]は,Weightのi番目の成分を意味する.ベクトル同士に対して「*」によるかけ算をおこなうことで,どのような演算結果が返されるか?

S16$Weight * S16$Height
S16$Weight[1] * S16$Height[1]

lm関数について

lm関数は,回帰モデルを求める関数である.以下のformulaでは回帰モデルの形を指定する.formulaにはいくつかの種類があるので,以下のWebサイトを参考にしてほしい.

参考ウェブサイト:回帰分析と重回帰分析

formulaを「y ~ x」と指定することで,モデル式 \(y = a + bx + \epsilon\)( \(\epsilon\)は誤差項)の回帰モデルを求めることができる.被説明変数 y と説明変数 x は,ベクトルで指定する.

lm(formula)

たとえば,以下は体重を被説明変数,身長を説明変数とした場合である.

lm(S16$Weight~S16$Height)

formulaを「y ~ x+z」と指定することで,モデル式 \(y = a + bx + cz + \epsilon\)の回帰モデルを求めることができる.

◆レポート提出課題(全10問)

1.read.csv関数を使って,以下のURLのデータを読み込み,オブジェクト「S16」に代入(付値)するプログラムを示しなさい.

2.「*」はかけ算の演算子である.ベクトル同士に対して「*」によるかけ算をおこなうことで,どのような演算結果が返されるか?以下を実行することにより,考察せよ.

S16$Weight * S16$Height
S16$Weight[1] * S16$Height[1]
S16$Weight[2] * S16$Height[2]

3.sum関数を用いて\(\sum x, \sum y, \sum z, \sum x^2, \sum z^2,\sum xy,\sum xz,\sum zy\)を求めるプログラムを示しなさい.

4.重回帰パラメータを導出する正規方程式に,上記の値を当てはめた連立方程式を書け(解かなくて良い).

5.lm関数を用いて,被説明変数に「身長」を,説明変数に「体重」と「胸囲」を置いた重回帰モデルを作成せよ.また,パラメータa,b,cを明示したうえで,重回帰式を示せ.

6.lm関数の結果をオブジェクト「result」に代入(付値)するプログラムを示しなさい.

7.summary(result)を実行せよ.また,以下の表を参考に,結果について以下の点を考察せよ.

  • モデル式とデータの入力情報について日本語で説明せよ.
  • 残差とは,回帰式と実データによって求められる偏差\(y_i – \overline{y}_{xz}\)のことである.その要約統計量を示しているのが,Residualsである.これより,残差についてわかることや,気が付いたことを説明せよ.
  • 重決定係数の値を求めよ.また,sqrt関数を用いることにより,重相関係数を求めよ(sqrt関数については「?sqrt」と入力して,ヘルプを見ること).
R上の表記説明
Call:モデル式とデータの入力情報.
Residuals:残差の四分位数に関する情報.
Coefficients:Estimate: 回帰係数の推計値\(\overline{y}_{xz}\)
Std. Error:回帰係数の標準誤差
t value:回帰係数t値.推定値と標準誤差の比.
Pr(>|t|):回帰係数のp値.
Residual standard error:残差の標準誤差.誤差項の標準偏差の推定値.
Multiple R-squared:決定係数(※通常の/演習でもとめている決定係数のこと)
Adjusted R-squared:自由度調整済み決定係数
F-statistic:F 値.回帰式が意味があるかどうかを検定する統計量。

p-value:
(F 検定に基づく)p 値.回帰式が意味が無い(全ての説明変数の係数がゼロである)確率
summary(result)の結果の読み取り方

8.lm関数を用いて,被説明変数に「身長」を,説明変数に「体重」のみを置いた単回帰モデルを作成せよ.この結果は「result2」に代入せよ.また,パラメータa,bを明示したうえで,単回帰式を示せ.

9.8の単回帰モデルについて,summary(result2)を実行することにより,重回帰分析の時と同様に結果を考察せよ.

10.単回帰モデルの結果と重回帰モデルのsummary関数の実行結果をもとに結果の比較を行い,どちらのほうが「身長」をよく説明するモデルか,考察せよ.

◆レポートの作成方法と提出方法

作成方法

  • wordファイルで作成すること.
  • 以下の内容を必ず含むこと
    • 講義名,学籍番号,名前,提出日
    • 問の番号と問の内容,解答
    • 図を貼りつける場合,図の番号と図のキャプション

提出方法

  • moodle上で以下2点を提出すること.
    • wordファイル
    • 課題を解く際に作成したRのプログラム

自転車駐輪場の利用利用金収入額と駅の平均乗降客数の分析

教科書p.60[事例5-b]の分析

回帰式を求め,データと共に描画します.

これを描くRのコードは以下です.

# coded by F.Nakai
# 2021.01.20

#データの読み込み
SIU <- read.csv("https://fnakai.web.nitech.ac.jp/data/13_LP/StationIncomeAndUsers.csv")

#Rのlm関数による回帰パラメータの導出
result <- lm(SIU$Income~SIU$UsersNum)
result

##問2 手計算の過程の再現
#yの偏差
yy <- (SIU$Income-mean(SIU$Income))
#xの偏差
xx <- (SIU$UsersNum-mean(SIU$UsersNum))

#xの分散とxとyの共分散を求める
Sxy <- sum(xx*yy)*(1/length(xx))
Sxx <- sum(xx^2)*(1/length(xx))

# b=(xとyの共分散/xの分散)
b <- Sxy/Sxx

# a= (yの平均)-b*(xの平均) 
a <- (mean(SIU$Income))-(b*mean(SIU$UsersNum))

#教科書p.61の図を描く
plot(SIU$UsersNum,SIU$Income,
     xlim = c(0,300),
     ylim = c(0,250),
     main = "自転車駐輪場の利用利用金収入額と駅の平均乗降客数",
     xlab = "駅の平均乗降客数(百人/日)",
     ylab = "自転車駐輪場の利用料金収入額(万円/月)"
)

#回帰直線を重ねる
abline(a = a, b = b)

##発展編:きれいに図を描くggplotパッケージの利用

#パッケージのインストール
install.packages("ggplot2")
install.packages("ggpmisc")

#パッケージの読み込み
library(ggplot2)
library(ggpmisc)

#回帰のモデル式を保存しておく
my.formula <- y ~ x

#図を描く
ggplot(data = SIU, mapping = aes(x = UsersNum, y = Income)) +
  geom_point() +
  labs(title = "自転車駐輪場の利用利用金収入額と駅の平均乗降客数") +
  xlab("駅の平均乗降客数(百人/日)") +
  ylab("自転車駐輪場の利用料金収入額(万円/月)") +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE, color="black", formula = my.formula, size = 0.5)+
  stat_poly_eq(formula = my.formula, 
               aes(label = paste(..eq.label..)), 
               parse = TRUE) +
  theme_gray()

#問2 Table3の計算
table3 <-cbind(SIU,yy,xx,yy^2,xx^2,xx*yy)
table3

##問3 決定係数を求める
#yの予測値からの偏差平方和の平均(配布資料の2.2式)
Sxy2 <- (1/length(SIU$Income))*sum((SIU$Income - predict(result))^2)
Sxy2

#yの分散
Syy <- sum(yy^2)*(1/length(yy))
Syy

#決定係数
r2 <- 1-(Sxy2/Syy)

#相関係数
r <- sqrt(r2)

#Rのlm関数によって求めた決定係数との比較:summary関数を用いた表示
summary(result)

#Rのcor関数によって求めた相関係数との比較
cor(SIU$Income,SIU$UsersNum)

多変量データを扱う

      ◆データの読み込み

      多変量解析で用いるデータを,外部(パソコン上や,インターネット上)からRへ読み込む.今回の演習で扱うデータは,CSV(Comma-Separated Values)と呼ばれるデータ形式のものを利用する.CSVは,Excelやテキストエディタなどのソフトウェアで開くことが出来る.

      CSV形式のデータファイルは,read.csv関数で読み込むことが出来る.read.csvのあとの括弧()の中にある「file」や「header = TRUE」のことを引数(ひきすう)という.引数は複数指定することができ,左から「第一引数」「第二引数」という.第一引数は,読み込みたいファイルが置いてあるフォルダのパスとファイル名をつなげて指定する.

      read.csv(file, header = TRUE, sep = ",",...)

      ◆コロナウイルス感染者データの入手と読み込み

      都道府県別新型コロナウイルス感染者数マップ URL では,発表された症例一覧のCSVファイルが入手できる.

      コロナウイルス感染者データの読み込みと表示

      このCSVデータにはいくつかの欠損値などがあり扱いにくい.そこで,今回の演習用に,加工したデータを以下のURLに置いた.

      https://fnakai.web.nitech.ac.jp/data/COVID-19_modified.csv

      このデータをread.csv関数を使って読み込み,読み込んだファイルを「COVID-19」というオブジェクトに代入(付値)せよ.URL指定による方法でも,ローカルフォルダに保存してから読み込む方法でも,どちらでも良い.

      第一引数には,ファイルのパス/ファイル名を” “(ダブルクオーテーション)でくくって指定する.また,第二引数には,header= TRUEを指定する.これは「1行目は列名を意味しますよ」ということを示すものである.第3引数にはsep=”,”を指定する.sepという引数は,読み込むファイルが ,(カンマ)で区切られていることを示すものである.

      1.オンライン上のデータを読み込む方法(URLを指定してアクセス)

      ※今回の演習ではRstudio Cloudを使う都合上,こちらの方法で読み込む

      第一引数には,以下のURLを” “(ダブルクオーテーション)でくくって指定する.

      https://fnakai.web.nitech.ac.jp/data/COVID-19_modified.csv

      2.ローカルフォルダのデータを読み込む方法(フォルダのパスを指定してアクセス)

      データをダウンロードし,EIP2020に保存する.そして,そのデータをread.csv関数を使って読み込む.

      第一引数には,(1)フォルダのフルパス”C:/Users/Fuko Nakai/Desktop/EIP2020/NagoyaCity_population_en.csv “を指定する方法,(2)ファイル名”NagoyaCity_population_en.csv “だけ指定する方法がある.(2)は,CWDを省略した表記である.今回はCWDがEIP2020に指定されているため省略可能であるが,CWDとは異なるパスのファイルを読み込む場合には,フルパスを指定する.

      以下は,一度ローカルフォルダに保存してから読み込むケースである.colClassesという引数を用いて各列のデータ型を指定する.

      COVID19 <- read.csv("COVID-19_modified.csv",
                          sep=",",
                          colClasses = c("integer","integer","character","character","character","Date"))

      オブジェクト名を入力し,実行すると,読み込んだデータの中身を見ることができる.「COVID19」と入力し,中身を確認せよ.また,以下の表に示す列が現れることを確認せよ.

      COVID19
      データの列名内容
      ID本CSV内におけるデータの通し番号
      Age公表された年代。10歳刻み。 10歳未満の場合「0-10」と記載。実際の年齢が公表された場合にも年代に丸めて記載。不明もしくは非公表の場合には「不明」と記載。
      Hospital Pref受診都道府県
      Residential Pref居住都道府県(居住地)
      Gender公表された性別、男性、女性。 不明もしくは非公表の場合には「不明」と記載。
      Confirmation DatePCR検査における陽性が判明した日。 日付は日本時間(JST)。自治体の公表日(報道発表資料公開日)や、報道された日付ではない。
      演習用データ COVID-19_modified.csv の定義(引用:新型コロナウイルス感染者数マップ

      各列のデータの表示

      読み込んだデータの型がdata.frameという形式の場合,「オブジェクト名$列名」と指定することで,各列のデータのみを表示できる.今回は,既にdata.frame型になっている.たとえば,居住都道府県(居住地)のみのデータを表示したければ,以下のように入力すればよい.

      COVID19$Residential.Pref

      ◆データの要約

      膨大なデータの概観を得るためには,データの特徴を少数の値に要約する(まとめる)ことがもとめられる.ここでは,table関数,summary関数を用いたデータの要約をおこなう.

      問1 table関数を用いた都道府県別感染者数の集計

      COVID-19は,1つのレコードに感染者1名分のデータが入っている.これを,各列ごとに集計してみよう.まずは,下記のコードを入力してみてほしい.どのような結果が出てきただろうか?

      table(COVID19$Residential.Pref)
      ?table

      実行結果を踏まえて,以下の1,2に答えよ.

      1. Age,Hospital Pref,Residential Pref,Genderの4つの列について,table関数による集計を実行し,そのコードを示せ.
      2. 「?table」は,table関数がどのような関数を表示するコマンドである.これを参考に,table関数がどのような関数か,説明せよ.

      問2 summary関数を用いた要約

      下記のコードを実行せよ.

      summary(COVID19)

      列ごとに,「Min.」「1st Qu.」「Median」「Mean」「3rd Qu.」「Max.」「NA’s」という項目が表示される.これは,上から順に「最小値」「第 1 四分位点」「中央値」「平均」「第 3 四分位点」「最大値」「NAの数」を意味する.「NA」は,データの欠損を意味する.また,数値データ以外の列は,「Length」「Class」「Mode」となっており,それぞれ「データの数(行数),データ型,データ型」を意味する.

      summary関数の実行結果を踏まえて,以下の1,2に答えよ.

      1. 今回演習で扱うCOVID19データには何人分のデータが入っているか.
      2. Ageの列の平均値,最小値,最大値を求めよ.
      3. 今回Ageの列において,欠損しているデータはいくつあるか.

      ◆量的データの可視化

      膨大なデータの概観を把握するもう一つの手法は「可視化」である.量的データの可視化法として最も代表的なのがヒストグラム(度数分布)である.度数分布はhist関数を用いて表示することができる.

      問3 hist関数を用いた年代データの可視化

      年代に関する度数分布を表示する下記のコードを実行せよ.

      hist(x = COVID19$Age,
           xlab = "XXX",
           ylab = "XXX",
           main = "XXX",
           col = "red",
           border = "gray",
           breaks=seq(0,100,by=10))

      各引数の意味は以下の通りである.

      xヒストグラムを描く変数(ベクトル形式)
      breaksヒストグラムの横軸の範囲と分割の幅(ベクトル形式)
      xlab,横軸の名称(文字列)
      ylab縦軸の名称(文字列)
      mainヒストグラム上のタイトル(文字列)
      colヒストグラムの棒の色を指定(文字列または数値)
      border棒の枠の色を指定(文字列または数値)
      breaksベクトル形式で,階級分けの下端値から上端値の値を指定する.例えばbreaks=seq(0,100,by=10)は「0歳代から100歳代までの年齢を10歳刻みのヒストグラムに分けよ」という意味である.
      各引数の意味

      実行結果を踏まえて,以下の問いに答えよ.

      1. 縦軸に「頻度」,横軸に「年代」,タイトルに「年代のヒストグラム」と指定し,棒とその枠の色を「black」としたヒストグラムを描け.また,グラフは,形式を「PNG」にして,EIP2020フォルダに保存し,レポートに貼りつけること(図番号とキャプションもつけること).
      2. 一番感染者数が多い年代と,少ない年代を答えよ.また,描いたヒストグラムからわかることを説明せよ.

      問4 hist関数を用いた陽性判明日データの可視化

      陽性が判明した日に関する度数分布を表示する下記のコードを実行せよ.

      hist(x = COVID19$Confirmation.Date,
           breaks = length(unique(COVID19$Confirmation.Date)),
           xlab = "XXX",
           ylab = "XXX",
           main = "XXX",
           col = "gray",
           freq = TRUE)

      上記を実行した結果、以下のようなエラー文「’x’ must be numeric(オブジェクト’x’は実数型でなければいけない)」が出る場合がある。これは、COVID19$Confirmation.Dateというオブジェクトが、単なる文字列として扱われているためである。この場合、

      > hist(x = COVID19$Confirmation.Date,
      +      breaks = length(unique(COVID19$Confirmation.Date)),
      +      xlab = "XXX",
      +      ylab = "XXX",
      +      main = "XXX",
      +      col = "gray",
      +      freq = TRUE)
      Error in hist.default(x = COVID19$Confirmation.Date, breaks = length(unique(COVID19$Confirmation.Date)),  : 
        'x' must be numeric

      hist()関数を実行する「前」に以下のコードを実行する。as.Date()は、()内の日付に関するオブジェクトをcharacter型からDate型へ変換する関数である。変換した後の結果を、オブジェクトCOVID19$Confirmation.Dateに改めて代入する。

      COVID19$Confirmation.Date <- as.Date(COVID19$Confirmation.Date)

      hist関数の実行結果を踏まえて,以下の問いに答えよ.

      1. 縦軸,横軸,タイトルに適当(データを的確に表すよう)な文字列が入るよう,引数を指定したヒストグラムを描け.また,グラフは,形式を「PNG」にして,EIP2020フォルダに保存し,レポートに貼りつけること(図番号とキャプションもつけること).
      2. summary関数を用いて,いつからいつまでの期間のデータが収録されているかを示せ.
      3. 陽性が判明した日について,ヒストグラムからわかることを3つ以上挙げよ.

      ◆質的データの可視化

      質的データを可視化するには,各カテゴリに属する割合を図示するのが良い.ここでは円グラフを扱う.円グラフは,pie関数を用いて表示することができる.

      5 pie関数を用いた性別データの可視化

      性別に関する度数分布を表示する下記のコードを実行せよ.

      x <- table(COVID19$Gender)
      pie(x, clockwise = TRUE)

      第一引数には,各カテゴリの個体数をベクトルで指定する必要がある.第二引数のclockwiseには,カテゴリの順序を反時計回りに配置するかをTRUE,FALSEで指定する.

      実行結果を踏まえて,以下の問いに答えよ.

      1. グラフを描き,形式を「PNG」にして,EIP2020フォルダに保存し,レポートに貼りつけること(図番号とキャプションもつけること).
      2. 今回は「x <- table(COVID19$Gender)」を実行してから,「pie(x, clockwise = TRUE)」を実行している.この二行がどのような操作を意味するのか,なぜ一行目が必要であるのかを説明せよ.
      3. 円グラフからわかることを説明せよ.

      ◆データの抽出

      ここまでは,全国のデータを扱ってきたが,愛知県のデータだけを取り出して,そのデータの特性を見る.

      下記のコードを入力してみてほしい.急にややこしいのが出てきた,と思うかもしれないが,一つ一つ丁寧に見てみよう.

      AichiCOVID19 <- COVID19[which(COVID19$Residential.Pref=="Aichi"),]

      まず,以下のコードを実行してみてほしい.FALSEやTRUEという文字列が現れる.これは,COVID19$Residential.Prefの値が”Aichi”であればTRUE,そうでなければ”FALSE”を返している.

      COVID19$Residential.Pref=="Aichi"

      次に,以下のコードを実行してみる.which関数である.これは,which( )の中の論理式を満たす行の番号を抽出する関数である.つまりCOVID19$Residential.Prefの値が”Aichi”である行の番号一覧が返されていることになる.

      which(COVID19$Residential.Pref=="Aichi")

      最後に,以下のコードを実行してみる.以下のようにオブジェクト名[ ]と表記しているのは,[ ]の数字を使ってオブジェクトの何番目のデータを取り出すかを指定するものである.今回は行列データを扱っているのでオブジェクト名[i,j]でi行j列のデータを指定することが出来る.iもjも,ベクトルで指定することが出来る.which関数の結果はベクトル形式になっているから,そのベクトルの数値に該当する行のデータを返すことになる.

      COVID19[which(COVID19$Residential.Pref=="Aichi"),]

      問6 データの抽出

      全国のデータから岐阜県(”Gifu”)のデータだけを取り出し,オブジェクト名「GifuCOVID19」に付値せよ.また,そのコードを示せ.

      ◆発展課題(+10点)

      以下の問題は,できる範囲で取り組んでください.

      1. 好きな都道府県(愛知,岐阜,あるいはそれ以外)を選び,hist関数を用いた年代データと陽性判明日データ,pie関数を用いた性別データの可視化をおこない,保存し,レポートに張り付けよ.
      2. 結果について,全国を対象に分析した問1~問5の結果と比較し,考察しなさい.

      演習作業用フォルダの準備方法

      書いたプログラムや分析するデータを保存するフォルダを用意しておきます.演習の際には,このフォルダをワーキングディレクトリ(作業用フォルダ)として指定します.

      自分がどのフォルダで作業しているかを意識することは非常に大事です.適当なところに保存して「保存したプログラムがどこにあるかわからない,消えた」といったことにならないように,最後の確認までしっかりやりましょう.

      ◆演習作業用フォルダ(ワーキングディレクトリ:WD)の準備

      デスクトップ上で右クリック→「新規作成」→「フォルダー」から,新しいフォルダを作成し,「EIP2021」という名前にしておきます.

      EIP2021フォルダを右クリック→プロパティを選択すると,「場所:XXXXX」という項目が見られます.このパスの末尾に「\EIP2021」を足した文字列で示されるパスをWDに指定することにします.

      フォルダのプロパティ

      この場合のWDは以下の通りです.今後,演習で使うファイル類はすべてこれに入れてください.後でまとめて提出してもらいます.

      C:\Users\ユーザー名\Desktop\EIP2021

      ◆Rでワーキングディレクトリを指定する

      R上で,getwd()と入力して実行してみてください.以下のような値が返ってきます.getwd関数は,現在のワーキングディレクトリ(CWD: Current Working Directory)を返す関数です.

      > getwd()
      [1] "D:/documents"

      CWDを,先ほど作成した演習作業用フォルダ「C:\Users\ユーザー名\Desktop\EIP2021」に変更します.※ユーザー名は各自のフォルダのパスに依存しますので注意してください.フォルダの階層の区切りは「/」を使います.

      setwd("C:/Users/Fuko Nakai/Desktop/EIP2021")

      もう一度,getwd()と入力して実行してみてください.以下のような値が返ってきたら成功です.

      > getwd()
      [1] "C:/Users/Fuko Nakai/Desktop/EIP2021"

      ◆演習用エディタの立ち上げとスクリプトの保存

      「RConsole」と書かれたウィンドウ(コンソール画面)が現れていることを確認します.ここに直接コマンドを入力することもできますが,ここに打ち込んだ内容は,Rを終了すると消えてしまいます.そこで「Rエディタ」を立ち上げ,スクリプトファイルを新規作成し,保存することにします.

      Rstudioの場合

      左上に「Untitled」と書かれたタブ(Rエディタ)が出ていることを確認して,左上のタブの「File」→「Save as」をクリックします.

      先ほど作成したEIP2021内に「yyyymmdds_Name.R」という名前で保存します.「.R」というのはR拡張子です.ファイル名は何でも良いですが,後で誰が見ても内容のわかる内容にしましょう.yyyymmddには作成年月日を,Nameには自分の名前アルファベット表記で入れてください.

      保存したら,EIP2021の中に「yyyymmdds_Name.R」が出現していることを確認しましょう.演習のプログラムはこのファイル上に書き込み,上書き保存してください.

      図1 Rstudioの画面

      R Guiの場合

      メニューバーから「ファイル(File)」→「新しいスクリプト」を選択します.その後,下のような「無題 – Rエディタ」がでます.「ファイル」→「保存」で保存の操作ができます.

      Rの場合の表示画面

      以上で作業は終了です.