gesuidouimportance

https://rpubs.com/fumi/569286 # 参考資料 ★★★★　ロジスティック回帰分析　浅野 正彦 http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/15_logit.html https://www.asanoucla.com/

★★★ 1.回帰分析の説明変数や目的変数は正規分布していなくてもよいか？ https://toukeier.hatenablog.com/entry/2019/09/08/224346

回帰分析の説明変数や目的変数は正規分布していなくてもよいか？という質問に対しては、「していなくてよい」が答えだ。 正規分布している必要があるのは残差なのだ。説明変数や目的変数ではない。これを間違えてはいけない。 残差が正規分布している必要がある理由は、分散分析のF検定の検定統計量の分母に該当するから。F分布は、分母・分子とも正規分布していることが前提だからだ。 https://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1283089632

2.統計ブログランキング https://blog.with2.net/rank4465-0.html

3.三重大学奥村先生サイト https://oku.edu.mie-u.ac.jp/~okumura/stat/first.html

4.rの逆引き辞書最新版 https://books.google.co.jp/books?id=YeaoDwAAQBAJ&pg=PA22&dq=r&hl=ja&sa=X&ved=0ahUKEwilhYWLofHnAhUm-2EKHRCKCYMQ6AEIYDAG#v=onepage&q=r&f=false

5.Rでロジスティック回帰分析 そのまま使える自作関数 https://to-kei.net/r-beginner/r-5-logistic/

4.自分のサイト - http://rpubs.com/fumi/569286 # 始めに文法の基本： - “文頭に#+スペースで大みだしができる。” - “以降を箇条書きで書いていく場合は、ここに書いたとおりスペース+”-“+スペースで書ける” - このローカルデータは、“C:\721540.—Rmd”であり - webには、http://rpubs.com/fumi/562118 - バックアップとして、上記のデータのダウンロードを C:\721540

#　厚木の緊急度推定における機械学習法と線形回帰分析の比較検討 - 厚木データのインポート

library(car)

## Loading required package: carData

library(caret)

## Warning: package 'caret' was built under R version 3.6.2

## Loading required package: lattice

library(cluster)
library(dummies)

## dummies-1.5.6 provided by Decision Patterns

library(data.table)

## Warning: package 'data.table' was built under R version 3.6.2

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:data.table':
## 
##     between, first, last

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(e1071)

## Warning: package 'e1071' was built under R version 3.6.2

