R Logistic Regression (ロジスティック回帰)

Program R Mathematics Logistic Regression ロジスティック回帰 
#ロジスティクス回帰分析は、1つのカテゴリカル変数を目的変数
#とし、その目的変数を説明変数で説明するモデルによる分析方法。

#より実践的なのはこちら

################目的変数が2値の場合################
##単項多重ロジスティック回帰
#個票データが与えられている場合
library(vcd)
head(SpaceShuttle)
glm(formula = Fail ~ Temperature, data = SpaceShuttle, family = binomial)

#集計データが与えられている場合
library(MASS)
head(housing) #集計データ
glm(formula = Sat ~ Cont, data = housing, weight = Freq, family = binomial)

#カテゴリーごとにYesとNoの人数がわかる場合
dose <- c(1:5)
effect <- c(2, 3, 9, 15, 18)
effect.prop <- effect/20
plot(dose, effect.prop)
glm(cbind(effect, 20 - effect)~dose,family=binomial)

##多重ロジスティック回帰
library(vcd)
result <- glm(Fail ~ Temperature + Pressure, data = SpaceShuttle, family=binomial)
result
summary(result) #係数に関する検定や信頼区間の情報を知りたいとき

#Residual Devienceは対数尤度
#1 observation delited due to missingnessは欠損値の数を示す
#複数の説明変数が考えられる場合には、どの説明変数を用いたモデルが
#良いのか、ということが問題となる。AICがその指標としてよく用いられる。
#   AIC = -2(最大対数尤度)+2(モデルに含まれるパラメータ数)

##ステップワイズ法
#多重ロジスティック回帰においてモデルの選択を考える場合には、用いる
#変数の組み合わせに対応して統計モデルが決まる。そのため、説明変数の個数が
#大きくなる組み合わせが指数関数的に大きくなり、それらをすべて解析して
#モデルを選択することは難しい。
#ステップワイズ法には、説明変数の数が0の場合からスタートして、
#だんだん説明変数を増やしていく変数増加法と、すべての説明変数を用いた
#モデルからスタートして、だんだん変数を減らしていく変数減少法がある。

step(glm(Fail ~ Temperature + Pressure, data = SpaceShuttle,         
         family = binomial))

#欠損値がある場合は、AICではなくBIC(Bayesian Information Criteria)を
#用いることもある。
#BIC = -2(最大対数尤度)+2log(有効データ数)×(パラメータ数)

step(glm(Fail ~ Temperature + Pressure, data = SpaceShuttle,
         family = binomial), k = log(23))

#23は有効データ数
#最後に表示されるAICはAIC、途中のAICはBICとなる

################目的変数がカテゴリカルの場合################
##カテゴリー間に順序が無い場合
summary(iris)
library(nnet)
multinom(formula = Species ~ Sepal.Length, data = iris)

#あるひとつのカテゴリを基準として、他のカテゴリとの対数オッズ
#を推定する。

##カテゴリー間に順序がある場合
summary(housing)
head(housing)
Ctable <- xtabs(Freq~Cont+Sat, data=housing)
mosaicplot(Ctable,col=TRUE)

library(MASS)
polr(Sat ~ Cont, data = housing, weight = Freq)

#家主の満足度(Sat)をLow+Medium,High/Low,Medium+Highで分けた
#場合で、グループ間のオッズを表すパラメータを推定する。

##条件付きロジスティック回帰分析
head(infert)
infert[74,]
Ctable <- xtabs(~case+parity+education,data=infert)
Ctable

library(survival)
clogit(case ~ spontaneous+induced+strata(stratum),data=infert)
infert
head(infert)

################他のデータで練習################
#1
class(Titanic)
TN1 <- data.frame(Titanic)
TN2 <- data.frame(lapply(TN1, function(i) rep(i, TN1[,5])))
TN3 <- TN2[,-5]

head(TN3)
TN4 <- glm(Survived ~ Class+Sex+Age, data=TN3, family=binomial)
summary(TN4)
cbind(TN1, SurviveRate = predict(TN4, Titanic, type ="response"))
exp(TN4$coefficient)

#2
library(vcd)
head(JointSports)
nrow(JointSports)
xtabs(~opinion+year, data = JointSports)
JS <- polr(formula = opinion ~ year + grade + gender, data = JointSports)
summary(JS)

##############SASデータを取り込んで実験##############
install.packages("sas7bdat")
library(sas7bdat)

#SASデータはこちらから(http://nakhirot.hatenablog.com/entry/2013/05/01/220451)
resp <- read.sas7bdat("resp_score.sas7bdat")

names(resp)
delete <- c(18:26)

resp1 <- resp[,-delete]
resp.LR <- glm(resp~.,family=binomial,data=resp1)
summary(resp.LR)
resp.LR1 <- summary(step(resp.LR)) #説明変数が多すぎるのでstepで回数削減
exp(resp.LR1$coefficients) #オッズを確認