課題名
rm(list = ls()) # 全オブジェクト削除
# ここにRコードを記入する。
#例）
plot(1:9)
d <- read.csv(file = 'https://stats.dip.jp/01_ds/data/dirty_data2.csv')
library(DT)
datatable(d)
d1 <- na.omit(d)
summary(d1) # 要約統計値 NAがなくなった。
## 会社名 気温 湿度 風速
## Length:18 Min. :12.00 Min. :-50.00 Min. :-1.000
## Class :character 1st Qu.:15.50 1st Qu.: 30.00 1st Qu.: 1.000
## Mode :character Median :32.40 Median : 50.00 Median : 2.000
## Mean :28.05 Mean : 38.89 Mean : 2.228
## 3rd Qu.:34.11 3rd Qu.: 51.00 3rd Qu.: 4.000
## Max. :50.40 Max. : 60.00 Max. : 5.000
## 全天日射量MJ
## Min. :-1.200
## 1st Qu.: 0.340
## Median : 1.215
## Mean : 2.791
## 3rd Qu.: 2.250
## Max. :23.400
datatable(d1)
# 各気象変数の許容範囲に入っているか否かを示す論理式を作成
is.tp <- -20 <= d1$気温 & d1$気温 <= 50
is.hm <- 0 <= d1$湿度 & d1$湿度 <= 100
is.ws <- 0 <= d1$風速 & d1$風速 <= 70
is.mj <- 0<= d1$全天日射量 & d1$全天日射量 <= 4.00
summary(is.tp)
## Mode FALSE TRUE
## logical 2 16
summary(is.hm)
## Mode FALSE TRUE
## logical 1 17
summary(is.ws)
## Mode FALSE TRUE
## logical 1 17
summary(is.mj)
## Mode FALSE TRUE
## logical 4 14
# 適切な範囲のレコードだけを抽出する。
d2 <- d1[is.tp & is.hm & is.ws & is.mj,]
datatable(d2)
d3 <- d2
d3$気温 <- round(d2$気温, 1)
d3$湿度 <- round(d2$湿度, 0)
d3$風速 <- round(d2$風速, 0)
d3$全天日射量 <- round(d2$全天日射量, 2)
datatable(d3)
#d3$会社名
# ソニー株式会社を意味する表記
#LAB.SONY <- c("東芝","Toshiba","TOSHIBA","toshiba","東芝株式会社","東芝（株） " , "東芝(株)" )
# 富士通株式会社を意味する表記
#LAB.FUJITSU <- c("トヨタ","トヨタ)
d <- read.csv(file = 'https://stats.dip.jp/01_ds/data/gender.csv')
library(DT)
datatable(d, caption = '性別データ')
d$性別 <- ifelse(d$性別 == '男性', 1, 0)
x <- d$身長
y <- d$性別
fit <- glm(y ~ x, data = d, family = 'binomial')
#summary(fit)
library(sjPlot)
## #refugeeswelcome
tab_model(fit, show.aic = T)
|
|
Dependent variable
|
|
Predictors
|
Odds Ratios
|
CI
|
p
|
|
(Intercept)
|
0.00
|
0.00 – 0.00
|
<0.001
|
|
x
|
1.84
|
1.53 – 2.37
|
<0.001
|
|
Observations
|
135
|
|
R2 Tjur
|
0.849
|
|
AIC
|
45.452
|
matplot(x = x, y = y, pch = 1,
main = 'ロジスティック回帰分析',
xlab = '身長',
ylab = '男性(1)，女性(0)')
grid()
xp <- seq(min(d$身長), max(d$身長),1)
yp <- predict(fit, type = 'response', newdata = data.frame(x = xp))
matlines(x = xp, y = yp, lty = 1, col = 2, lwd = 2)
library(latex2exp)
text(x = 3.8, y = 0.4, adj = 0,
TeX('$\\leftarrow \\hat{p}=\\frac{1}{1+e^{-(b_0 + b_1 x)}}$'))
legend('topleft', pch = c(1, NA), lty = c(NA, 1), col = c(1, 2),
legend = c('男性(1)，女性(0)', 'ロジスティック回帰モデル'))
yp <- predict(fit, type = 'response', newdata = data.frame(x = 170))
sprintf('男性である確率が%2.1f％であるため女性であると予測される', yp * 100)
## [1] "男性である確率が17.2％であるため女性であると予測される"