回帰分析をRで実施する。Bostonデータセットを使って実践。
その後、train、validation、testセットに分けて、MSEを計算する。
library(MASS) # Boston データセットを使う
library(tidyverse) # ggplot2とdiplyrを使う線形回帰分析 Regression 重回帰・単回帰
以下の形で、回帰分析のオブジェクトを作る。
mylm <- lm(data=データフレーム, outcome ~ predictor_1 + predictor_2)
outcomeは目的変数y、predictor_1は説明変数１、predictor_2は説明変数２とする。
今回は、ＭＡＳＳの中にあるＢostonデータセットを使用する。Bostonの中には、変数medv(median value of owner-occupied homes in $1000s)と変数lstat(lower status of the population (percent).)がある。
medvをｙとして、lstatをｘとして式を定義する。このときに、Boston\(medv ~ Boston\)lstat とすると、うまくいかない。
mylm <- lm(data=Boston, medv ~ lstat)coef()を使うと、Interceptとcoefficientsを得ることができる。
coef(mylm)## (Intercept) lstat
## 34.5538409 -0.9500494summary() を使うと、Multiple R-squared、Adjusted R-squared、Intercept、coefficients等など、様々な情報を得ることができる。
summary(mylm)##
## Call:
## lm(formula = medv ~ lstat, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.168 -3.990 -1.318 2.034 24.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.55384 0.56263 61.41 <2e-16 ***
## lstat -0.95005 0.03873 -24.53 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared: 0.5441, Adjusted R-squared: 0.5432
## F-statistic: 601.6 on 1 and 504 DF, p-value: < 2.2e-16predict()を使うと、さきほどの回帰分析のモデルを使って目的変数を予測することできる。
predict(回帰モデル,説明変数)
これで得られるものは、目的変数を予想したもの。 特に意味はないが、得られた回帰モデルを使って、説明変数から目的変数を予測してみる。
predicted_value <- predict(mylm,Boston[,13,drop=F])
head(predicted_value)## 1 2 3 4 5 6
## 29.82260 25.87039 30.72514 31.76070 29.49008 29.60408以下のように説明変数にdrop=Fが必要なのは、説明変数がデータフレームである必要があるから。 Boston$lstatだと、ベクターになってしまう。
新たな説明変数を使って、予測してみたい。列の名前は、モデルの説明変数の名前と同じにしなければならない。
pred_dat <- data.frame(seq(1, 40, length=1000))
names(pred_dat) <- "lstat"
y_pred_new <- predict(mylm, pred_dat)
head(y_pred_new)## 1 2 3 4 5 6
## 33.60379 33.56670 33.52961 33.49252 33.45544 33.4183595%信頼区間を得る方法。
y_pred_95 <- predict(mylm, newdata = pred_dat[,1,drop=F], interval = 'confidence')
head(y_pred_95)## fit lwr upr
## 1 33.60379 32.56402 34.64356
## 2 33.56670 32.52947 34.60394
## 3 33.52961 32.49491 34.56432
## 4 33.49252 32.46035 34.52470
## 5 33.45544 32.42578 34.48509
## 6 33.41835 32.39122 34.44547グラフにしたいので、説明変数の列を加える。
y_pred_95 <- predict(mylm, newdata = pred_dat[,1,drop=F], interval = 'confidence')
y_pred_95 <- data.frame(y_pred_95,pred_dat[,1,drop=F])
head(y_pred_95)## fit lwr upr lstat
## 1 33.60379 32.56402 34.64356 1.000000
## 2 33.56670 32.52947 34.60394 1.039039
## 3 33.52961 32.49491 34.56432 1.078078
## 4 33.49252 32.46035 34.52470 1.117117
## 5 33.45544 32.42578 34.48509 1.156156
## 6 33.41835 32.39122 34.44547 1.195195分かりやすいように名前を変える
names(y_pred_95) <- c("medv","lwr","upr","lstat")ggplot() +
geom_ribbon(data=y_pred_95,aes(x=lstat,ymin=lwr,ymax=upr),color="blue",fill="blue") +
geom_point(data=Boston,aes(x=lstat,y=medv)) +
geom_line(data=y_pred_95,aes(x=lstat,y=medv),color="green") +
theme_classic()
MSE(Mean Squared Error)の関数を自作する。これは、予測値と実際の数値の違いを足し合わせたもの。
mse <- function(y_true, y_pred) {
vec_length = length(y_true)
mse_value =0
for (i in 1:vec_length){
mse_value = mse_value + (y_true[i]-y_pred[i])^2
}
mse_value = mse_value/vec_length
return(mse_value)
}さきほど予測したmedvの値と、実際のmedvの値についてMSEを求める。
mse(Boston$medv,predicted_value)## 1
## 38.48297Bostonのデータセット502行を、 train データセット（５割）と、validationデータセット（３割）と、testデータセット（２割）に分別する。
splits <- c(rep("train",253),rep("validation",152),rep("test",101))
head(splits)## [1] "train" "train" "train" "train" "train" "train"Bostonのデータセットに、上で作った１列を追加する。名前は、boston_masterとする。
boston_master <- Boston %>%
mutate(splits=sample(splits))boson_masterから、train,validation,test に分ける。
boston_train <- boston_master %>% filter(splits=="train")
boston_valid <- boston_master %>% filter(splits=="validation")
boston_test <- boston_master %>% filter(splits=="test")trainセットのみで、regressionを学習させる。
model_1 <- lm(data=boston_train, medv~lstat)
summary(model_1)##
## Call:
## lm(formula = medv ~ lstat, data = boston_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.137 -4.350 -1.369 1.784 23.615
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.12775 0.81584 44.28 <2e-16 ***
## lstat -1.09719 0.05916 -18.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.363 on 251 degrees of freedom
## Multiple R-squared: 0.5781, Adjusted R-squared: 0.5765
## F-statistic: 344 on 1 and 251 DF, p-value: < 2.2e-16train内のMSEを計算する。
pred_train <- predict(model_1,boston_train[,13,drop=F])
model_1_mse_train <- mse(boston_train$medv,pred_train)
model_1_mse_train## 1
## 40.16525trainで学習させたモデルを使って、validationを予測する。そして、MSEを計算する。
pred_valid <- predict(model_1, boston_valid[,13,drop=F])
model_1_mse_valid <- mse(boston_valid$medv,pred_valid)
model_1_mse_valid## 1
## 33.58956他の数式を使って、regressionを学習させる。
model_2 <- lm(data=boston_train, medv~age+tax)
pred_train2 <- predict(model_2,boston_train[,c(7,10),drop=F])
model_2_mse_train <- mse(boston_train$medv,pred_train2)
print(paste("model_2_mse_train:",as.character(model_2_mse_train)))## [1] "model_2_mse_train: 75.5125419055823"pred_valid2 <- predict(model_2,boston_valid[,c(7,10),drop=F])
model_2_mse_valid <- mse(boston_valid$medv,pred_valid2)
print(paste("model_2_mse_valid:",as.character(model_2_mse_valid)))## [1] "model_2_mse_valid: 40.5666012367619"数式によるMSEの比較してみる。
print(paste("model_1_mse_train:",as.character(model_1_mse_train)))## [1] "model_1_mse_train: 40.1652537382155"print(paste("model_1_mse_valid:",as.character(model_1_mse_valid)))## [1] "model_1_mse_valid: 33.5895601321273"print(paste("model_2_mse_train:",as.character(model_2_mse_train)))## [1] "model_2_mse_train: 75.5125419055823"print(paste("model_2_mse_valid:",as.character(model_2_mse_valid)))## [1] "model_2_mse_valid: 40.5666012367619"MSEはモデル1の方がモデル2よりも低かったので、モデル1を選択する。
今度は、最終的にtestセットでMSEを計算する。
pred_test <- predict(model_1,boston_test[,13,drop=F])
model_1_mse_test <- mse(boston_test$medv,pred_test)
print(paste("model_1_mse_test:",as.character(model_1_mse_test)))## [1] "model_1_mse_test: 47.5681944974811"K-fold cross validationをRで実装する。
kfold_validate <- function(formula,dataset,k){
if (k>1){
d_length <- nrow(dataset)
indexs <- rep(1:k, length.out = d_length)
indexs <- data.frame(indexs)
names(indexs) <- "indexs"
outcome <- all.vars(formula)[1]
predictors <- all.vars(formula)[-1]
MSE_value = 0
for (i in 1:k){
new_data <- cbind(dataset,indexs)
dataset_train <- new_data %>% filter(indexs != i)
dataset_valid <- new_data %>% filter(indexs == i)
mylm <- lm(data = dataset_train,formula)
pred <- predict(mylm,dataset_valid[,predictors])
MSE_value = MSE_value + mse(dataset_valid[,outcome],pred)
}
}else{
print("k should be higher 1")
}
return(MSE_value/k)
}解説 outcome <- all.vars(formula)[1] predictors <- all.vars(formula)[-1]
は、以下の式から左辺と右辺をとりだしている。 formula <- outcome ~ predictor_1 + predictor_2
dataset_valid[,outcome] は、上記の左辺と同じ列を取り出している。
dataset_valid[,predictors] は、上記の右辺と同じ列を取り出している。
試しに、9-fold cross validationを実行してみる。下のように、２つのモデルを比べてみる。 medv ~ lstat + age + tax medv ~ lstat + I(lstat^2) + age + tax
boston_train_validation <- boston_master %>% filter(splits=="train"|splits=="validation")
formula1 <- medv ~ lstat +age +tax
results1 <- kfold_validate(formula1,boston_train_validation,9)
print(paste("MSE_k-fold cross-validation by formula1:",results1))## [1] "MSE_k-fold cross-validation by formula1: 36.4827657175701"mylm <- lm(data = boston_test,formula1)
pred <- predict(mylm,boston_test[,c("lstat","age","tax")])
print(paste("MSE test set by formula1:",mse(boston_test[,"medv"],pred)))## [1] "MSE test set by formula1: 40.4491711116041"formula2 <- medv ~ lstat + I(lstat^2) + age + tax
results2 <- kfold_validate(formula2,boston_train_validation,9)
print(paste("MSE_k-fold cross-validation by formula2:",results2))## [1] "MSE_k-fold cross-validation by formula2: 28.1007057996556"mylm <- lm(data = boston_test,formula2)
pred <- predict(mylm,boston_test[,c("lstat","age","tax")])
print(paste("MSE test set by formula2:",mse(boston_test[,"medv"],pred)))## [1] "MSE test set by formula2: 26.3088757875499"Formula 2 のほうがＭＳＥ的にいいスコアを出した。