ロジスティック曲線
kagaya fumitaka
2020/2/28
# 結論：ロジスティック回帰は本解析には向かない
参考資料 https://to-kei.net/r-beginner/r-5-logistic/
ロジスティック回帰と変数選択 https://oku.edu.mie-u.ac.jp/~okumura/stat/140921.html
「ロジスティック回帰」は機械学習モデルか？ https://assign-navi.jp/magazine/consultant/data-scientist/c033.html
-ロジスティック回帰は機械学習か？
機械学習とは 最近、機械学習という言葉がIT業界でもよく用いられるようになったと感じています。
名前だけ見ると、工場にあるような機械が学習するような印象がありますが、英語のMachine Learningの直訳表現に過ぎません。英語のMachineにはコンピュータやサーバなどの意味もあり、Machine Learningはコンピュータに学ばせる、コンピュータ学習などといった意味になります。
ではコンピュータに何を学ばせるのでしょうか。それは、データに潜む構造的なパターンや特性です。ただし、構造そのものを学習してくれるわけではありません。あくまで、ある構造を仮定したときに、その構造を決定する特性値を学習します。
例で示しましょう。前回の線形回帰分析も、見方によれば機械学習のひとつといえます。線形回帰分析の構造を仮定しなければなりません。
その構造とは、「(正規分布に従うと考えられる)数値データが、他の変数の一次結合(ある値を掛けて足し合わせること)によって表現される」というものです。では、この構造を決定する特性値とは、何でしょうか。それは、一次結合する際に変数に掛ける値となります。
前回のコラムでは、最高気温とアイスクリームの売上の関係について解説しました。売上は最高気温と直線関係にあることから、次のような式(つまり構造)を仮定したわけです。
(アイスクリームの売上)＝(特性値1)×(最高気温)円 + (特性値2)円 最高気温という変数に特性値1を掛けたものと、特性値2を足し合わせています。なお、特性値2は一般的に切片と呼ばれますが、これは“1“という変数に特性値2が掛けられたものと捉えてください。こうすることで、既存の変数では説明できない特徴を表現することができます。
構造を決めたら、あとはデータが構造に最も当てはまるような特性値1と特性値2を、データを用いてコンピュータが学習するというわけです。コンピュータが特性値を学習すると、データが持つ構造(数式)が決まります。
また、新しいデータが生成されたときにその構造に当てはめることで、結果を予測することができるというわけです。
この考え方が機械学習の本質となります。世の中のデータに潜む構造には様々なものがあります。
はじめに人間が構造を仮定する必要があることから、これまで多くの構造が考案されてきました。それらは一般的にはモデルと呼ばれています。線形回帰モデル、ロジスティック回帰モデル、ニューラルネットワークモデルなど、現在では多くのモデルが存在します。また、モデルを決定する特性値のことをパラメータと呼びます。
広義の機械学習は上記の理解で良いのですが、狭義には線形回帰モデルは機械学習に含まれません。実は線形回帰モデルはパラメータがデータの関数として一意に決まります(データを決まった数式に入れるだけでパラメータが求まってしまいます)。
一方、ロジスティック回帰の場合は、データを入れるだけでパラメータが簡単に求まるような数式はありません。そのため、まず初めにランダムなパラメータでモデルを作り、そのモデルによる予測値と実際のデータとの誤差を計算します。
そして、少しずつパラメータを変化させ、誤差が最も小さくなるようなパラメータを探索します。パラメータの探索は、誤差が小さくなるように探索する方法のほか、実際のデータが得られる確率が最も大きくなるような方法もあります。
このパラメータ探索のことを“学習”と呼ぶことがあるため、パラメータ探索をしていない線形回帰モデルは機械学習には含まれないとする見方があります。
- ロジスティック回帰とは
# ロジスティック回帰による病気分析
logi_fun <- function(data,file,disease){
ans <- glm(data$Y~.,data=data,family=binomial)
s.ans <- summary(ans)
coe <- s.ans$coefficient
RR <- exp(coe[,1])
RRlow <- exp(coe[,1]-1.96*coe[,2])
RRup <- exp(coe[,1]+1.96*coe[,2])
N <- nrow(data)
aic <- AIC(ans)
result <- cbind(coe,RR,RRlow,RRup,aic,N)
colnames(result)[6:7] <- c("RR95%CI.low","RR95%CI.up")
if(nrow(result)>=2){
result[2:nrow(result),8:9] <- ""
}
write.table(disease,file,append=T,quote=F,sep=",",row.names=F,col.names=F)
write.table(matrix(c("",colnames(result)),nrow=1),file,append=T,quote=F,sep=",",row.names=F,col.names=F)
write.table(result,file,append=T,quote=F,sep=",",row.names=T,col.names=F)
write.table("",file,append=T,quote=F,sep=",",row.names=F,col.names=F)
}df <- read.csv("sample-data.csv",header=T,row.names=1)dat <- df[,c(5,1,2,6)]#５列が目的変数。1,2,6列目が説明変数。colnames(dat)[1] <- "Y"結果は、C:\721540~-data.csvにある。
#その2 http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/15_logit.html
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("ggplot2")
library("reshape")## Warning: package 'reshape' was built under R version 3.6.2##
## Attaching package: 'reshape'## The following object is masked from 'package:dplyr':
##
## renamelibrary("reshape2")##
## Attaching package: 'reshape2'## The following objects are masked from 'package:reshape':
##
## colsplit, melt, recastlibrary("stargazer")##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazer## ロジスティック関数を定義する
## 引数：x = 数値ベクトル
## 返り値：y
invlogit <- function(x){
return( y <- 1 / (1 + exp(-x)) )
}- ロジスティック曲線を書く http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/15_logit.html
x <- seq(-5, 5, length = 100)
plot(x, invlogit(x), type = "l", lwd = 2, col = "darkblue",
ylab = "invlogit(x)", main = "Logistic Function")
abline(h = c(0, 1), col = "gray")
- ロジスティック関数による当落判定の実行 - data：logit.csv→C:\721540~
# データフレーム logit に含まれている変数を表示する。
logit <- read.csv("logit.csv", na = ".")
names(logit)## [1] "wlsmd" "previous" "expm"# データフレーム logit を表示する。
logit## wlsmd previous expm
## 1 1 0 10.0
## 2 1 1 10.0
## 3 1 3 8.9
## 4 1 5 7.7
## 5 1 7 5.4
## 6 1 4 3.0
## 7 1 5 2.0
## 8 1 4 7.0
## 9 0 0 2.2
## 10 0 0 4.3
## 11 0 5 4.0
## 12 0 2 3.0
## 13 0 4 4.0
## 14 0 1 2.2
## 15 0 3 4.0-2.2 帰無仮説と対立仮説の設定 -ここでの帰無仮説は次のとおり。
・H01：当選回数は小選挙区での当選確率に影響を与えない ・H02：選挙費用は小選挙区での当選確率に影響を与えない
-対立仮説は次のとおり。
・Ha1：当選回数が多いほど、小選挙区での当選確率が高くなる ・Ha2：選挙費用が多いほど、小選挙区での当選確率が高くなる
-2.3 説明変数と応答変数の散布図 ・ダウンロードしたデータ logit を使い、候補者の小選挙区における当落 (wlsmd) を縦軸に、それまでの当選回数 (previous) と選挙費用 (expm) を横軸にした散布図を表示する。
-以下の二つのグラフを見ると、当選回数と当落は関係が低く、選挙費用と当落は大いに関係がある。
# 当選回数 (previous) と当落 (wlsmd)の散布図
ggplot(logit, aes(previous, wlsmd)) + geom_point() +
stat_smooth(method = lm, se = FALSE) +
theme_bw() +
geom_jitter(width = 0.05, height = 0.05) #jitterで重複したデータを散らす
# 選挙費用 (expm) と当落 (wlsmd)の散布図
ggplot(logit, aes(expm, wlsmd)) + geom_point() +
stat_smooth(method = lm, se = FALSE)+
theme_bw() +
geom_jitter(width = 0.05, height = 0.05)#jitterで重複したデータを散らす
-ロジスティック回帰式の推定値を求める
・ロジスティック回帰分析を行うには glm() 関数を利用する。 ・ glm() 関数はロジスティク回帰以外にも様々な一般化線形モデルで利用される。 ・ロジスティック回帰に使うには、family を以下のように指定する。
model_1 <- glm(wlsmd ~ previous + expm, data = logit,family = binomial(link = "logit"))
summary(model_1)##
## Call:
## glm(formula = wlsmd ~ previous + expm, family = binomial(link = "logit"),
## data = logit)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5741 -0.3781 0.2013 0.3943 1.4948
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.3811 3.5147 -1.816 0.0694 .
## previous 0.8085 0.5851 1.382 0.1670
## expm 0.8088 0.4000 2.022 0.0431 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 20.728 on 14 degrees of freedom
## Residual deviance: 10.384 on 12 degrees of freedom
## AIC: 16.384
##
## Number of Fisher Scoring iterations: 6- 当選確率の予測値 pを求める。 -・Rでは predict() 関数を使うと、予測値を計算してくれる。
predict(model_1, type = "response")## 1 2 3 4 5 6
## 0.846455457 0.925228189 0.962419494 0.979953464 0.974574115 0.327276800
## 7 8 9 10 11 12
## 0.327213998 0.925168964 0.009934946 0.051995861 0.710296138 0.088057029
## 13 14 15
## 0.522058216 0.022027721 0.327339608- 1から8までが当選した実績があり、9から15までが落選した実績がある。 - これからすると11を除き、0.3以上と以下に分かれ、概ね当落実績どおりの結果となっている。
stargazer(model_1, style = "ajps", # AJPS仕様
dep.var.labels = "Win or Lose", # 応答変数の名前
title = "Logistic Regression Analysis on Winning a Seat", # タイトル
type = "html")| Win or Lose | |
| previous | 0.809 |
| (0.585) | |
| expm | 0.809** |
| (0.400) | |
| Constant | -6.381* |
| (3.515) | |
| N | 15 |
| Log Likelihood | -5.192 |
| AIC | 16.384 |
| p < .01; p < .05; p < .1 | |
-・previousは有意ではないが、expm は 5%で統計的に有意であることがわかる。
-・expmの係数 (0.809)の意味: 候補者の選挙費用が 100万円増える → 当選する log odds が 0.809 増える ・これだとわかりにくい
-ひとつの解決策: ・log odds → predict()関数を使って当選確率 に変換 ・二つの説明変数（previous と expm) の値を変えて当選確率を予測す
厚木市下水道管路で試行する
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(caret)## Warning: package 'caret' was built under R version 3.6.2## Loading required package: latticelibrary(cluster)
library(dummies)## dummies-1.5.6 provided by Decision Patternslibrary(data.table)## Warning: package 'data.table' was built under R version 3.6.2##
## Attaching package: 'data.table'## The following objects are masked from 'package:reshape2':
##
## dcast, melt## The following object is masked from 'package:reshape':
##
## melt## The following objects are masked from 'package:dplyr':
##
## between, first, lastlibrary(dplyr)
library(e1071)## Warning: package 'e1071' was built under R version 3.6.2library(epitools)
library(effects)## Warning: package 'effects' was built under R version 3.6.3## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car## Use the command
## lattice::trellis.par.set(effectsTheme())
## to customize lattice options for effects plots.
## See ?effectsTheme for details.library(ggplot2)
library(ggthemes)## Warning: package 'ggthemes' was built under R version 3.6.2library(randomForest)## Warning: package 'randomForest' was built under R version 3.6.3## randomForest 4.6-14## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:ggplot2':
##
## marginlibrary(ranger)## Warning: package 'ranger' was built under R version 3.6.2##
## Attaching package: 'ranger'## The following object is masked from 'package:randomForest':
##
## importancelibrary(rgl)
library(rattle)## Warning: package 'rattle' was built under R version 3.6.2## Rattle: A free graphical interface for data science with R.
## バージョン 5.3.0 Copyright (c) 2006-2018 Togaware Pty Ltd.
## 'rattle()' と入力して、データを多角的に分析します。##
## Attaching package: 'rattle'## The following object is masked from 'package:ranger':
##
## importance## The following object is masked from 'package:randomForest':
##
## importancelibrary(readr)## Warning: package 'readr' was built under R version 3.6.2library(rpart.plot)## Warning: package 'rpart.plot' was built under R version 3.6.2## Loading required package: rpartlibrary(rpart)
library(readr)
library(reshape)
library(rsconnect)## Warning: package 'rsconnect' was built under R version 3.6.2library(reshape2)
library(tidyr)##
## Attaching package: 'tidyr'## The following object is masked from 'package:reshape2':
##
## smiths## The following objects are masked from 'package:reshape':
##
## expand, smithslibrary(xtable)
library(nnet)## Warning: package 'nnet' was built under R version 3.6.2library(stargazer)Rによる条件付きロジスティック回帰分析
https://to-kei.net/r-beginner/r-5-logistic/
ロジスティック関数を R で定義してみる。
## ロジスティック関数を定義する
## 引数：x = 数値ベクトル
## 返り値：y
invlogit <- function(x){
return( y <- 1 / (1 + exp(-x)) )
}確認のため、(−5<x<5) の範囲で logistic(x)=logit(x−1)) の値を図示する。
x <- seq(-5, 5, length = 100)
plot(x, invlogit(x), type = "l", lwd = 2, col = "darkblue",
ylab = "invlogit(x)", main = "Logistic Function")
abline(h = c(0, 1), col = "gray")
-Windows では を指定しては駄目．例えば，setwd(“c:”) とするとエラーが出る．setwd(“c:/usr”) とするか setwd(“c:\usr”) こと．
#setwd("C:/Users/721540/Documents/practice/ロジスティック回帰/~/practice ")
#gesui = read_csv("gesuidou.csv")
#getwd()下水道データ読み込み# 基本統計量表示 gesui # 教科書ではlogit
#setwd("C:/Users/721540/Documents/practice/ロジスティック回帰/~/practice ")
#gesui = read_csv("gesuidou.csv")
#gesui = read_csv("osui.csv")
gesui = read_csv("enbi.csv")## Parsed with column specification:
## cols(
## OBJECTID = col_double(),
## sys_name = col_double(),
## slope = col_double(),
## uedokaburi = col_double(),
## masuhonsuu = col_double(),
## long = col_double(),
## kubun = col_double(),
## did = col_double(),
## kouhou = col_double(),
## nendo = col_double(),
## ekijyouka = col_double(),
## kyouyounensuu = col_double(),
## kansyu = col_double(),
## kei = col_double(),
## kinkyuudo = col_double(),
## taisyo = col_double()
## )gesui<- data.frame(gesui) # 教科書ではlogit
#gesui <- gesui[-1] #OBJECTID列をデータから削除
gesui <- gesui[-1:-2] #OBJECTID,sys_name列をデータから削除
gesui <- gesui[-13]stargazer(as.data.frame(gesui),type = "html")| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| slope | 282 | 3.309 | 2.017 | 0.000 | 1.900 | 4.100 | 9.900 |
| uedokaburi | 282 | 4.218 | 2.570 | 1.009 | 2.462 | 5.397 | 13.385 |
| masuhonsuu | 282 | 1.284 | 1.765 | 0 | 0 | 2 | 11 |
| long | 282 | 31.300 | 15.309 | 0.970 | 21.325 | 40.492 | 96.820 |
| kubun | 282 | 1.209 | 0.407 | 1 | 1 | 1 | 2 |
| did | 282 | 0.766 | 0.424 | 0 | 1 | 1 | 1 |
| kouhou | 282 | 0.337 | 0.473 | 0 | 0 | 1 | 1 |
| nendo | 282 | 1,988.486 | 5.204 | 1,976 | 1,989 | 1,991 | 2,006 |
| ekijyouka | 282 | 0.202 | 0.402 | 0 | 0 | 0 | 1 |
| kyouyounensuu | 282 | 27.514 | 5.204 | 10 | 25 | 27 | 40 |
| kansyu | 282 | 2.000 | 0.000 | 2 | 2 | 2 | 2 |
| kei | 282 | 390.248 | 162.287 | 200 | 250 | 600 | 900 |
| taisyo | 282 | 0.312 | 0.464 | 0 | 0 | 1 | 1 |
ggplot(gesui,aes(slope))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
ggplot(gesui,aes(long))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
ggplot(gesui,aes(uedokaburi))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
ggplot(gesui,aes(sitadokaburi))+ geom_histogram(aes(y = ..density..), bins = 10, colour = “gray”, fill = “blue”)+ theme_classic()
- 桝本数
ggplot(gesui,aes(masuhonsuu ))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
- 施工年度
ggplot(gesui,aes(nendo))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
- 管路径
ggplot(gesui,aes(kei))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
- 液状化
ggplot(gesui,aes(ekijyouka))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
- 緊急度
ggplot(gesui,aes(taisyo))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
- 散布図-plot
head(gesui)## slope uedokaburi masuhonsuu long kubun did kouhou nendo ekijyouka
## 1 1.22 1.054575 1 3.39 2 1 0 2004 1
## 2 2.50 1.533001 0 7.78 2 1 0 1988 0
## 3 4.71 1.414000 0 5.02 2 1 0 1992 1
## 4 1.10 1.544714 3 13.17 2 1 0 1992 1
## 5 1.80 4.412133 1 5.56 2 1 1 1992 1
## 6 8.90 1.738222 0 15.72 2 1 0 1992 1
## kyouyounensuu kansyu kei taisyo
## 1 12 2 200 0
## 2 28 2 250 1
## 3 24 2 250 0
## 4 24 2 250 0
## 5 24 2 250 1
## 6 24 2 250 0plot(gesui[,4],gesui[,1],type="p",xlab="Length of long", ylab="Length of slope", main="延長と傾斜の散布図",cex=0.8,col="red")
#plot(gesui[,4],gesui[,2],type="l",xlab="Length of slope", ylab="Length of long",cex=2,col="red")
#参照頁はcol=”red”でエラー
plot(gesui[,4],gesui[,2],type="p",xlab="Length of long", ylab="Length of slope", main="延長と土被りの散布図",cex=0.8,col="red")
関数 plot(x, y, xlim=range(x), ylim=range(y), type=“p”, main, xlab, ylab……) x 横軸のデータ y 縦軸のデータ xlim 横軸の範囲 ylim 縦軸の範囲 type “p” は、点を描く “l” は、 点を通過する線のみを描く “b” は、点を描き、線で点を連結 “c” は、点を飛ばして線を描く “h” は、各点から横軸までの垂直線を描く “n” は、枠(軸)のみを描く main 散布図の上部にタイトルを付ける xlab x軸(横軸)の名前を付ける ylab y軸(縦軸)の名前を付ける cex 点のサイズ、指定しない場合は１
対散布図（散布図行列） 用いる変数が2以上で、かつ変数がそれほど多くない量的データ場合は、本格的な解析を行う前 に、すべての変数を組み合わせた散布図について考察を行うことで、データ間の関連性を視覚的に把握することができる。1つの画面に複数の変数を組み合わせた散布図を対散布図、あるいは散布図行列という。対散布図の作成は、関数 pairs を用いる。gesui のデータを用いて対散布図の作成について説明する。iris データの四つの変数を用いた最もシンプルな対散布図は次のコマンドで作成できる。
pairs(gesui[1:7], pch = 21, bg = c("red", "green3", "blue", "orenge","black", "yellow", "skyblue")[unclass(gesui$kinkyuudo)])
pairs(gesui[1:13], pch = 21)
- Rのpairs関数の出力をカスタマイズする
- https://entertainment-lab.blogspot.com/2013/01/r_27.html
upper <- function(x, y, …){ # 右上の散布図を設定 oldpar <- par(usr = c(0, 1, 0, 1)) # 3行下でテキストの位置を(0.5, 0.5)にするための調整 v <- abs(cor(x, y)) sz <- 1 + v * 2 # テキストのサイズを直線性の強さで調整 text(0.5, 0.5, sprintf(“%.3f”, v), cex = sz) par(oldpar) }
lower <- function(x, y, …){ # 左下の散布図を設定 points(x, y, pch = 21, bg = c(“red”, “green3”, “blue”)[unclass(iris$Species)]) if( abs(cor(x, y)) > 0.8 ){ # 直線性の強いデータの組だけ回帰式を作成 lmobj <- lm(y~x) abline(lmobj) } }
pairs(gesui[,2:12], upper.panel=upper, lower.panel=lower)
各変数の正規分布性の確認(qqnorm()関数）s悦明変数は正規分布しているか？
https://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1283089632 - https://sites.google.com/site/webtextofr/graph - 正規分布ならばQQプロットは直線になるはずです．
まずSLOPE
qqnorm(gesui$slope)
対数正規性を調べる http://cse.naro.affrc.go.jp/takezawa/r-tips/r/63.html https://tutorialmore.com/questions-1478863.htm
x <- gesui$slope
xy <- qqplot(qt(ppoints(50), df=5), x)
abline(lsfit(xy$x, xy$y), lty=2)
plot(qlnorm(ppoints(x)), sort(x))
次に桝本数：カテゴリー変数であり正規分布にはならない
qqnorm(gesui$masu)
- 上土被り
qqnorm(gesui$uedokaburi)
下土被り ＃ qqnorm(gesui$sitadokaburi)
桝本数
qqnorm(gesui$masu)
- 延長
qqnorm(gesui$long)
シャピロ・ウィルク検定
-敢えて検定するまでもなく今までの分布形状を見ると、左側に偏った形状をしている。 https://tkstock.site/2017/04/12/2017-04-12-000000_1/ で傾斜に対してシャピロ・ウィルク検定を行うと
shapiro.test(gesui$slope)##
## Shapiro-Wilk normality test
##
## data: gesui$slope
## W = 0.86861, p-value = 8.629e-15Shapiro-Wilk normality test
data: gesui$slope
W = 0.89138, p-value < 2.2e-16
結果の W は本検定の検定統計量を示す．この検定では，p値<2.2e-16i以下であるため，有意水準が5%で帰無仮説が棄却され (p ≤ α)，データXの分布は正規分布に従うとはいえない (実質的に従わない) と判断できる．本検定にて計算されるp値はデータ数が4個以上のときは近似的に求められているものである．
同様にして管路延長も正規分布でない
shapiro.test(gesui$long)##
## Shapiro-Wilk normality test
##
## data: gesui$long
## W = 0.97596, p-value = 0.0001102- Shapiro-Wilk normality test
-data: gesui$long - W = 0.92612, p-value = 1.365e-14
同様に桝も正規分でない
shapiro.test(gesui$masu)##
## Shapiro-Wilk normality test
##
## data: gesui$masu
## W = 0.73761, p-value < 2.2e-16Shapiro-Wilk normality test
data: gesui$masu W = 0.73461, p-value < 2.2e-16
# 説明変数slope等は正規分布する必要はない。必要なのは誤差分散が正規分布しているかいなかである
slopeの確立密度関数
https://oku.edu.mie-u.ac.jp/~okumura/stat/first.html x=slope
x <- gesui$slope
curve(dnorm(x), 6, -6)
slope等の非正規なデータの最尤法フィッティング
下準備：http://tips-r.blogspot.com/2014/05/repseq.html
rep(1, 10) #「○を○個だけrepeat」って感じ## [1] 1 1 1 1 1 1 1 1 1 1slope等の非正規なデータの最尤法フィッティング
-https://oku.edu.mie-u.ac.jp/~okumura/stat/141115.html
slopeの平均 -平均3.509，標準偏差2.128の正規分布の密度関数 と指数関数 とを50：10で混ぜ合わせた簡単な関数を考えます：
何をやりたいのか途中でわからなくなったので中止 x=478 y=gesui$slope x= 1:20 y = c(11,4,13,10,4,8,6,16,7,12,10,13,6,5,1,4,2,0,0,1)
z = rep(x, y) r = glm(y ~ dnorm(x,10,3) + exp(-x/10) - 1, family=poisson(link=“identity”)) summary(r) curve(dnorm(x,10,3))
目的変数・説明変数共に正規性は関係ないこととしてロジスティック回帰式の推定値を求める
-http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/15_logit.html
#0305より以下のファクター化について検証する。
gesui$taisyo <- as.factor(gesui$taisyo)
gesui$kansyu <- as.factor(gesui$kansyu)
gesui$kubun <- as.factor(gesui$kubun)
gesui$did <- as.factor(gesui$did)
gesui$ekijyouka <- as.factor(gesui$ekijyouka)
#gesui$kinkyuudo <- as.factor(gesui$kinkyuudo)
sapply(gesui, class)## slope uedokaburi masuhonsuu long kubun
## "numeric" "numeric" "numeric" "numeric" "factor"
## did kouhou nendo ekijyouka kyouyounensuu
## "factor" "numeric" "numeric" "factor" "numeric"
## kansyu kei taisyo
## "factor" "numeric" "factor"dim(gesui)## [1] 282 13summary(gesui)## slope uedokaburi masuhonsuu long kubun
## Min. :0.000 Min. : 1.009 Min. : 0.000 Min. : 0.97 1:223
## 1st Qu.:1.900 1st Qu.: 2.462 1st Qu.: 0.000 1st Qu.:21.32 2: 59
## Median :2.685 Median : 3.402 Median : 1.000 Median :30.06
## Mean :3.309 Mean : 4.218 Mean : 1.284 Mean :31.30
## 3rd Qu.:4.100 3rd Qu.: 5.397 3rd Qu.: 2.000 3rd Qu.:40.49
## Max. :9.900 Max. :13.385 Max. :11.000 Max. :96.82
## did kouhou nendo ekijyouka kyouyounensuu kansyu
## 0: 66 Min. :0.0000 Min. :1976 0:225 Min. :10.00 2:282
## 1:216 1st Qu.:0.0000 1st Qu.:1989 1: 57 1st Qu.:25.00
## Median :0.0000 Median :1991 Median :25.00
## Mean :0.3369 Mean :1988 Mean :27.51
## 3rd Qu.:1.0000 3rd Qu.:1991 3rd Qu.:27.00
## Max. :1.0000 Max. :2006 Max. :40.00
## kei taisyo
## Min. :200.0 0:194
## 1st Qu.:250.0 1: 88
## Median :250.0
## Mean :390.2
## 3rd Qu.:600.0
## Max. :900.0names(gesui)## [1] "slope" "uedokaburi" "masuhonsuu" "long"
## [5] "kubun" "did" "kouhou" "nendo"
## [9] "ekijyouka" "kyouyounensuu" "kansyu" "kei"
## [13] "taisyo"stargazer(as.data.frame(gesui),type = "html")| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| slope | 282 | 3.309 | 2.017 | 0.000 | 1.900 | 4.100 | 9.900 |
| uedokaburi | 282 | 4.218 | 2.570 | 1.009 | 2.462 | 5.397 | 13.385 |
| masuhonsuu | 282 | 1.284 | 1.765 | 0 | 0 | 2 | 11 |
| long | 282 | 31.300 | 15.309 | 0.970 | 21.325 | 40.492 | 96.820 |
| kouhou | 282 | 0.337 | 0.473 | 0 | 0 | 1 | 1 |
| nendo | 282 | 1,988.486 | 5.204 | 1,976 | 1,989 | 1,991 | 2,006 |
| kyouyounensuu | 282 | 27.514 | 5.204 | 10 | 25 | 27 | 40 |
| kei | 282 | 390.248 | 162.287 | 200 | 250 | 600 | 900 |
model_1 <- glm(taisyo ~ slope + uedokaburi + masuhonsuu + long + kubun + did + kouhou + ekijyouka + kyouyounensuu + kei, data = gesui,family = binomial(link = "logit"))
model_1##
## Call: glm(formula = taisyo ~ slope + uedokaburi + masuhonsuu + long +
## kubun + did + kouhou + ekijyouka + kyouyounensuu + kei, family = binomial(link = "logit"),
## data = gesui)
##
## Coefficients:
## (Intercept) slope uedokaburi masuhonsuu long
## -5.6129170 0.1093328 -0.2013296 0.0507254 0.0156753
## kubun2 did1 kouhou ekijyouka1 kyouyounensuu
## -0.9457226 0.3048199 0.3984623 -0.4591354 0.1744175
## kei
## -0.0006761
##
## Degrees of Freedom: 281 Total (i.e. Null); 271 Residual
## Null Deviance: 350.1
## Residual Deviance: 294.5 AIC: 316.5summary(model_1)##
## Call:
## glm(formula = taisyo ~ slope + uedokaburi + masuhonsuu + long +
## kubun + did + kouhou + ekijyouka + kyouyounensuu + kei, family = binomial(link = "logit"),
## data = gesui)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9704 -0.7882 -0.5798 0.8758 2.4113
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.6129170 1.2987147 -4.322 1.55e-05 ***
## slope 0.1093328 0.0703765 1.554 0.1203
## uedokaburi -0.2013296 0.0915007 -2.200 0.0278 *
## masuhonsuu 0.0507254 0.0862614 0.588 0.5565
## long 0.0156753 0.0111230 1.409 0.1588
## kubun2 -0.9457226 0.4558858 -2.074 0.0380 *
## did1 0.3048199 0.4019810 0.758 0.4483
## kouhou 0.3984623 0.5226390 0.762 0.4458
## ekijyouka1 -0.4591354 0.4017207 -1.143 0.2531
## kyouyounensuu 0.1744175 0.0326481 5.342 9.17e-08 ***
## kei -0.0006761 0.0013735 -0.492 0.6226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 350.10 on 281 degrees of freedom
## Residual deviance: 294.49 on 271 degrees of freedom
## AIC: 316.49
##
## Number of Fisher Scoring iterations: 4-上記表の見方は以下のサイトに掲載されている。 https://oku.edu.mie-u.ac.jp/~okumura/stat/logistic.html https://oku.edu.mie-u.ac.jp/~okumura/stat/140921.html 係数は非常に有意なもの（***）からまったく有意でないもの（無印）まで多様である。 -5％で統計的に優位なパラメータは、nendo,kansyu,5％ kei とdidは1％で有意となった。即ち影響があった。
あまり関係なさそうな変数を手作業で除いていってもよいが，自動化するには，step() という関数を使う： これは，デフォルトではフルモデルから1個ずつ変数を外すことを試み，外すと一番AICが良く（小さく）なるものを実際に外す。外すことによってAICが改善しなくなるまで続ける。
y=gesui$taisyo
x=gesui$slope + gesui$uedokaburi + gesui$masuhonsuu + gesui$long + gesui$kubun + gesui$did + gesui$kouhou + gesui$ekijyouka + gesui$kyouyounensuu + gesui$kei ## Warning in Ops.factor(gesui$slope + gesui$uedokaburi + gesui$masuhonsuu + : '+'
## not meaningful for factors## Warning in Ops.factor(gesui$slope + gesui$uedokaburi + gesui$masuhonsuu + : '+'
## not meaningful for factors## Warning in Ops.factor(gesui$slope + gesui$uedokaburi + gesui$masuhonsuu + : '+'
## not meaningful for factorsplot(x, y, xlim=c(0,100),ylim=c(-0.2,1.2), yaxt="n")
axis(2, at=c(0,1), labels=c(0,1), las=1)
curve(1 / (1 + exp(-(0.08319 * x - 3.56668))), add=TRUE)
fitted(model_1)## 1 2 3 4 5 6
## 0.008747289 0.189617095 0.084170525 0.073914071 0.054624257 0.138664030
## 7 8 9 10 11 12
## 0.744858292 0.714723090 0.249135324 0.532773990 0.536150649 0.215311302
## 13 14 15 16 17 18
## 0.494906618 0.104345177 0.529591222 0.553280843 0.220533228 0.095490621
## 19 20 21 22 23 24
## 0.139721161 0.168287094 0.152825444 0.125061957 0.121563863 0.255344535
## 25 26 27 28 29 30
## 0.238946325 0.275273368 0.608818934 0.821614171 0.222291128 0.228683384
## 31 32 33 34 35 36
## 0.325626394 0.038040378 0.179718332 0.611933921 0.391937046 0.218474222
## 37 38 39 40 41 42
## 0.376775496 0.228687903 0.271550546 0.256379263 0.215755218 0.089377430
## 43 44 45 46 47 48
## 0.104327601 0.080403010 0.033682282 0.145489541 0.072316571 0.213318910
## 49 50 51 52 53 54
## 0.221105840 0.235561780 0.073417024 0.217798721 0.253053634 0.373172727
## 55 56 57 58 59 60
## 0.189583217 0.695717932 0.501657793 0.524527940 0.502664060 0.538682781
## 61 62 63 64 65 66
## 0.483944686 0.656554818 0.662540868 0.627207502 0.632437504 0.346992254
## 67 68 69 70 71 72
## 0.632107038 0.676727247 0.610354071 0.577782854 0.140144824 0.158531391
## 73 74 75 76 77 78
## 0.144705311 0.142581013 0.081228520 0.146423894 0.082523053 0.111057135
## 79 80 81 82 83 84
## 0.083000093 0.088704317 0.165344782 0.419623118 0.624289160 0.267284790
## 85 86 87 88 89 90
## 0.262207644 0.160867383 0.307699694 0.205629459 0.193141445 0.233980381
## 91 92 93 94 95 96
## 0.195879617 0.193474660 0.240330255 0.236470059 0.286368438 0.316158552
## 97 98 99 100 101 102
## 0.235911757 0.284777581 0.322506108 0.284985081 0.308084134 0.244404141
## 103 104 105 106 107 108
## 0.224616780 0.372733887 0.621812838 0.126697484 0.514177052 0.360517182
## 109 110 111 112 113 114
## 0.183410426 0.094096563 0.064697010 0.615375893 0.249466211 0.229250030
## 115 116 117 118 119 120
## 0.260825738 0.200998806 0.229004096 0.243966295 0.130535943 0.100301256
## 121 122 123 124 125 126
## 0.060678521 0.079772818 0.069604077 0.071900431 0.102049606 0.117285449
## 127 128 129 130 131 132
## 0.154895264 0.205517532 0.178687266 0.127021777 0.253226357 0.063879560
## 133 134 135 136 137 138
## 0.179656640 0.211319835 0.155127282 0.155126073 0.655417278 0.229234577
## 139 140 141 142 143 144
## 0.250583359 0.245594489 0.247638202 0.283636553 0.237196439 0.305006731
## 145 146 147 148 149 150
## 0.152353300 0.212226561 0.191385158 0.240933531 0.257469154 0.248192510
## 151 152 153 154 155 156
## 0.207712667 0.249731678 0.206776275 0.325251427 0.361129660 0.244674344
## 157 158 159 160 161 162
## 0.316531084 0.286435970 0.293013835 0.278280882 0.023911144 0.239733974
## 163 164 165 166 167 168
## 0.267795612 0.305503046 0.251019369 0.381593325 0.132044548 0.095571628
## 169 170 171 172 173 174
## 0.081312430 0.179464449 0.163124647 0.165627459 0.115482421 0.211675603
## 175 176 177 178 179 180
## 0.168807435 0.218136478 0.227079438 0.221112139 0.340395990 0.852927112
## 181 182 183 184 185 186
## 0.854027121 0.866446456 0.852757988 0.326648654 0.345982187 0.256619686
## 187 188 189 190 191 192
## 0.231776100 0.253967128 0.874677074 0.572561260 0.856482295 0.863587912
## 193 194 195 196 197 198
## 0.862004577 0.858130454 0.870311860 0.876667812 0.607897579 0.627289434
## 199 200 201 202 203 204
## 0.300543115 0.255194516 0.410938456 0.410144433 0.377602758 0.389072842
## 205 206 207 208 209 210
## 0.426793763 0.424918330 0.454160752 0.426678512 0.389830290 0.507267875
## 211 212 213 214 215 216
## 0.831958278 0.877802482 0.878487566 0.601149816 0.417719955 0.475615938
## 217 218 219 220 221 222
## 0.425729269 0.537967414 0.353806261 0.292441840 0.134989368 0.196225566
## 223 224 225 226 227 228
## 0.240081233 0.339263213 0.154546091 0.273145176 0.209704048 0.840282540
## 229 230 231 232 233 234
## 0.831434726 0.470452686 0.382577889 0.372431630 0.429243582 0.327376787
## 235 236 237 238 239 240
## 0.269640759 0.254085926 0.161330633 0.332568503 0.403237990 0.338927937
## 241 242 243 244 245 246
## 0.353738591 0.294675167 0.230834343 0.152133388 0.118316485 0.212365898
## 247 248 249 250 251 252
## 0.203740159 0.140841209 0.088763951 0.077698344 0.091951396 0.087568078
## 253 254 255 256 257 258
## 0.106363967 0.062444115 0.074111007 0.289716731 0.288198996 0.269040841
## 259 260 261 262 263 264
## 0.210010864 0.320802279 0.380977704 0.279845278 0.374249869 0.278665995
## 265 266 267 268 269 270
## 0.301789629 0.294891362 0.303153597 0.230440177 0.254690206 0.226299818
## 271 272 273 274 275 276
## 0.323454938 0.296175946 0.339922998 0.281993243 0.266168883 0.385620457
## 277 278 279 280 281 282
## 0.305320838 0.330433095 0.278436915 0.369152326 0.322467902 0.285397697#predition = predict(model_1, gesui)
#preditionr2 = step(model_1)## Start: AIC=316.49
## taisyo ~ slope + uedokaburi + masuhonsuu + long + kubun + did +
## kouhou + ekijyouka + kyouyounensuu + kei
##
## Df Deviance AIC
## - kei 1 294.73 314.73
## - masuhonsuu 1 294.83 314.83
## - did 1 295.07 315.07
## - kouhou 1 295.07 315.07
## - ekijyouka 1 295.83 315.83
## - long 1 296.48 316.48
## <none> 294.49 316.49
## - slope 1 296.89 316.89
## - kubun 1 299.06 319.06
## - uedokaburi 1 300.20 320.20
## - kyouyounensuu 1 330.52 350.52
##
## Step: AIC=314.73
## taisyo ~ slope + uedokaburi + masuhonsuu + long + kubun + did +
## kouhou + ekijyouka + kyouyounensuu
##
## Df Deviance AIC
## - masuhonsuu 1 295.10 313.10
## - kouhou 1 295.12 313.12
## - did 1 295.71 313.71
## - ekijyouka 1 296.10 314.10
## - long 1 296.58 314.58
## <none> 294.73 314.73
## - slope 1 297.37 315.37
## - kubun 1 299.07 317.07
## - uedokaburi 1 300.44 318.44
## - kyouyounensuu 1 331.65 349.65
##
## Step: AIC=313.1
## taisyo ~ slope + uedokaburi + long + kubun + did + kouhou + ekijyouka +
## kyouyounensuu
##
## Df Deviance AIC
## - kouhou 1 295.37 311.37
## - did 1 296.45 312.45
## - ekijyouka 1 296.58 312.58
## <none> 295.10 313.10
## - long 1 297.54 313.54
## - slope 1 297.71 313.71
## - kubun 1 299.30 315.30
## - uedokaburi 1 301.01 317.01
## - kyouyounensuu 1 331.72 347.72
##
## Step: AIC=311.37
## taisyo ~ slope + uedokaburi + long + kubun + did + ekijyouka +
## kyouyounensuu
##
## Df Deviance AIC
## - did 1 296.54 310.54
## - ekijyouka 1 296.64 310.64
## <none> 295.37 311.37
## - slope 1 297.88 311.88
## - long 1 298.62 312.62
## - kubun 1 299.96 313.96
## - uedokaburi 1 302.53 316.53
## - kyouyounensuu 1 332.12 346.12
##
## Step: AIC=310.54
## taisyo ~ slope + uedokaburi + long + kubun + ekijyouka + kyouyounensuu
##
## Df Deviance AIC
## - ekijyouka 1 297.54 309.54
## <none> 296.54 310.54
## - long 1 299.14 311.14
## - slope 1 299.51 311.51
## - kubun 1 300.55 312.55
## - uedokaburi 1 303.60 315.60
## - kyouyounensuu 1 334.79 346.79
##
## Step: AIC=309.54
## taisyo ~ slope + uedokaburi + long + kubun + kyouyounensuu
##
## Df Deviance AIC
## <none> 297.54 309.54
## - long 1 299.73 309.73
## - slope 1 300.62 310.62
## - kubun 1 303.69 313.69
## - uedokaburi 1 304.34 314.34
## - kyouyounensuu 1 336.30 346.30r3 <- glm(taisyo ~ long + slope + kubun + uedokaburi +kyouyounensuu , data = gesui,family = binomial(link = "logit"))
summary(r3)##
## Call:
## glm(formula = taisyo ~ long + slope + kubun + uedokaburi + kyouyounensuu,
## family = binomial(link = "logit"), data = gesui)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9399 -0.7789 -0.6057 0.9437 2.4390
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.68164 1.05957 -5.362 8.22e-08 ***
## long 0.01444 0.00980 1.474 0.1406
## slope 0.12110 0.06872 1.762 0.0780 .
## kubun2 -0.99217 0.41994 -2.363 0.0181 *
## uedokaburi -0.17575 0.07165 -2.453 0.0142 *
## kyouyounensuu 0.17621 0.03183 5.537 3.08e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 350.10 on 281 degrees of freedom
## Residual deviance: 297.54 on 276 degrees of freedom
## AIC: 309.54
##
## Number of Fisher Scoring iterations: 4r4 = step(r3)## Start: AIC=309.54
## taisyo ~ long + slope + kubun + uedokaburi + kyouyounensuu
##
## Df Deviance AIC
## <none> 297.54 309.54
## - long 1 299.73 309.73
## - slope 1 300.62 310.62
## - kubun 1 303.69 313.69
## - uedokaburi 1 304.34 314.34
## - kyouyounensuu 1 336.30 346.30round(fitted(r3))## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 0 0 0 0 0 0 1 1 0 1 1 0 0 0 1 1 0 0 0 0
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 1 1 1 1 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 1 1 1 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
## 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
## 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 281 282
## 0 0上記の結果、long + slope + kubun + uedokaburi +kyouyounensuu 以上に変数を減らしても一番AICが良く（小さく）なるものなく予測精度は向上しないことが判ったためこのモデルで予測行う
このモデルで予測を行うと
plot(fitted(r3), gesui$taisyo, ylim=c(-0.2,1.2), yaxt="n")
axis(2, c(0,1))
predict(r3, type = "response")## 1 2 3 4 5 6 7
## 0.01047846 0.16879968 0.11401981 0.08372471 0.05108327 0.19071246 0.74906640
## 8 9 10 11 12 13 14
## 0.68063511 0.27796436 0.51119454 0.51481085 0.19379670 0.45301715 0.09245345
## 15 16 17 18 19 20 21
## 0.57207984 0.63884529 0.18262554 0.10916858 0.15875300 0.24449432 0.21849889
## 22 23 24 25 26 27 28
## 0.18562137 0.17462883 0.36113327 0.37720948 0.25508842 0.58387732 0.80970806
## 29 30 31 32 33 34 35
## 0.21481036 0.18859096 0.28354515 0.03994467 0.18303967 0.60676926 0.39955801
## 36 37 38 39 40 41 42
## 0.20538023 0.36463775 0.22070669 0.25024810 0.24632609 0.20149183 0.07782016
## 43 44 45 46 47 48 49
## 0.09166838 0.07268942 0.03494991 0.13446824 0.07445423 0.20549213 0.20754535
## 50 51 52 53 54 55 56
## 0.23335847 0.06369104 0.22308868 0.25566349 0.29838230 0.19765038 0.74167183
## 57 58 59 60 61 62 63
## 0.56467723 0.57595189 0.56727367 0.60887350 0.53895526 0.61679878 0.62073619
## 64 65 66 67 68 69 70
## 0.55567832 0.58898039 0.37053157 0.60373737 0.62721636 0.65490105 0.65177254
## 71 72 73 74 75 76 77
## 0.11682593 0.21293549 0.19019545 0.18827376 0.10363686 0.19399950 0.10497207
## 78 79 80 81 82 83 84
## 0.12811778 0.10680785 0.10307365 0.21559270 0.44854580 0.51946341 0.28010578
## 85 86 87 88 89 90 91
## 0.27463383 0.19558906 0.33094610 0.25690841 0.23949363 0.21166381 0.16395516
## 92 93 94 95 96 97 98
## 0.16928500 0.36760039 0.31038637 0.29833999 0.34130224 0.32808508 0.38443820
## 99 100 101 102 103 104 105
## 0.40824345 0.38175007 0.42102845 0.26396380 0.24611706 0.35011620 0.61159972
## 106 107 108 109 110 111 112
## 0.15074207 0.38993641 0.30972225 0.15796125 0.06917479 0.05703547 0.53896247
## 113 114 115 116 117 118 119
## 0.14482871 0.13423114 0.12850508 0.16926398 0.18379403 0.22668980 0.10547853
## 120 121 122 123 124 125 126
## 0.08933132 0.05910159 0.07632067 0.06800294 0.06989630 0.09491246 0.10821304
## 127 128 129 130 131 132 133
## 0.13844257 0.18458510 0.22075542 0.16090577 0.28552508 0.06056554 0.21693757
## 134 135 136 137 138 139 140
## 0.26088322 0.18951499 0.19252436 0.66720469 0.24028542 0.26090958 0.25544388
## 141 142 143 144 145 146 147
## 0.25857343 0.29758306 0.24936915 0.31995163 0.16982674 0.22833630 0.20622749
## 148 149 150 151 152 153 154
## 0.26110605 0.27124021 0.26043966 0.21741948 0.26251313 0.21971837 0.40213293
## 155 156 157 158 159 160 161
## 0.41794544 0.31686772 0.42854003 0.38663479 0.37596648 0.35175942 0.03513619
## 162 163 164 165 166 167 168
## 0.33120709 0.35404968 0.39594788 0.34387471 0.42818965 0.11941883 0.09138043
## 169 170 171 172 173 174 175
## 0.09978569 0.15554224 0.14468338 0.17780560 0.17370011 0.27710503 0.24894879
## 176 177 178 179 180 181 182
## 0.23686607 0.23990679 0.23494549 0.44006899 0.82927930 0.83604872 0.84851029
## 183 184 185 186 187 188 189
## 0.84104979 0.30274850 0.30741005 0.30481720 0.28081724 0.30211356 0.85238251
## 190 191 192 193 194 195 196
## 0.52654058 0.83361176 0.84046639 0.83804841 0.83297113 0.84556643 0.85552636
## 197 198 199 200 201 202 203
## 0.57797475 0.60565917 0.34056741 0.30615692 0.36407987 0.36434306 0.34552714
## 204 205 206 207 208 209 210
## 0.35591810 0.35693669 0.36163698 0.40546227 0.37790110 0.34020963 0.46606426
## 211 212 213 214 215 216 217
## 0.81865998 0.83703191 0.84351787 0.57136959 0.34325686 0.42167812 0.38336344
## 218 219 220 221 222 223 224
## 0.48316961 0.25288383 0.24222144 0.12930172 0.20423849 0.19235036 0.27579963
## 225 226 227 228 229 230 231
## 0.11847738 0.22110575 0.17025689 0.84766280 0.84027180 0.49233486 0.41153129
## 232 233 234 235 236 237 238
## 0.40192213 0.34537627 0.29039970 0.21399922 0.22840300 0.13874704 0.25632263
## 239 240 241 242 243 244 245
## 0.39570333 0.25002291 0.29938411 0.29601484 0.23152607 0.15640746 0.16137823
## 246 247 248 249 250 251 252
## 0.24333245 0.23508524 0.16591162 0.08885907 0.07645103 0.08909981 0.07930001
## 253 254 255 256 257 258 259
## 0.09364770 0.08958573 0.10009410 0.32323505 0.26357040 0.23485055 0.20849195
## 260 261 262 263 264 265 266
## 0.28650043 0.31965635 0.23024594 0.31003537 0.24075995 0.26056884 0.25506425
## 267 268 269 270 271 272 273
## 0.26242362 0.18719726 0.20591597 0.25320326 0.30952665 0.24618817 0.27036681
## 274 275 276 277 278 279 280
## 0.26346087 0.24967733 0.31582744 0.28818514 0.28197780 0.24147386 0.32968203
## 281 282
## 0.26183732 0.23376474round(predict(r3, type = "response")) #四捨五入し、1になるか0かをわかり易く表示## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 0 0 0 0 0 0 1 1 0 1 1 0 0 0 1 1 0 0 0 0
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 1 1 1 1 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 1 1 1 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
## 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
## 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
## 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 281 282
## 0 0-結果は全て0であり、taisyo(劣化がある管路)は一つもないという結果になった。？？？
回帰係数を個別に検定する。
stargazer(r3,
style = "ajps", # AJPS仕様
dep.var.labels = "影響あり or 影響なし", # 応答変数の名前
title = "厚木のテストデータを使ったパラメータの影響度判定", # タイトル
type ="html")| 影響あり or 影響なし | |
| long | 0.014 |
| (0.010) | |
| slope | 0.121* |
| (0.069) | |
| kubun2 | -0.992** |
| (0.420) | |
| uedokaburi | -0.176** |
| (0.072) | |
| kyouyounensuu | 0.176*** |
| (0.032) | |
| Constant | -5.682*** |
| (1.060) | |
| N | 282 |
| Log Likelihood | -148.772 |
| AIC | 309.543 |
| p < .01; p < .05; p < .1 | |
-5％で統計的に優位なパラメータは、nendo,kansyu, kei とdidは1％で有意となった。即ち影響があった
http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/15_logit.html