axjack's blog

### axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz ###

統計学実践ワークブック 第18章 質的回帰の問18.1-3のロジスティック回帰分析/プロビットモデルをRで実行

第18章 質的回帰

  • データの引用元
    • 統計学実践ワークブック pp.152-153
# 問18.1 ####
# データ
LIs    <- c(8,10,12,14,16,18,20,22,24,26,28,32,34,38)
Kanja  <- c(2,2,3,3,3,1,3,2,1,1,1,1,1,3)
Kankai <- c(0,0,0,0,0,1,2,1,0,1,1,0,1,2)

# 加工
LI     <- rep(LIs,Kanja)
resp   <- mapply( (\(x,y) c(rep(1,x),rep(0,y-x))), Kankai, Kanja ) |> unlist() 
df     <- data.frame(LI,resp)

# ロジスティック回帰
glm.logistic <- glm(df$resp ~ ., family=binomial(link="logit"), data = df)
summary(glm.logistic)

# [1]の確率
predict(glm.logistic, newdata = data.frame(LI=30), type = "response")



# 問18.2
# データ
Smokings  <- c(0,1,0,1,0,1,0,1)
Obesitys  <- c(0,0,1,1,0,0,1,1)
Snorings  <- c(0,0,0,0,1,1,1,1)
reisu     <- c(60,17,8,2,187,85,51,23)
kouketsus <- c(5,2,1,0,35,13,15,8)

# 加工
smoking   <- rep(Smokings, reisu)
obesity   <- rep(Obesitys, reisu)
snoring   <- rep(Snorings, reisu)
resp      <- mapply( (\(x,y) c(rep(1,x),rep(0,y-x))), kouketsus, reisu ) |> unlist() 
df2       <- data.frame(smoking,obesity,snoring,resp)

# ロジスティック回帰
glm.logistic2 <- glm(resp ~., family = binomial(link = "logit"), data = df2)
summary(glm.logistic2)

# [1]の確率
predict(glm.logistic2, newdata = data.frame(smoking=1,obesity=1,snoring=1),type="response")



# 問18.3
# プロビットモデル
# データは問18.2と同じ
glm.probit <- glm(resp ~., family = binomial(link = "probit"), data = df2)
summary(glm.probit)

# [1]の確率
predict(glm.probit, newdata = data.frame(smoking=1,obesity=1,snoring=1),type="response")

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.