library(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.2

library(randomForest)

## 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':
## 
##     margin

library(ranger)

## Warning: package 'ranger' was built under R version 3.6.2

## 
## Attaching package: 'ranger'

## The following object is masked from 'package:randomForest':
## 
##     importance

library(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':
## 
##     importance

library(readr)

## Warning: package 'readr' was built under R version 3.6.2

library(rpart.plot)

## Warning: package 'rpart.plot' was built under R version 3.6.2

## Loading required package: rpart

library(rpart)
library(readr)
library(reshape)

## Warning: package 'reshape' was built under R version 3.6.2

## 
## Attaching package: 'reshape'

## The following object is masked from 'package:dplyr':
## 
##     rename

## The following object is masked from 'package:data.table':
## 
##     melt

library(rsconnect)

## Warning: package 'rsconnect' was built under R version 3.6.2

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following objects are masked from 'package:reshape':
## 
##     colsplit, melt, recast

## The following objects are masked from 'package:data.table':
## 
##     dcast, melt

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:reshape2':
## 
##     smiths

## The following objects are masked from 'package:reshape':
## 
##     expand, smiths

library(xtable)
library(nnet)

## Warning: package 'nnet' was built under R version 3.6.2

library(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

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")

## Parsed with column specification:
## cols(
##   OBJECTID = col_double(),
##   slope = col_double(),
##   long = col_double(),
##   uedokaburi = col_double(),
##   sitadokaburi = col_double(),
##   masuhonsuu = col_double(),
##   nendo = col_double(),
##   kei = col_double(),
##   kubun = col_double(),
##   did = col_double(),
##   kouhou = col_double(),
##   ekijyouka = col_double(),
##   kansyu = col_double(),
##   kinkyuudo = col_double(),
##   taisyo = col_double()
## )

gesui<- data.frame(gesui) # 教科書ではlogit
gesui <- gesui[-1] #OBJECTID列をデータから削除
dim(gesui)

## [1] 478  14

summary(gesui)

##      slope            long         uedokaburi      sitadokaburi      
##  Min.   :0.000   Min.   :1.004   Min.   : 1.020   Min.   : 0.000685  
##  1st Qu.:2.000   1st Qu.:2.152   1st Qu.: 2.675   1st Qu.: 1.844750  
##  Median :2.800   Median :3.183   Median : 3.722   Median : 2.693033  
##  Mean   :3.509   Mean   :3.646   Mean   : 4.138   Mean   : 2.983303  
##  3rd Qu.:4.700   3rd Qu.:4.717   3rd Qu.: 4.900   3rd Qu.: 3.820872  
##  Max.   :9.710   Max.   :9.940   Max.   :13.863   Max.   :11.196442  
##    masuhonsuu         nendo           kei             kubun       
##  Min.   : 0.000   Min.   :1976   Min.   :0.2000   Min.   :0.0000  
##  1st Qu.: 0.000   1st Qu.:1978   1st Qu.:0.2500   1st Qu.:1.0000  
##  Median : 1.000   Median :1986   Median :0.6000   Median :1.0000  
##  Mean   : 1.554   Mean   :1985   Mean   :0.5779   Mean   :0.8285  
##  3rd Qu.: 2.000   3rd Qu.:1991   3rd Qu.:0.8000   3rd Qu.:1.0000  
##  Max.   :14.000   Max.   :1992   Max.   :1.6500   Max.   :1.0000  
##       did             kouhou         ekijyouka          kansyu     
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.000  
##  1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.000  
##  Median :1.0000   Median :1.0000   Median :0.0000   Median :1.000  
##  Mean   :0.8264   Mean   :0.5021   Mean   :0.3033   Mean   :0.705  
##  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.000  
##    kinkyuudo        taisyo       
##  Min.   :0.00   Min.   :0.00000  
##  1st Qu.:0.00   1st Qu.:0.00000  
##  Median :2.00   Median :0.00000  
##  Mean   :1.51   Mean   :0.05858  
##  3rd Qu.:3.00   3rd Qu.:0.00000  
##  Max.   :3.00   Max.   :1.00000

names(gesui)

##  [1] "slope"        "long"         "uedokaburi"   "sitadokaburi" "masuhonsuu"  
##  [6] "nendo"        "kei"          "kubun"        "did"          "kouhou"      
## [11] "ekijyouka"    "kansyu"       "kinkyuudo"    "taisyo"

stargazer(as.data.frame(gesui),type = "html")

Statistic	N	Mean	St. Dev.	Min	Pctl(25)	Pctl(75)	Max
slope	478	3.509	2.128	0	2	4.7	10
long	478	3.646	1.844	1.004	2.152	4.717	9.940
uedokaburi	478	4.138	2.096	1.020	2.675	4.900	13.863
sitadokaburi	478	2.983	1.686	0.001	1.845	3.821	11.196
masuhonsuu	478	1.554	2.090	0	0	2	14
nendo	478	1,985.094	6.015	1,976	1,978	1,991	1,992
kei	478	0.578	0.353	0.200	0.250	0.800	1.650
kubun	478	0.828	0.377	0	1	1	1
did	478	0.826	0.379	0	1	1	1
kouhou	478	0.502	0.501	0	0	1	1
ekijyouka	478	0.303	0.460	0	0	1	1
kansyu	478	0.705	0.457	0	0	1	1
kinkyuudo	478	1.510	1.462	0	0	3	3
taisyo	478	0.059	0.235	0	0	0	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(kinkyuudo))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
  theme_classic()

• 散布図-plot

head(gesui)

##   slope  long uedokaburi sitadokaburi masuhonsuu nendo  kei kubun did kouhou
## 1  0.00 9.940   3.758000     3.258000          0  1983 1.10     1   1      0
## 2  5.00 1.232   2.331392     2.255000          0  1982 1.10     1   1      1
## 3  5.00 6.500   4.654844     4.513411          0  1981 0.60     1   1      0
## 4  2.70 5.606   4.584913     3.037999          0  1976 0.60     1   1      0
## 5  4.71 5.020   1.414000     1.403298          0  1992 0.25     0   1      1
## 6  1.10 1.317   1.544714     1.153690          3  1992 0.25     0   1      1
##   ekijyouka kansyu kinkyuudo taisyo
## 1         1      1         0      0
## 2         1      1         3      0
## 3         1      1         3      0
## 4         1      1         3      0
## 5         1      0         0      0
## 6         1      0         0      0

plot(gesui[,2],gesui[,3])

• 散布図-plot 2 https://www1.doshisha.ac.jp/~mjin/R/Chap_07/07.html参照

#plot(gesui[,2],gesui[,3],type="l",xlab="Length of slope", ylab="Length of long",cex=2,col="red")
#参照頁はcol=”red”でエラー
plot(gesui[,2],gesui[,3],type="p",xlab="Length of slope", ylab="Length of long", 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[2:7], pch = 21, bg = c("red", "green3", "blue", "orenge","black", "yellow", "skyblue")[unclass(gesui$kinkyuudo)])

pairs(gesui[2:14], 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:14], 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.89138, p-value < 2.2e-16

• Shapiro-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.92612, p-value = 1.365e-14

• 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.73461, p-value < 2.2e-16

Shapiro-Wilk normality test

data: gesui$masu W = 0.73461, p-value < 2.2e-16

#　説明変数slope等は正規分布する必要はない。必要なのは誤差分散が正規分布しているかいなかである

-https://toukeier.hatenablog.com/entry/2019/09/08/224346

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 1

slope等の非正規なデータの最尤法フィッティング

-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

model_1 <- glm(taisyo ~ slope + long + uedokaburi + sitadokaburi + masuhonsuu + nendo + kei + kubun + did + kouhou + ekijyouka + kansyu, data = gesui,family = binomial(link = "logit"))
summary(model_1)

## 
## Call:
## glm(formula = taisyo ~ slope + long + uedokaburi + sitadokaburi + 
##     masuhonsuu + nendo + kei + kubun + did + kouhou + ekijyouka + 
##     kansyu, family = binomial(link = "logit"), data = gesui)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0921  -0.3082  -0.1684  -0.1077   3.9648  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -487.54742  149.09610  -3.270  0.00108 **
## slope          -0.21202    0.14667  -1.446  0.14830   
## long           -0.05377    0.15201  -0.354  0.72357   
## uedokaburi      0.28273    0.22129   1.278  0.20137   
## sitadokaburi   -0.47435    0.32603  -1.455  0.14569   
## masuhonsuu      0.12920    0.09966   1.296  0.19484   
## nendo           0.24178    0.07482   3.231  0.00123 **
## kei             2.37393    1.18788   1.998  0.04567 * 
## kubun          -1.65915    1.02478  -1.619  0.10544   
## did             2.79188    1.37822   2.026  0.04279 * 
## kouhou          0.25207    0.83722   0.301  0.76335   
## ekijyouka      -0.51234    0.70686  -0.725  0.46857   
## kansyu          3.32820    1.12433   2.960  0.00307 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 213.22  on 477  degrees of freedom
## Residual deviance: 162.56  on 465  degrees of freedom
## AIC: 188.56
## 
## Number of Fisher Scoring iterations: 7

-5％で統計的に優位なパラメータは、nendo,kansyu,5％　kei とdidは1％で有意となった。即ち影響があった。 次にこのモデルで予測を行うと

predict(model_1, type = "response")

##            1            2            3            4            5            6 
## 0.0725040033 0.0448721395 0.0043359577 0.0043658239 0.0130653901 0.0564017744 
##            7            8            9           10           11           12 
## 0.0198760232 0.0197128511 0.0065984796 0.0056155123 0.0054487253 0.0084206018 
##           13           14           15           16           17           18 
## 0.0143467270 0.0478749683 0.1678590705 0.0225336689 0.0208140027 0.0468409812 
##           19           20           21           22           23           24 
## 0.0080280480 0.0158301977 0.0001825539 0.0004996235 0.0067460629 0.0076642463 
##           25           26           27           28           29           30 
## 0.0064056898 0.0088376293 0.0092112471 0.0162465375 0.0123458863 0.0759057171 
##           31           32           33           34           35           36 
## 0.0205716535 0.0594941187 0.0060755214 0.0254147495 0.0811381107 0.0685212456 
##           37           38           39           40           41           42 
## 0.0662464858 0.0256417461 0.0904094057 0.0406520082 0.0101422898 0.0590159071 
##           43           44           45           46           47           48 
## 0.0139437091 0.0076023350 0.0087942592 0.0107405036 0.0008981441 0.0006837599 
##           49           50           51           52           53           54 
## 0.0009911884 0.0017749492 0.0020705528 0.0026745912 0.0202732585 0.0096139654 
##           55           56           57           58           59           60 
## 0.0144259534 0.0082155720 0.0090044035 0.2387274029 0.2176800984 0.2399581367 
##           61           62           63           64           65           66 
## 0.2086613044 0.2619393540 0.1908378323 0.2406194286 0.2277373197 0.1910085787 
##           67           68           69           70           71           72 
## 0.2378551530 0.2145186525 0.2193836251 0.1948299764 0.2313721594 0.0042542114 
##           73           74           75           76           77           78 
## 0.1254535992 0.0108279575 0.0188440829 0.0870726801 0.1549223403 0.0062184516 
##           79           80           81           82           83           84 
## 0.0066360033 0.0048286502 0.0127358289 0.0077619329 0.0100315213 0.0080884635 
##           85           86           87           88           89           90 
## 0.0864833754 0.2867198764 0.2777075572 0.1862056803 0.2486370847 0.0432980664 
##           91           92           93           94           95           96 
## 0.1009961108 0.0613813573 0.2445284817 0.3235728639 0.2878900289 0.0831270544 
##           97           98           99          100          101          102 
## 0.3253913482 0.2480980612 0.2165351043 0.4223100152 0.0824799783 0.1703627928 
##          103          104          105          106          107          108 
## 0.2529973652 0.1488999480 0.2756728926 0.0601561487 0.2399946782 0.0759691641 
##          109          110          111          112          113          114 
## 0.0045572807 0.0048734246 0.0252358351 0.0120204071 0.0575735593 0.0406485690 
##          115          116          117          118          119          120 
## 0.0536589515 0.0059906654 0.0098145236 0.1031670754 0.1296090984 0.1183500780 
##          121          122          123          124          125          126 
## 0.0690381988 0.0121073276 0.0084329711 0.0068277360 0.0181762670 0.0649880355 
##          127          128          129          130          131          132 
## 0.1480533068 0.0740524188 0.1100594829 0.1426727507 0.1426886747 0.0866141651 
##          133          134          135          136          137          138 
## 0.0826255902 0.0687580866 0.1048154654 0.1230285505 0.0996387940 0.1329356066 
##          139          140          141          142          143          144 
## 0.0259353042 0.0264191366 0.0956589931 0.1988794572 0.0393625100 0.0487149985 
##          145          146          147          148          149          150 
## 0.2108789630 0.1333228686 0.0760057483 0.0923711475 0.0353665683 0.0074607034 
##          151          152          153          154          155          156 
## 0.0085853870 0.0078553146 0.0128380299 0.0028488474 0.0320860981 0.1252757781 
##          157          158          159          160          161          162 
## 0.1070639226 0.1015773908 0.1712923182 0.0059077965 0.0010494281 0.0080403704 
##          163          164          165          166          167          168 
## 0.0355546160 0.0512487555 0.0395152336 0.0508355507 0.0064204582 0.0267268913 
##          169          170          171          172          173          174 
## 0.0115934594 0.1561917278 0.0099442882 0.0041572188 0.0095734813 0.0143429846 
##          175          176          177          178          179          180 
## 0.0144698727 0.0096954342 0.0118498594 0.0375668817 0.0405393382 0.0219199620 
##          181          182          183          184          185          186 
## 0.0201527573 0.0188125897 0.0071695145 0.0139421307 0.0200428422 0.0118629136 
##          187          188          189          190          191          192 
## 0.0190098200 0.0264565573 0.0384937500 0.1430715927 0.1591828350 0.0114030642 
##          193          194          195          196          197          198 
## 0.0097386343 0.0044666464 0.0084423941 0.0029784795 0.0060360527 0.0082936808 
##          199          200          201          202          203          204 
## 0.0085772823 0.0082366694 0.0194814447 0.0401224721 0.0189724059 0.0214461220 
##          205          206          207          208          209          210 
## 0.0514377642 0.0466354495 0.0302665281 0.0528450019 0.0188630114 0.0255588074 
##          211          212          213          214          215          216 
## 0.0481057472 0.0059699287 0.0234600340 0.0417339728 0.2814484996 0.0014024426 
##          217          218          219          220          221          222 
## 0.0016169707 0.0065803736 0.0130435736 0.0077240425 0.0269447637 0.0523778487 
##          223          224          225          226          227          228 
## 0.0178763798 0.0168355533 0.0188920365 0.0228099015 0.0350618544 0.0677421067 
##          229          230          231          232          233          234 
## 0.0767755379 0.0292888552 0.0474828727 0.0193091782 0.0234620064 0.0263130132 
##          235          236          237          238          239          240 
## 0.0366193333 0.0044323639 0.0040263671 0.0035128966 0.0068011451 0.0118366056 
##          241          242          243          244          245          246 
## 0.0087150848 0.0141883302 0.0131342345 0.0165399908 0.0094172446 0.0216135531 
##          247          248          249          250          251          252 
## 0.0261084202 0.0395784870 0.0242008747 0.0162974769 0.0204533699 0.0240486939 
##          253          254          255          256          257          258 
## 0.0225613318 0.0288900530 0.1065444974 0.0160576801 0.0121704170 0.0147266052 
##          259          260          261          262          263          264 
## 0.0087297806 0.0144408675 0.0227424852 0.0281021664 0.0227872180 0.0238113365 
##          265          266          267          268          269          270 
## 0.0330708761 0.0268897692 0.0392351606 0.0378249996 0.0497442402 0.0456391347 
##          271          272          273          274          275          276 
## 0.0288289184 0.0327382893 0.0260656676 0.0147353040 0.0132537050 0.0280756812 
##          277          278          279          280          281          282 
## 0.0930888909 0.0677233105 0.0535207464 0.0380424990 0.0350487844 0.1108822276 
##          283          284          285          286          287          288 
## 0.0150880903 0.0174576079 0.0647802033 0.1019953611 0.0630478427 0.0261417532 
##          289          290          291          292          293          294 
## 0.0007070790 0.0006541019 0.0182165443 0.0051960757 0.0357289946 0.0076962713 
##          295          296          297          298          299          300 
## 0.0038627135 0.0058964253 0.0126039070 0.0155120738 0.0029603869 0.0004693318 
##          301          302          303          304          305          306 
## 0.0072338219 0.0106449588 0.0067767190 0.0003859820 0.0002566383 0.0004331259 
##          307          308          309          310          311          312 
## 0.0002694527 0.0002018924 0.0074665085 0.0096965495 0.0003303121 0.0003518156 
##          313          314          315          316          317          318 
## 0.0001919613 0.0057404719 0.0124933981 0.0089949108 0.0096876559 0.0099867038 
##          319          320          321          322          323          324 
## 0.0101024741 0.0106626120 0.0110233934 0.0108392819 0.0111731442 0.0109244550 
##          325          326          327          328          329          330 
## 0.0080019150 0.0023988757 0.0124355588 0.0045263245 0.0024838625 0.0024269251 
##          331          332          333          334          335          336 
## 0.0037844377 0.0030008229 0.0028314493 0.0066004009 0.0048578909 0.0088360030 
##          337          338          339          340          341          342 
## 0.0033892364 0.0060898713 0.0075970032 0.0035238533 0.0099132484 0.0041163738 
##          343          344          345          346          347          348 
## 0.0066471869 0.0087983480 0.0118002251 0.0139560795 0.0128485720 0.0150826394 
##          349          350          351          352          353          354 
## 0.0148116991 0.0143985678 0.0152712715 0.0154876031 0.0112017658 0.0005040848 
##          355          356          357          358          359          360 
## 0.0004143871 0.0004363385 0.0004556623 0.0003827507 0.0005882375 0.0003385529 
##          361          362          363          364          365          366 
## 0.0005365494 0.0004995020 0.0006108586 0.0003635722 0.0003628278 0.0004142098 
##          367          368          369          370          371          372 
## 0.0005997298 0.0004122594 0.0006081718 0.0238978724 0.0346017753 0.0758008507 
##          373          374          375          376          377          378 
## 0.0094261340 0.0108249972 0.0003922537 0.0037270305 0.0086156663 0.0008182783 
##          379          380          381          382          383          384 
## 0.0015595264 0.0016119119 0.0010085383 0.0003845786 0.0004403294 0.0005447982 
##          385          386          387          388          389          390 
## 0.0527375996 0.0363100706 0.0188992282 0.0118427120 0.0777915611 0.0541306907 
##          391          392          393          394          395          396 
## 0.0219089149 0.0162650560 0.0834911937 0.0343928965 0.0252144178 0.0379676109 
##          397          398          399          400          401          402 
## 0.0340900826 0.0326463795 0.0398526961 0.0152039635 0.0018760312 0.0015179910 
##          403          404          405          406          407          408 
## 0.0060611827 0.0113828755 0.0181360540 0.0377117012 0.0149664598 0.0136507462 
##          409          410          411          412          413          414 
## 0.2138529896 0.2017812669 0.2245981159 0.2829840230 0.2029949162 0.2572739823 
##          415          416          417          418          419          420 
## 0.2201157484 0.3128836546 0.1964682294 0.2676889485 0.1776080089 0.2982749504 
##          421          422          423          424          425          426 
## 0.2431801667 0.2130357567 0.2757494497 0.2137064539 0.2651818765 0.2181662382 
##          427          428          429          430          431          432 
## 0.2519675901 0.2583012289 0.1076541317 0.1402065791 0.2826394452 0.3810175935 
##          433          434          435          436          437          438 
## 0.4377475894 0.2403659805 0.3430171535 0.3236465042 0.3769300343 0.4492017843 
##          439          440          441          442          443          444 
## 0.4176103920 0.0156141041 0.0162957238 0.0112647831 0.2153376565 0.0119578190 
##          445          446          447          448          449          450 
## 0.0134157681 0.0004188285 0.0005780872 0.0002755386 0.0011469468 0.0218676916 
##          451          452          453          454          455          456 
## 0.1494798809 0.1533309152 0.1104114466 0.0081089730 0.0092154705 0.0135578085 
##          457          458          459          460          461          462 
## 0.0176130650 0.0036123043 0.0030766300 0.0095421608 0.0081161226 0.0082444649 
##          463          464          465          466          467          468 
## 0.0079809811 0.0030202710 0.0043399718 0.0004520883 0.0096981917 0.0109732813 
##          469          470          471          472          473          474 
## 0.0043349919 0.0055452555 0.0068387213 0.0062069761 0.0038810816 0.0067169680 
##          475          476          477          478 
## 0.0379043157 0.0099554614 0.0036894795 0.0035481054

-結果は最大で0.47であり,全て0.5以下になっている。

回帰係数を個別に検定する。

stargazer(model_1, 
          style = "ajps", # AJPS仕様
          dep.var.labels = "影響あり or 影響なし", # 応答変数の名前
          title = "厚木のテストデータを使ったパラメータの影響度判定", # タイトル
          type ="html")

厚木のテストデータを使ったパラメータの影響度判定

	影響あり or 影響なし
slope	-0.212
	(0.147)
long	-0.054
	(0.152)
uedokaburi	0.283
	(0.221)
sitadokaburi	-0.474
	(0.326)
masuhonsuu	0.129
	(0.100)
nendo	0.242***
	(0.075)
kei	2.374**
	(1.188)
kubun	-1.659
	(1.025)
did	2.792**
	(1.378)
kouhou	0.252
	(0.837)
ekijyouka	-0.512
	(0.707)
kansyu	3.328***
	(1.124)
Constant	-487.547***
	(149.096)
N	478
Log Likelihood	-81.280
AIC	188.559
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

それぞれについて、調べる対象となる値のベクトルを作る。 必要な値は次のとおり。

最小値 最大値 0 1 平均値 - 1/2 平均値 + 1/2 平均値 - sd/2 平均値 + sd/2 sd(x)：不偏標準偏差

slope.mean <- mean(gesui$slope)  # slopeの平均値
slope.sd <- sd(gesui$slope)      # slopeのsd
# slopeの8つの値
slope.vec <- c(min(gesui$slope), max(gesui$slope), 0, 1,
             slope.mean - 1/2, slope.mean + 1/2,
              slope.mean - slope.sd/2, slope.mean + slope.sd/2)


long.mean <- mean(gesui$long)  # longの平均値
long.sd <- sd(gesui$long)      # longのsd
# longの8つの値
long.vec <- c(min(gesui$long), max(gesui$long), 0, 1,
             long.mean - 1/2, long.mean + 1/2,
              long.mean - long.sd/2, long.mean + long.sd/2)

uedokaburi.mean <- mean(gesui$uedokaburi)  # uedokaburiの平均値
uedokaburi.sd <- sd(gesui$uedokaburi)      # uedokaburiのsd
# uedokaburiの8つの値
uedokaburi.vec <- c(min(gesui$uedokaburi), max(gesui$uedokaburi), 0, 1,
             uedokaburi.mean - 1/2, uedokaburi.mean + 1/2,
              uedokaburi.mean - uedokaburi.sd/2, uedokaburi.mean + uedokaburi.sd/2)

sitadokaburi.mean <- mean(gesui$sitadokaburi)  # sitadokaburiの平均値
sitadokaburi.sd <- sd(gesui$sitadokaburi)      # sitadokaburiのsd
# sitadokaburiの8つの値
sitadokaburi.vec <- c(min(gesui$sitadokaburi), max(gesui$sitadokaburi), 0, 1,
             sitadokaburi.mean - 1/2, sitadokaburi.mean + 1/2,
              sitadokaburi.mean - sitadokaburi.sd/2, sitadokaburi.mean + sitadokaburi.sd/2)

masuhonsuu.mean <- mean(gesui$masuhonsuu)  # masuhonsuuの平均値
masuhonsuu.sd <- sd(gesui$masuhonsuu)      # masuhonsuuのsd
# masuhonsuuの8つの値
masuhonsuu.vec <- c(min(gesui$masuhonsuu), max(gesui$masuhonsuu), 0, 1,
             masuhonsuu.mean - 1/2, masuhonsuu.mean + 1/2,
              masuhonsuu.mean - masuhonsuu.sd/2, masuhonsuu.mean + masuhonsuu.sd/2)




nendo.mean <- mean(gesui$nendo)  # nendoの平均値
nendo.sd <- sd(gesui$nendo)      # nendoのsd
# nendoの8つの値
nendo.vec <- c(min(gesui$nendo), max(gesui$nendo), 0, 1,
             nendo.mean - 1/2, nendo.mean + 1/2,
              nendo.mean - nendo.sd/2, nendo.mean + nendo.sd/2)


kei.mean <- mean(gesui$kei)  # keiの平均値
kei.sd <- sd(gesui$kei)      # keiのsd
# keiの8つの値
kei.vec <- c(min(gesui$kei), max(gesui$kei), 0, 1,
             kei.mean - 1/2, kei.mean + 1/2,
              kei.mean - kei.sd/2, kei.mean + kei.sd/2)


kubun.mean <- mean(gesui$kubun)  # kubunの平均値
kubun.sd <- sd(gesui$kubun)      # kubunのsd
# kubunの8つの値
kubun.vec <- c(min(gesui$kubun), max(gesui$kubun), 0, 1,
             kubun.mean - 1/2, kubun.mean + 1/2,
              kubun.mean - kubun.sd/2, kubun.mean + kubun.sd/2)


did.mean <- mean(gesui$did)  # didの平均値
did.sd <- sd(gesui$did)      # didのsd
# didの8つの値
did.vec <- c(min(gesui$did), max(gesui$did), 0, 1,
             did.mean - 1/2, did.mean + 1/2,
              did.mean - did.sd/2, did.mean + did.sd/2)


kouhou.mean <- mean(gesui$kouhou)  # kouhouの平均値
kouhou.sd <- sd(gesui$kouhou)      # kouhouのsd
# kouhouの8つの値
kouhou.vec <- c(min(gesui$kouhou), max(gesui$kouhou), 0, 1,
             kouhou.mean - 1/2, kouhou.mean + 1/2,
              kouhou.mean - kouhou.sd/2, kouhou.mean + kouhou.sd/2)


ekijyouka.mean <- mean(gesui$ekijyouka)  # ekijyoukaの平均値
ekijyouka.sd <- sd(gesui$ekijyouka)      # ekijyoukaのsd
# ekijyoukaの8つの値
ekijyouka.vec <- c(min(gesui$ekijyouka), max(gesui$ekijyouka), 0, 1,
             ekijyouka.mean - 1/2, ekijyouka.mean + 1/2,
              ekijyouka.mean - ekijyouka.sd/2, ekijyouka.mean + ekijyouka.sd/2)



kansyu.mean <- mean(gesui$kansyu
                   )  # kansyunの平均値
kansyu.sd <- sd(gesui$kansyu)      # kansyuのsd
## kansyuの8つの値
kansyu.vec <- c(min(gesui$kansyu), max(gesui$kansyu), 0, 1,
              kansyu.mean - 1/2, kansyu.mean + 1/2,
              kansyu.mean - kansyu.sd/2, kansyu.mean + kansyu.sd/2)

k <- 8  # 調べる値の数

次に、nendo, kansyu それぞれについて、他方を平均値に固定した上で、8つの値を取ったときの緊急度2以上（5年以内で破損する）の確率の予測値を計算する。

関数repは、かなり柔軟な方法でベクトルを繰り返すために使用 predict関数は予測モデルmodel_1でnewdataで指定したデータに対して実行する。 type = “response”を加えると、ロジスティック関数に代入すると応答変数の予測の確率値が求まる。 https://www1.doshisha.ac.jp/~mjin/R/47/47.html

## slopeの値を変えて当選確率を予測する：その他のパラメータは平均値に固定
wl.slope <- predict(model_1, type = "response",
newdata = list( slope = slope.vec,
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))



## longの値を変えて当選確率を予測する：その他のパラメータは平均値に固定
wl.long <- predict(model_1, type = "response",
newdata = list( long = long.vec,
slope = rep(slope.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))


## uedokaburiの値を変えて当選確率を予測する：その他のパラメータは平均値に固定
wl.uedokaburi <- predict(model_1, type = "response",
newdata = list( uedokaburi = uedokaburi.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))

## sitadokaburiの値を変えて当選確率を予測する：他平均値に固定
wl.sitadokaburi <- predict(model_1, type = "response",
newdata = list( sitadokaburi = sitadokaburi.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))


## masuhonsuuの値を変えて当選確率を予測する：他平均値に固定

wl.masuhonsuu <- predict(model_1, type = "response",
newdata = list( masuhonsuu = masuhonsuu.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))

## nendoの値を変えて当選確率を予測する：他平均値に固定
wl.nendo <- predict(model_1, type = "response",
newdata = list( nendo = nendo.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))


## keiの値を変えて当選確率を予測する：他平均値に固定
wl.kei <- predict(model_1, type = "response",
newdata = list( kei = kei.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))


## kubunの値を変えて当選確率を予測する：他平均値に固定
wl.kubun <- predict(model_1, type = "response",
newdata = list( kubun = kubun.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))

## didの値を変えて当選確率を予測する：他平均値に固定
wl.did <- predict(model_1, type = "response",
newdata = list( did = did.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))


## kouhouの値を変えて当選確率を予測する：他平均値に固定
wl.kouhou <- predict(model_1, type = "response",
newdata = list( kouhou = kouhou.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
ekijyouka = rep(ekijyouka.mean, k),
kansyu = rep(kansyu.mean, k)
))


## ekijyoukaの値を変えて当選確率を予測する：他平均値に固定
wl.ekijyouka <- predict(model_1, type = "response",
newdata = list( ekijyouka = ekijyouka.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
kansyu = rep(kansyu.mean, k)
))


## kansyuの値を変えて当選確率を予測する：他平均値に固定
wl.kansyu <- predict(model_1, type = "response",
newdata = list( kansyu = kansyu.vec, 
slope = rep(slope.mean, k),
long = rep(long.mean, k),
uedokaburi = rep(uedokaburi.mean, k),
sitadokaburi = rep(sitadokaburi.mean, k),
masuhonsuu = rep(masuhonsuu.mean, k),
nendo = rep(nendo.mean, k),
kei = rep(kei.mean, k),
kubun = rep(kubun.mean, k),
did = rep(did.mean, k),
kouhou = rep(kouhou.mean, k),
ekijyouka = rep(ekijyouka.mean, k)
))

結果をまとめる

res <- rbind(wl.slope, wl.long,wl.uedokaburi ,wl.sitadokaburi ,wl.masuhonsuu ,wl.nendo ,wl.kei ,wl.kubun ,wl.did
             ,wl.kouhou , wl.ekijyouka,wl.kansyu)
res

##                           1            2            3            4           5
## wl.slope        0.039045092 0.0051585421 3.904509e-02 3.182284e-02 0.021014992
## wl.long         0.021769612 0.0135772231 2.294937e-02 2.177419e-02 0.019447322
## wl.uedokaburi   0.007931464 0.2318693507 5.956147e-03 7.886951e-03 0.016485429
## wl.sitadokaburi 0.073612763 0.0003922474 7.363492e-02 4.713295e-02 0.023890042
## wl.masuhonsuu   0.015548655 0.0879136046 1.554865e-02 1.765514e-02 0.017777435
## wl.nendo        0.002137239 0.0929976989 2.220446e-16 2.220446e-16 0.016820686
## wl.kei          0.007810324 0.1974566699 4.872546e-03 4.995806e-02 0.005856925
## wl.kubun        0.070911765 0.0143166218 7.091177e-02 1.431662e-02 0.042382326
## wl.did          0.001918319 0.0303980120 1.918319e-03 3.039801e-02 0.004757666
## wl.kouhou       0.016727160 0.0214198893 1.672716e-02 2.141989e-02 0.016735836
## wl.ekijyouka    0.022055865 0.0133313940 2.205587e-02 1.333139e-02 0.024336957
## wl.kansyu       0.001844431 0.0490070597 1.844431e-03 4.900706e-02 0.003642672
##                           6           7          8
## wl.slope        0.017068568 0.023622502 0.01517326
## wl.long         0.018448092 0.019884571 0.01804185
## wl.uedokaburi   0.021754811 0.014153379 0.02530720
## wl.sitadokaburi 0.015001832 0.027993139 0.01277790
## wl.masuhonsuu   0.020179681 0.016588540 0.02162029
## wl.nendo        0.021323344 0.009244214 0.03841589
## wl.kei          0.059506201 0.012537251 0.02852190
## wl.kubun        0.008352028 0.025725151 0.01392074
## wl.did          0.072335719 0.011243858 0.03173901
## wl.kouhou       0.021430946 0.017803899 0.02014976
## wl.ekijyouka    0.014723757 0.021260698 0.01687047
## wl.kansyu       0.092524553 0.008951257 0.03963461

最後に、変化量を調べるために差をとる。

table <- cbind(res[,2] - res[,1],
                  res[,4] - res[,3],
                  res[,6] - res[,5],
                  res[,8] - res[,7])
colnames(table) <- c("min -> max", "0 -> 1", "-+1/2", "-+sd/2")

table

##                   min -> max       0 -> 1         -+1/2       -+sd/2
## wl.slope        -0.033886550 -0.007222252 -0.0039464244 -0.008449246
## wl.long         -0.008192389 -0.001175174 -0.0009992296 -0.001842723
## wl.uedokaburi    0.223937887  0.001930805  0.0052693824  0.011153816
## wl.sitadokaburi -0.073220515 -0.026501967 -0.0088882107 -0.015215235
## wl.masuhonsuu    0.072364950  0.002106488  0.0024022468  0.005031754
## wl.nendo         0.090860460  0.000000000  0.0045026575  0.029171678
## wl.kei           0.189646346  0.045085510  0.0536492764  0.015984651
## wl.kubun        -0.056595144 -0.056595144 -0.0340302978 -0.011804407
## wl.did           0.028479693  0.028479693  0.0675780529  0.020495148
## wl.kouhou        0.004692729  0.004692729  0.0046951103  0.002345858
## wl.ekijyouka    -0.008724471 -0.008724471 -0.0096132003 -0.004390224
## wl.kansyu        0.047162629  0.047162629  0.0888818816  0.030683349

#nendo
#exp.vec = seq(min(nendo), max(nendo), length = n())
#exp.vec

以上の結果から変数をnendo,kansyuに絞り込み再度予測モデルを作成し予測→図化する。

logit <- gesui[-1:-5]
logit <- logit[-2:-6]
logit <- logit[-3]
#kansyuがカテゴリー変数に定義
logit$kansyu <- as.factor(logit$kansyu)
logit

##     nendo kansyu taisyo
## 1    1983      1      0
## 2    1982      1      0
## 3    1981      1      0
## 4    1976      1      0
## 5    1992      0      0
## 6    1992      0      0
## 7    1977      1      0
## 8    1992      0      0
## 9    1992      0      0
## 10   1977      1      0
## 11   1977      1      0
## 12   1977      1      0
## 13   1977      1      0
## 14   1983      1      0
## 15   1989      1      1
## 16   1977      1      0
## 17   1977      1      0
## 18   1977      1      0
## 19   1982      1      0
## 20   1982      1      0
## 21   1982      0      0
## 22   1982      0      0
## 23   1982      1      0
## 24   1982      1      0
## 25   1982      1      0
## 26   1982      1      0
## 27   1982      1      0
## 28   1978      1      0
## 29   1978      1      0
## 30   1978      1      1
## 31   1978      1      0
## 32   1978      1      0
## 33   1978      1      0
## 34   1978      1      0
## 35   1978      1      0
## 36   1978      1      0
## 37   1978      1      0
## 38   1979      1      0
## 39   1979      1      0
## 40   1978      1      0
## 41   1978      1      0
## 42   1977      1      0
## 43   1977      1      0
## 44   1977      1      0
## 45   1991      0      0
## 46   1991      0      0
## 47   1992      0      0
## 48   1991      0      0
## 49   1991      0      0
## 50   1991      0      0
## 51   1991      0      0
## 52   1991      0      0
## 53   1992      1      0
## 54   1992      1      0
## 55   1992      1      0
## 56   1992      1      0
## 57   1991      0      0
## 58   1991      1      0
## 59   1991      1      0
## 60   1991      1      0
## 61   1991      1      0
## 62   1991      1      1
## 63   1991      1      0
## 64   1991      1      0
## 65   1991      1      0
## 66   1991      1      0
## 67   1991      1      0
## 68   1991      1      0
## 69   1991      1      1
## 70   1991      1      1
## 71   1991      1      0
## 72   1991      0      0
## 73   1991      1      1
## 74   1991      0      0
## 75   1991      1      0
## 76   1991      1      0
## 77   1991      1      0
## 78   1991      0      0
## 79   1991      0      0
## 80   1991      0      0
## 81   1991      0      0
## 82   1991      0      0
## 83   1991      0      0
## 84   1991      0      0
## 85   1981      1      0
## 86   1981      1      0
## 87   1981      1      0
## 88   1981      1      0
## 89   1981      1      0
## 90   1979      1      0
## 91   1979      1      0
## 92   1979      1      0
## 93   1979      1      0
## 94   1979      1      0
## 95   1981      1      0
## 96   1981      1      0
## 97   1981      1      1
## 98   1981      1      0
## 99   1981      1      1
## 100  1981      1      1
## 101  1981      1      0
## 102  1981      1      0
## 103  1981      1      1
## 104  1981      1      0
## 105  1981      1      1
## 106  1983      1      1
## 107  1983      1      0
## 108  1983      1      1
## 109  1977      1      0
## 110  1976      1      0
## 111  1980      1      0
## 112  1982      1      0
## 113  1981      1      0
## 114  1981      1      0
## 115  1981      1      0
## 116  1976      1      0
## 117  1976      1      0
## 118  1981      1      1
## 119  1981      1      0
## 120  1981      1      0
## 121  1982      1      0
## 122  1976      1      0
## 123  1976      1      0
## 124  1976      1      0
## 125  1981      1      0
## 126  1981      1      0
## 127  1982      1      0
## 128  1982      1      0
## 129  1982      1      0
## 130  1982      1      0
## 131  1982      1      0
## 132  1982      1      0
## 133  1981      1      0
## 134  1981      1      0
## 135  1982      1      0
## 136  1982      1      0
## 137  1982      1      0
## 138  1982      1      0
## 139  1982      1      0
## 140  1982      1      0
## 141  1982      1      0
## 142  1982      1      0
## 143  1983      1      0
## 144  1983      1      0
## 145  1983      1      0
## 146  1982      1      0
## 147  1982      1      0
## 148  1982      1      0
## 149  1990      0      0
## 150  1990      0      0
## 151  1990      0      0
## 152  1990      0      0
## 153  1981      1      0
## 154  1981      1      0
## 155  1992      1      0
## 156  1992      1      0
## 157  1992      1      0
## 158  1992      1      0
## 159  1992      1      0
## 160  1976      1      0
## 161  1976      1      0
## 162  1976      1      0
## 163  1992      0      0
## 164  1992      0      0
## 165  1992      0      0
## 166  1992      0      0
## 167  1992      0      0
## 168  1977      1      0
## 169  1976      1      1
## 170  1990      1      0
## 171  1982      1      0
## 172  1982      1      0
## 173  1982      1      0
## 174  1982      1      0
## 175  1982      1      0
## 176  1982      1      0
## 177  1982      1      0
## 178  1977      1      0
## 179  1977      1      0
## 180  1977      1      0
## 181  1977      1      0
## 182  1977      1      0
## 183  1977      1      0
## 184  1977      1      0
## 185  1977      1      0
## 186  1977      1      0
## 187  1977      1      0
## 188  1977      1      0
## 189  1986      1      0
## 190  1986      1      0
## 191  1991      1      0
## 192  1977      1      0
## 193  1977      1      0
## 194  1977      1      0
## 195  1977      1      0
## 196  1977      1      0
## 197  1977      1      0
## 198  1977      1      0
## 199  1977      1      0
## 200  1977      1      0
## 201  1977      1      0
## 202  1986      1      0
## 203  1979      1      0
## 204  1979      1      0
## 205  1983      1      0
## 206  1983      1      0
## 207  1982      1      0
## 208  1982      1      0
## 209  1982      1      0
## 210  1982      1      0
## 211  1982      1      0
## 212  1979      1      0
## 213  1981      1      0
## 214  1981      1      0
## 215  1989      1      0
## 216  1989      0      0
## 217  1989      1      0
## 218  1989      1      0
## 219  1989      1      0
## 220  1989      1      0
## 221  1979      1      0
## 222  1985      1      0
## 223  1977      1      0
## 224  1977      1      0
## 225  1977      1      0
## 226  1977      1      0
## 227  1977      1      0
## 228  1977      1      0
## 229  1981      1      0
## 230  1981      1      0
## 231  1977      1      0
## 232  1982      1      0
## 233  1982      1      0
## 234  1982      1      0
## 235  1982      1      0
## 236  1982      1      0
## 237  1982      1      0
## 238  1982      1      0
## 239  1982      1      0
## 240  1982      1      0
## 241  1982      1      0
## 242  1982      1      0
## 243  1982      1      0
## 244  1982      1      0
## 245  1982      1      0
## 246  1977      1      0
## 247  1977      1      0
## 248  1977      1      0
## 249  1977      1      0
## 250  1977      1      0
## 251  1978      1      0
## 252  1978      1      0
## 253  1978      1      0
## 254  1978      1      0
## 255  1978      1      0
## 256  1978      1      0
## 257  1978      1      0
## 258  1978      1      0
## 259  1978      1      0
## 260  1978      1      0
## 261  1978      1      0
## 262  1978      1      0
## 263  1978      1      0
## 264  1978      1      0
## 265  1978      1      0
## 266  1978      1      1
## 267  1978      1      0
## 268  1978      1      0
## 269  1978      1      0
## 270  1978      1      0
## 271  1978      1      0
## 272  1978      1      0
## 273  1978      1      0
## 274  1977      1      0
## 275  1977      1      0
## 276  1977      1      0
## 277  1980      1      0
## 278  1978      1      0
## 279  1978      1      0
## 280  1978      1      0
## 281  1978      1      0
## 282  1977      1      0
## 283  1977      1      0
## 284  1977      1      0
## 285  1983      1      0
## 286  1983      1      1
## 287  1983      1      0
## 288  1980      1      0
## 289  1989      0      0
## 290  1989      0      0
## 291  1989      0      0
## 292  1989      0      0
## 293  1980      1      0
## 294  1989      0      0
## 295  1989      0      0
## 296  1989      0      0
## 297  1989      0      0
## 298  1989      0      0
## 299  1989      0      0
## 300  1989      0      0
## 301  1989      1      0
## 302  1989      0      0
## 303  1989      0      0
## 304  1989      0      1
## 305  1989      0      0
## 306  1989      0      0
## 307  1989      0      0
## 308  1989      0      0
## 309  1989      1      0
## 310  1989      1      0
## 311  1989      0      0
## 312  1989      0      0
## 313  1989      0      0
## 314  1977      1      0
## 315  1977      1      0
## 316  1977      1      0
## 317  1977      1      0
## 318  1977      1      0
## 319  1977      1      0
## 320  1977      1      0
## 321  1977      1      0
## 322  1977      1      0
## 323  1977      1      0
## 324  1991      0      0
## 325  1991      0      0
## 326  1991      0      0
## 327  1991      0      0
## 328  1991      0      0
## 329  1991      0      0
## 330  1991      0      0
## 331  1991      0      0
## 332  1991      0      0
## 333  1991      0      0
## 334  1991      0      0
## 335  1990      0      0
## 336  1990      0      0
## 337  1990      0      0
## 338  1990      0      0
## 339  1990      0      0
## 340  1991      0      0
## 341  1990      0      0
## 342  1990      0      0
## 343  1990      0      0
## 344  1990      0      0
## 345  1977      1      0
## 346  1977      1      0
## 347  1977      1      0
## 348  1977      1      0
## 349  1977      1      0
## 350  1977      1      0
## 351  1977      1      0
## 352  1977      1      0
## 353  1977      1      0
## 354  1991      0      0
## 355  1991      0      0
## 356  1991      0      0
## 357  1991      0      0
## 358  1991      0      0
## 359  1991      0      0
## 360  1991      0      0
## 361  1991      0      0
## 362  1991      0      0
## 363  1991      0      0
## 364  1991      0      0
## 365  1991      0      0
## 366  1991      0      0
## 367  1991      0      0
## 368  1991      0      0
## 369  1991      0      0
## 370  1992      1      0
## 371  1992      1      0
## 372  1992      1      0
## 373  1991      0      0
## 374  1991      0      0
## 375  1991      0      0
## 376  1991      0      0
## 377  1991      0      0
## 378  1991      0      0
## 379  1991      0      0
## 380  1991      0      0
## 381  1991      0      0
## 382  1991      0      0
## 383  1991      0      0
## 384  1991      0      0
## 385  1992      1      0
## 386  1992      1      0
## 387  1992      1      0
## 388  1992      1      0
## 389  1992      1      0
## 390  1992      1      0
## 391  1992      1      0
## 392  1992      1      0
## 393  1992      1      0
## 394  1992      1      0
## 395  1992      1      0
## 396  1992      1      0
## 397  1992      1      0
## 398  1992      1      0
## 399  1992      1      0
## 400  1992      1      0
## 401  1991      0      0
## 402  1991      0      0
## 403  1992      0      0
## 404  1992      1      0
## 405  1992      1      0
## 406  1992      1      0
## 407  1992      1      0
## 408  1992      1      0
## 409  1991      1      0
## 410  1991      1      1
## 411  1991      1      0
## 412  1991      1      1
## 413  1991      1      0
## 414  1991      1      1
## 415  1991      1      0
## 416  1991      1      1
## 417  1991      1      0
## 418  1991      1      0
## 419  1991      1      1
## 420  1991      1      0
## 421  1991      1      0
## 422  1991      1      0
## 423  1991      1      1
## 424  1991      1      0
## 425  1991      1      1
## 426  1991      1      0
## 427  1991      1      1
## 428  1991      1      0
## 429  1991      1      0
## 430  1991      1      0
## 431  1991      1      0
## 432  1991      1      1
## 433  1991      1      0
## 434  1991      1      1
## 435  1991      1      0
## 436  1991      1      0
## 437  1991      1      0
## 438  1991      1      0
## 439  1991      1      0
## 440  1991      0      0
## 441  1991      0      0
## 442  1991      0      0
## 443  1991      1      0
## 444  1992      0      0
## 445  1991      0      0
## 446  1991      0      0
## 447  1991      0      0
## 448  1991      0      0
## 449  1991      0      0
## 450  1991      1      0
## 451  1991      1      0
## 452  1991      1      0
## 453  1991      1      0
## 454  1991      0      0
## 455  1991      0      0
## 456  1991      0      0
## 457  1991      0      0
## 458  1991      0      0
## 459  1991      0      0
## 460  1991      0      0
## 461  1991      0      0
## 462  1991      0      0
## 463  1991      0      0
## 464  1991      0      0
## 465  1991      0      0
## 466  1991      0      0
## 467  1991      0      0
## 468  1991      0      0
## 469  1991      0      0
## 470  1991      0      0
## 471  1991      0      0
## 472  1991      0      0
## 473  1991      0      0
## 474  1991      0      0
## 475  1991      0      0
## 476  1991      0      0
## 477  1991      0      0
## 478  1991      0      0

model_2 <- glm(taisyo ~ nendo + kansyu, data = logit,
               family = binomial(link = "logit"))
summary(model_2)

## 
## Call:
## glm(formula = taisyo ~ nendo + kansyu, family = binomial(link = "logit"), 
##     data = logit)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6009  -0.3619  -0.2787  -0.1217   3.2008  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -218.46755   68.96078  -3.168  0.00153 **
## nendo          0.10727    0.03464   3.097  0.00196 **
## kansyu1        3.17478    1.03742   3.060  0.00221 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 213.22  on 477  degrees of freedom
## Residual deviance: 190.01  on 475  degrees of freedom
## AIC: 196.01
## 
## Number of Fisher Scoring iterations: 7

ロジスティック回帰式の精度-予測値を計算

predict(model_2, type = "response")

##           1           2           3           4           5           6 
## 0.070078832 0.063402850 0.057323644 0.034345442 0.008204210 0.008204210 
##           7           8           9          10          11          12 
## 0.038086251 0.008204210 0.008204210 0.038086251 0.038086251 0.038086251 
##          13          14          15          16          17          18 
## 0.038086251 0.070078832 0.125440738 0.038086251 0.038086251 0.038086251 
##          19          20          21          22          23          24 
## 0.063402850 0.063402850 0.002821886 0.002821886 0.063402850 0.063402850 
##          25          26          27          28          29          30 
## 0.063402850 0.063402850 0.063402850 0.042216685 0.042216685 0.042216685 
##          31          32          33          34          35          36 
## 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 
##          37          38          39          40          41          42 
## 0.042216685 0.046773282 0.046773282 0.042216685 0.042216685 0.038086251 
##          43          44          45          46          47          48 
## 0.038086251 0.038086251 0.007375892 0.007375892 0.008204210 0.007375892 
##          49          50          51          52          53          54 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.165192402 0.165192402 
##          55          56          57          58          59          60 
## 0.165192402 0.165192402 0.007375892 0.150926113 0.150926113 0.150926113 
##          61          62          63          64          65          66 
## 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 
##          67          68          69          70          71          72 
## 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 0.007375892 
##          73          74          75          76          77          78 
## 0.150926113 0.007375892 0.150926113 0.150926113 0.150926113 0.007375892 
##          79          80          81          82          83          84 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##          85          86          87          88          89          90 
## 0.057323644 0.057323644 0.057323644 0.057323644 0.057323644 0.046773282 
##          91          92          93          94          95          96 
## 0.046773282 0.046773282 0.046773282 0.046773282 0.057323644 0.057323644 
##          97          98          99         100         101         102 
## 0.057323644 0.057323644 0.057323644 0.057323644 0.057323644 0.057323644 
##         103         104         105         106         107         108 
## 0.057323644 0.057323644 0.057323644 0.070078832 0.070078832 0.070078832 
##         109         110         111         112         113         114 
## 0.038086251 0.034345442 0.051795092 0.063402850 0.057323644 0.057323644 
##         115         116         117         118         119         120 
## 0.057323644 0.034345442 0.034345442 0.057323644 0.057323644 0.057323644 
##         121         122         123         124         125         126 
## 0.063402850 0.034345442 0.034345442 0.034345442 0.057323644 0.057323644 
##         127         128         129         130         131         132 
## 0.063402850 0.063402850 0.063402850 0.063402850 0.063402850 0.063402850 
##         133         134         135         136         137         138 
## 0.057323644 0.057323644 0.063402850 0.063402850 0.063402850 0.063402850 
##         139         140         141         142         143         144 
## 0.063402850 0.063402850 0.063402850 0.063402850 0.070078832 0.070078832 
##         145         146         147         148         149         150 
## 0.070078832 0.063402850 0.063402850 0.063402850 0.006630644 0.006630644 
##         151         152         153         154         155         156 
## 0.006630644 0.006630644 0.057323644 0.057323644 0.165192402 0.165192402 
##         157         158         159         160         161         162 
## 0.165192402 0.165192402 0.165192402 0.034345442 0.034345442 0.034345442 
##         163         164         165         166         167         168 
## 0.008204210 0.008204210 0.008204210 0.008204210 0.008204210 0.038086251 
##         169         170         171         172         173         174 
## 0.034345442 0.137688675 0.063402850 0.063402850 0.063402850 0.063402850 
##         175         176         177         178         179         180 
## 0.063402850 0.063402850 0.063402850 0.038086251 0.038086251 0.038086251 
##         181         182         183         184         185         186 
## 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 
##         187         188         189         190         191         192 
## 0.038086251 0.038086251 0.094175730 0.094175730 0.150926113 0.038086251 
##         193         194         195         196         197         198 
## 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 
##         199         200         201         202         203         204 
## 0.038086251 0.038086251 0.038086251 0.094175730 0.046773282 0.046773282 
##         205         206         207         208         209         210 
## 0.070078832 0.070078832 0.063402850 0.063402850 0.063402850 0.063402850 
##         211         212         213         214         215         216 
## 0.063402850 0.046773282 0.057323644 0.057323644 0.125440738 0.005960242 
##         217         218         219         220         221         222 
## 0.125440738 0.125440738 0.125440738 0.125440738 0.046773282 0.085415035 
##         223         224         225         226         227         228 
## 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 
##         229         230         231         232         233         234 
## 0.057323644 0.057323644 0.038086251 0.063402850 0.063402850 0.063402850 
##         235         236         237         238         239         240 
## 0.063402850 0.063402850 0.063402850 0.063402850 0.063402850 0.063402850 
##         241         242         243         244         245         246 
## 0.063402850 0.063402850 0.063402850 0.063402850 0.063402850 0.038086251 
##         247         248         249         250         251         252 
## 0.038086251 0.038086251 0.038086251 0.038086251 0.042216685 0.042216685 
##         253         254         255         256         257         258 
## 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 
##         259         260         261         262         263         264 
## 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 
##         265         266         267         268         269         270 
## 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 0.042216685 
##         271         272         273         274         275         276 
## 0.042216685 0.042216685 0.042216685 0.038086251 0.038086251 0.038086251 
##         277         278         279         280         281         282 
## 0.051795092 0.042216685 0.042216685 0.042216685 0.042216685 0.038086251 
##         283         284         285         286         287         288 
## 0.038086251 0.038086251 0.070078832 0.070078832 0.070078832 0.051795092 
##         289         290         291         292         293         294 
## 0.005960242 0.005960242 0.005960242 0.005960242 0.051795092 0.005960242 
##         295         296         297         298         299         300 
## 0.005960242 0.005960242 0.005960242 0.005960242 0.005960242 0.005960242 
##         301         302         303         304         305         306 
## 0.125440738 0.005960242 0.005960242 0.005960242 0.005960242 0.005960242 
##         307         308         309         310         311         312 
## 0.005960242 0.005960242 0.125440738 0.125440738 0.005960242 0.005960242 
##         313         314         315         316         317         318 
## 0.005960242 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 
##         319         320         321         322         323         324 
## 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 0.007375892 
##         325         326         327         328         329         330 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         331         332         333         334         335         336 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.006630644 0.006630644 
##         337         338         339         340         341         342 
## 0.006630644 0.006630644 0.006630644 0.007375892 0.006630644 0.006630644 
##         343         344         345         346         347         348 
## 0.006630644 0.006630644 0.038086251 0.038086251 0.038086251 0.038086251 
##         349         350         351         352         353         354 
## 0.038086251 0.038086251 0.038086251 0.038086251 0.038086251 0.007375892 
##         355         356         357         358         359         360 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         361         362         363         364         365         366 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         367         368         369         370         371         372 
## 0.007375892 0.007375892 0.007375892 0.165192402 0.165192402 0.165192402 
##         373         374         375         376         377         378 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         379         380         381         382         383         384 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         385         386         387         388         389         390 
## 0.165192402 0.165192402 0.165192402 0.165192402 0.165192402 0.165192402 
##         391         392         393         394         395         396 
## 0.165192402 0.165192402 0.165192402 0.165192402 0.165192402 0.165192402 
##         397         398         399         400         401         402 
## 0.165192402 0.165192402 0.165192402 0.165192402 0.007375892 0.007375892 
##         403         404         405         406         407         408 
## 0.008204210 0.165192402 0.165192402 0.165192402 0.165192402 0.165192402 
##         409         410         411         412         413         414 
## 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 
##         415         416         417         418         419         420 
## 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 
##         421         422         423         424         425         426 
## 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 
##         427         428         429         430         431         432 
## 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 
##         433         434         435         436         437         438 
## 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 0.150926113 
##         439         440         441         442         443         444 
## 0.150926113 0.007375892 0.007375892 0.007375892 0.150926113 0.008204210 
##         445         446         447         448         449         450 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.150926113 
##         451         452         453         454         455         456 
## 0.150926113 0.150926113 0.150926113 0.007375892 0.007375892 0.007375892 
##         457         458         459         460         461         462 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         463         464         465         466         467         468 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         469         470         471         472         473         474 
## 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 0.007375892 
##         475         476         477         478 
## 0.007375892 0.007375892 0.007375892 0.007375892

stargazer(model_2, 
          style = "ajps", # AJPS仕様
          dep.var.labels = "Win or Lose", # 応答変数の名前
          title = "Logistic Regression Analysis on Winning a Seat", # タイトル
          type ="html")

Logistic Regression Analysis on Winning a Seat

	Win or Lose
nendo	0.107***
	(0.035)
kansyu1	3.175***
	(1.037)
Constant	-218.468***
	(68.961)
N	478
Log Likelihood	-95.006
AIC	196.013
p < .01; p < .05; p < .1

nendo、kansyu共に1%で統計的に有意であることがわかる

説明変数を複数の値に固定したグラフを作成する

それぞれについて、他方を平均値に固定した上で、8つの値を取ったときの予測値を計算する。

nendo.mean <- mean(logit\(nendo) # nendoの平均値 nendo.sd <- sd(logit\)nendo) # nendoのsd # nendoの8つの値 nendo.vec <- c(min(logit\(nendo), max(logit\)nendo), 0, 1, nendo.mean - 1/2, nendo.mean + 1/2, nendo.mean - nendo.sd/2, nendo.mean + nendo.sd/2) nendo.vec

kansyu.mean <- mean(logit\(kansyu) # kansyunの平均値 kansyu.sd <- sd(logit\)kansyu) # kansyuのsd ## kansyuの8つの値 kansyu.vec <- c(min(logit\(kansyu), max(logit\)kansyu), 0, 1, kansyu.mean - 1/2, kansyu.mean + 1/2, kansyu.mean - kansyu.sd/2, kansyu.mean + kansyu.sd/2) kansyu.vec

k <- 8 # 調べる値の数

他方を平均値に固定した上で、8つの値を取ったときの当選確率の予測値を計算する。

wl.nendo <- predict(model_2, type = “response”, newdata = list( nendo = nendo.vec, kansyu = rep(kansyu.mean, k))) ## kansyuの値を変えて当選確率を予測する：nendoiousは平均値に固定 wl.kansyu <- predict(model_2, type = “response”, newdata = list( nendo = rep(nendo.mean, k), kansyu = kansyu.vec)) ## 結果をまとめる res <- rbind(wl.nendo, wl.kansyu)

mutate関数：データセットに新たに変数を追加する関数 https://kazutan.github.io/kazutanR/hands_on_170730/mutate.html

logit <- mutate(logit, exp.vec = seq(min(nendo), max(nendo), length = n()), p.hat.wl.prev0 = predict(model_2, type = “response”, list(nendo = rep(0, n()), kansyu = exp.vec)))

wl.exp.lab <- labs(x = “Campaign Money (Million yen)”, y = “Probability of Winning”) p.hat.0 <- ggplot(logit, aes(exp.vec, p.hat.wl.prev0)) + theme_bw() + geom_line(color = “red”) p.hat.0 + wl.exp.lab + ggtitle(“When previous = 0”)