# log cacluation
# j = 1
j <- log10(10)
j
# 10^j = (10)
10^j
k <- log2(10)
k
2^k
set.seed(101)
##########
# see youtube
# https://youtu.be/8nm0G-1uJzA
n.mut <- 23+117
n.norm <- 6+210
p.cancer.mut <- 23/(23+117)
p.cancer.norm <- 6/(6+210)
set.seed(1011)
c <- runif(n.mut, 0, 1)
# 0 = not cancer, 1 = cancer among mutant gene
mutant <- ifelse(c>=p.cancer.mut, 0, 1)
c <- runif(n.norm, 0, 1)
# 0 = not cancer, 1 = cancer among normal gene
normal <- ifelse(c>=p.cancer.norm, 0, 1)
# 0 = mutant; 1 = normal
gene <- c(rep(0, length(mutant)), rep(1, length(normal)))
# 0 = not cancer; 1 = cancer
cancer <- c(mutant, normal)
df <- as.data.frame(cbind(gene, cancer))
df
df$gene <- factor(df$gene, levels = c(0,1), labels = c("mutant", "norm"))
df$cancer <- factor(df$cancer, levels = c(0,1), labels = c("nocancer", "cancer"))
df
tab <- table(df)
tab
tab[1,2]
tab[1,1]
# p.c.m = p.cancer.mut the above
p.cancer.mutant <- tab[1,2]/(tab[1,1]+tab[1,2])
p.nocancer.mutant <- tab[1,1]/(tab[1,1]+tab[1,2])
p.cancer.mutant
1-p.cancer.mutant
p.nocancer.mutant
p.cancer.norm <- tab[2,2]/(tab[2,1]+tab[2,2])
p.nocancer.norm <- 1-p.cancer.norm
p.cancer.norm
p.nocancer.norm
odds(p.cancer.mutant)
odds(p.cancer.norm)
odds.ratio(p.cancer.mutant, p.cancer.norm)
###########################################
load("nsduh2019_adult_sub_rmph.RData")
# Shorter name
nsduh <- nsduh_adult_sub
tab <- table(nsduh$demog_sex, nsduh$mj_lifetime)
tab
table(nsduh$demog_sex)
table(nsduh$mj_lifetime)
# create functions for odd calculation
odds <- function(p) p/(1-p)
odds.ratio <- function(p1, p2) odds(p1)/odds(p2)
# log odds를 구하는 function
logit <- function(p) log(p/(1-p))
# probability를 구하는 function
# log(p/(1-p)) = x
# p/(1-p) = e^x
# p = e^x * (1-p)
# p = e^x - e^x*p
# p + p*e^x = e^x
# p*(1+e^x) = e^x
# p = e^x / (1 +e^x)
# 따라서 아래는 p를 구하는 평션
ilogit <- function(x) exp(x)/(1+exp(x))
# exp() is the exponential function
# 위의 펑션으로 구한 값
PM.yes <- tab[1,2]/(tab[1,2]+tab[1,1])
PF.yes <- tab[2,2]/(tab[2,2]+tab[2,1])
OM <- odds(PM.yes)
OF <- odds(PF.yes)
OR.mvsf <- odds.ratio(PM.yes, PF.yes)
# outputs
round(c(PM.yes, PF.yes, OM, OF, OR.mvsf), 5)
library(dplyr)
# female을 reference로 바꾸는 작업 relevel
nsduh <- nsduh %>%
mutate(demog_sex = relevel(demog_sex, ref = "Female"))
# generalized linear modeling (glm) for logistic regression
# note: IV가 종류측정 변인
fit.ex6.2 <- glm(mj_lifetime ~ demog_sex,
family = binomial, data = nsduh)
summary(fit.ex6.2)
# 아래는 다른 펑션으로 아웃풋이 다름름
# library(lessR)
# fit.ex6.2b <- Logit(mj_lifetime ~ demog_sex, data=nsduh)
# 절편 해석
# log(p) = a + b * demog_sex 에서x
# demog_sexFemale 일 때, b = 0 이므로
# log(p) = a = -0.13504
a = -0.13504
p.femal.yes <- exp(a)/(1 + exp(a)) # 위의 절편에 대한 p 값 계산
p.femal.yes
PF.yes
# x = 1일 때, log(x) = a + b
summary(fit.ex6.2)$coefficient # coefficient 값들 출력
a <- summary(fit.ex6.2)$coefficient[1,1] # 절편 값 출력
b <- summary(fit.ex6.2)$coefficient[2,1] # 절편 값 출력
a
b
adash <- a + b
adash
# or
exp(adash) / (1 + exp(adash))
ilogit(adash)
# 위 값은 PM.yes 값과 같다
PM.yes
# 한편 위에서의 b, coefficient 값은 b 인데
# see http://commres.net/wiki/logistic_regression#odds_ratio_in_logistic
# 아래는 odds ratio 의 값이 된다. (위의 위키문서 읽을 것)
# 즉, x가 한 unit 증가할 때의 odds값의 ratio를 말한다.
exp(b)
# [1] 1.444613
# X변인의 한 unit이 증가함은 female에서 male이 되는 것
# 그 때의 odds ratio는 exp(b), 1.444613
# female에 비해서 male의 me가 약 1.45배 정도 된다는 것것
# 위의 b값에 대한 CI을 구하기 위해서 confint 펑션을 사용한다
confint(fit.ex6.2)
# b값에 대한 것만 알고 싶으므로
confint(fit.ex6.2)[2, , drop=F]
# 그리고 이 값들의 실제 odds ratio값을 보려면
exp(confint(fit.ex6.2)[2,])
# 위는 female에서 male이 될면, 즉, 0 -> 1 이 될때의
# 마리화나 경험의 (me) odds ratio는 CI 범위는
# 1.126 ~ 1.855 즉, 13% 에서 18% 정도 증가한다는 뜻이다
# install.packages("car")
library(car)
# coefficient probability test
car::Anova(fit.ex6.2, type = 3, test.statistic = "Wald")
########################################
########################################
########################################
# numeric IV
fit.ex6.3 <- glm(mj_lifetime ~ alc_agefirst,
family = binomial, data = nsduh)
summary(fit.ex6.3)
round(summary(fit.ex6.3)$coef, 4)
ilogit(coef(fit.ex6.3)["(Intercept)"])
or <- exp(coef(fit.ex6.3)["alc_agefirst"])
1-or
# 0.9952 = prob of starting marijuana
# when the age is zero (intercept이므로)
# age = 0 에서 추정하는 것은 이상함
summary(nsduh$alc_agefirst)
# 위의 아웃풋에서 Mean값이 약 17이므로 17을
# 기준으로 하여 다시 보면
nsduh <- nsduh %>%
mutate(calc_agefirst = alc_agefirst - 17)
fit.ex6.3.centered <- glm(mj_lifetime ~ calc_agefirst,
family = binomial, data = nsduh)
summary(fit.ex6.3.centered)
# 우선 나이가 17살일 때 (x=0)
# 절편 값인 0.5207이 logit 값이 된다.
# 이 때, prob 를 구하면 아래와 같다
p <- ilogit(coef(fit.ex6.3.centered)["(Intercept)"])
p
# 아니면
a <- fit.ex6.3.centered$coefficients["(Intercept)"]
exp(a)/(1+exp(a))
# 독립변인인 나이에 대한 해석
# b coefficient
# 17살일 때를 기준으로 한살씩 증가할 때마다의
# 마리화나 경험/비경험의 Odds ratio는 -0.2835
# 이를 수치화하면 (odds ratio는 exp(b)이다)
exp(coef(fit.ex6.3.centered)["calc_agefirst"])
# 이는 17이후에 한살씩 알콜처음 경험을 늦추면
# 마리화나 경험율 대 미경험 odds ratio가
# 0.247 낮아진다고 할 수 있다 (0.7531 증가는)
# 1-0.7531로 보는 것
# 그리고 이에 대한 CI를 보면 아래와 같고
confint(fit.ex6.3.centered)[2, , drop = F]
# 이를 승비로 (odds ratio) 보면
exp(confint(fit.ex6.3.centered)[2, , drop = F])
#################################
#################################
# 1년이 아니라 3년일 경우
fit.ex6.3.centered
coef(fit.ex6.3.centered)["calc_agefirst"]
coef(fit.ex6.3.centered)["calc_agefirst"]*3
exp(coef(fit.ex6.3.centered)["calc_agefirst"]*3)
# CI 의 경우 아래와 같고
confint(fit.ex6.3.centered)[2,]*3
# 이에 해당하는 값은
exp(confint(fit.ex6.3.centered)[2,]*3)
#################################
#################################
## Multiple regression
#################################
fit.ex6.3.adj <- glm(mj_lifetime ~ alc_agefirst +
demog_age_cat6 + demog_sex +
demog_income,
family = binomial, data = nsduh)
# Regression coefficient table
round(summary(fit.ex6.3.adj)$coef, 4)
coef.values <- coef(fit.ex6.3.adj)
data.frame(coef.values)
logit.values <- exp(coef(fit.ex6.3.adj))
data.frame(logit.values)
confint(fit.ex6.3.adj)
exp(confint(fit.ex6.3.adj))
confint.logit.values <- exp(confint(fit.ex6.3.adj))
# logit.values 를 Adjusted Odds Ratio로 하고 출력
cbind("AdjOR" = logit.values, confint.logit.values)
# 절편은 coefficient 해석과 (odds ratio) 관계 없으므로 생략
cbind("AdjOR" = logit.values, confint.logit.values)[-1,]
# 소수점 3으로 제약
round(cbind("AdjOR" = logit.values,
confint.logit.values)[-1,],3)
# OR은 coefficient 값을 이야기하는 것을 다시 확인
# 또한 wald significant test
# P-values
car::Anova(fit.ex6.3.adj, type = 3, test.statistic = "Wald")
> # log cacluation
> # j = 1
> j <- log10(10)
> j
[1] 1
> # 10^j = (10)
> 10^j
[1] 10
> k <- log2(10)
> k
[1] 3.321928
> 2^k
[1] 10
>
> odds <- function(p) p/(1-p)
> odds.ratio <- function(p1, p2) odds(p1)/odds(p2)
> logit <- function(p) log(p/(1-p))
> ilogit <- function(x) exp(x)/(1+exp(x))
> # exp() is the exponential function
>
> ##########
> # see youtube
> # https://youtu.be/8nm0G-1uJzA
> n.mut <- 23+117
> n.norm <- 6+210
> p.cancer.mut <- 23/(23+117)
> p.cancer.norm <- 6/(6+210)
>
> set.seed(1011)
> c <- runif(n.mut, 0, 1)
> # 0 = not cancer, 1 = cancer among mutant gene
> mutant <- ifelse(c>=p.cancer.mut, 0, 1)
>
> c <- runif(n.norm, 0, 1)
> # 0 = not cancer, 1 = cancer among normal gene
> normal <- ifelse(c>=p.cancer.norm, 0, 1)
>
> # 0 = mutant; 1 = normal
> gene <- c(rep(0, length(mutant)), rep(1, length(normal)))
> # 0 = not cancer; 1 = cancer
> cancer <- c(mutant, normal)
>
> df <- as.data.frame(cbind(gene, cancer))
> head(df)
gene cancer
1 0 0
2 0 1
3 0 0
4 0 0
5 0 0
6 0 0
> df$gene <- factor(df$gene, levels = c(0,1), labels = c("mutant", "norm"))
> df$cancer <- factor(df$cancer, levels = c(0,1), labels = c("nocancer", "cancer"))
> head(df)
gene cancer
1 mutant nocancer
2 mutant cancer
3 mutant nocancer
4 mutant nocancer
5 mutant nocancer
6 mutant nocancer
> tab <- table(df)
> tab
cancer
gene nocancer cancer
mutant 121 19
norm 210 6
> tab[1,2]
[1] 19
> tab[1,1]
[1] 121
>
> # p.c.m = p.cancer.mut the above
> p.cancer.mutant <- tab[1,2]/(tab[1,1]+tab[1,2])
> p.nocancer.mutant <- tab[1,1]/(tab[1,1]+tab[1,2])
> p.cancer.mutant
[1] 0.1357143
> 1-p.cancer.mutant
[1] 0.8642857
> p.nocancer.mutant
[1] 0.8642857
>
> p.cancer.norm <- tab[2,2]/(tab[2,1]+tab[2,2])
> p.nocancer.norm <- 1-p.cancer.norm
> p.cancer.norm
[1] 0.02777778
> p.nocancer.norm
[1] 0.9722222
>
> odds(p.cancer.mutant)
[1] 0.1570248
> odds(p.cancer.norm)
[1] 0.02857143
> odds.ratio(p.cancer.mutant, p.cancer.norm)
[1] 5.495868
>
>
> ###########################################
>
> load("nsduh2019_adult_sub_rmph.RData")
>
> # Shorter name
> nsduh <- nsduh_adult_sub
> tab <- table(nsduh$demog_sex, nsduh$mj_lifetime)
> tab
No Yes
Male 206 260
Female 285 249
> table(nsduh$demog_sex)
Male Female
466 534
> table(nsduh$mj_lifetime)
No Yes
491 509
>
> # create functions for odd calculation
> odds <- function(p) p/(1-p)
> odds.ratio <- function(p1, p2) odds(p1)/odds(p2)
> # log odds를 구하는 function
> logit <- function(p) log(p/(1-p))
> # probability를 구하는 function
> # log(p/(1-p)) = x
> # p/(1-p) = e^x
> # p = e^x * (1-p)
> # p = e^x - e^x*p
> # p + p*e^x = e^x
> # p*(1+e^x) = e^x
> # p = e^x / (1 +e^x)
> # 따라서 아래는 p를 구하는 평션
> ilogit <- function(x) exp(x)/(1+exp(x))
> # exp() is the exponential function
>
>
> # 위의 펑션으로 구한 값
> PM.yes <- tab[1,2]/(tab[1,2]+tab[1,1])
> PF.yes <- tab[2,2]/(tab[2,2]+tab[2,1])
> OM <- odds(PM.yes)
> OF <- odds(PF.yes)
> OR.mvsf <- odds.ratio(PM.yes, PF.yes)
> # outputs
> round(c(PM.yes, PF.yes, OM, OF, OR.mvsf), 5)
[1] 0.55794 0.46629 1.26214 0.87368 1.44461
>
>
> library(dplyr)
다음의 패키지를 부착합니다: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
> # female을 reference로 바꾸는 작업 relevel
> nsduh <- nsduh %>%
+ mutate(demog_sex = relevel(demog_sex, ref = "Female"))
>
> # generalized linear modeling (glm) for logistic regression
> # note: IV가 종류측정 변인
> fit.ex6.2 <- glm(mj_lifetime ~ demog_sex,
+ family = binomial, data = nsduh)
> summary(fit.ex6.2)
Call:
glm(formula = mj_lifetime ~ demog_sex, family = binomial, data = nsduh)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.13504 0.08675 -1.557 0.11954
demog_sexMale 0.36784 0.12738 2.888 0.00388 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1386.0 on 999 degrees of freedom
Residual deviance: 1377.6 on 998 degrees of freedom
AIC: 1381.6
Number of Fisher Scoring iterations: 3
>
> # 아래는 다른 펑션으로 아웃풋이 다름름
> # library(lessR)
> # fit.ex6.2b <- Logit(mj_lifetime ~ demog_sex, data=nsduh)
>
> # 절편 해석
> # log(p) = a + b * demog_sex 에서x
> # demog_sexFemale 일 때, b = 0 이므로
> # log(p) = a = -0.13504
> a = -0.13504
> p.femal.yes <- exp(a)/(1 + exp(a)) # 위의 절편에 대한 p 값 계산
> p.femal.yes
[1] 0.4662912
> PF.yes
[1] 0.4662921
>
> # x = 1일 때, log(x) = a + b
> summary(fit.ex6.2)$coefficient # coefficient 값들 출력
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1350363 0.08674572 -1.556691 0.119543858
demog_sexMale 0.3678417 0.12737811 2.887794 0.003879538
> a <- summary(fit.ex6.2)$coefficient[1,1] # 절편 값 출력
> b <- summary(fit.ex6.2)$coefficient[2,1] # 절편 값 출력
> a
[1] -0.1350363
> b
[1] 0.3678417
> adash <- a + b
> adash
[1] 0.2328055
>
> # or
> exp(adash) / (1 + exp(adash))
[1] 0.5579399
> ilogit(adash)
[1] 0.5579399
> # 위 값은 PM.yes 값과 같다
> PM.yes
[1] 0.5579399
>
> # 한편 위에서의 b, coefficient 값은 b 인데
> # see http://commres.net/wiki/logistic_regression#odds_ratio_in_logistic
> # 아래는 odds ratio 의 값이 된다. (위의 위키문서 읽을 것)
> # 즉, x가 한 unit 증가할 때의 odds값의 ratio를 말한다.
> exp(b)
[1] 1.444613
> # [1] 1.444613
> # X변인의 한 unit이 증가함은 female에서 male이 되는 것
> # 그 때의 odds ratio는 exp(b), 1.444613
> # female에 비해서 male의 me가 약 1.45배 정도 된다는 것것
>
>
> # 위의 b값에 대한 CI을 구하기 위해서 confint 펑션을 사용한다
> confint(fit.ex6.2)
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
(Intercept) -0.3054835 0.03475989
demog_sexMale 0.1185985 0.61808060
> # b값에 대한 것만 알고 싶으므로
> confint(fit.ex6.2)[2, , drop=F]
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
demog_sexMale 0.1185985 0.6180806
> # 그리고 이 값들의 실제 odds ratio값을 보려면
> exp(confint(fit.ex6.2)[2,])
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
1.125918 1.855363
> # 위는 female에서 male이 될면, 즉, 0 -> 1 이 될때의
> # 마리화나 경험의 (me) odds ratio는 CI 범위는
> # 1.126 ~ 1.855 즉, 13% 에서 18% 정도 증가한다는 뜻이다다
>
> # install.packages("car")
> library(car)
필요한 패키지를 로딩중입니다: carData
다음의 패키지를 부착합니다: ‘car’
The following object is masked _by_ ‘.GlobalEnv’:
logit
The following object is masked from ‘package:dplyr’:
recode
> # coefficient probability test
> car::Anova(fit.ex6.2, type = 3, test.statistic = "Wald")
Analysis of Deviance Table (Type III tests)
Response: mj_lifetime
Df Chisq Pr(>Chisq)
(Intercept) 1 2.4233 0.11954
demog_sex 1 8.3394 0.00388 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> ########################################
> ########################################
> ########################################
> # numeric IV
> fit.ex6.3 <- glm(mj_lifetime ~ alc_agefirst,
+ family = binomial, data = nsduh)
> summary(fit.ex6.3)
Call:
glm(formula = mj_lifetime ~ alc_agefirst, family = binomial,
data = nsduh)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.3407 0.4747 11.25 <2e-16 ***
alc_agefirst -0.2835 0.0267 -10.62 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1141.45 on 842 degrees of freedom
Residual deviance: 968.44 on 841 degrees of freedom
(결측으로 인하여 157개의 관측치가 삭제되었습니다.)
AIC: 972.44
Number of Fisher Scoring iterations: 5
> round(summary(fit.ex6.3)$coef, 4)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.3407 0.4747 11.2510 0
alc_agefirst -0.2835 0.0267 -10.6181 0
>
>
> ilogit(coef(fit.ex6.3)["(Intercept)"])
(Intercept)
0.9952302
> or <- exp(coef(fit.ex6.3)["alc_agefirst"])
> 1-or
alc_agefirst
0.2468802
>
> # 0.9952 = prob of starting marijuana
> # when the age is zero (intercept이므로)
>
> # age = 0 에서 추정하는 것은 이상함
> summary(nsduh$alc_agefirst)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
3.00 15.00 17.00 17.49 19.00 45.00 157
>
> # 위의 아웃풋에서 Mean값이 약 17이므로 17을
> # 기준으로 하여 다시 보면
>
> nsduh <- nsduh %>%
+ mutate(calc_agefirst = alc_agefirst - 17)
> fit.ex6.3.centered <- glm(mj_lifetime ~ calc_agefirst,
+ family = binomial, data = nsduh)
> summary(fit.ex6.3.centered)
Call:
glm(formula = mj_lifetime ~ calc_agefirst, family = binomial,
data = nsduh)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.52065 0.07898 6.592 4.34e-11 ***
calc_agefirst -0.28353 0.02670 -10.618 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1141.45 on 842 degrees of freedom
Residual deviance: 968.44 on 841 degrees of freedom
(결측으로 인하여 157개의 관측치가 삭제되었습니다.)
AIC: 972.44
Number of Fisher Scoring iterations: 5
> # 우선 나이가 17살일 때 (x=0)
> # 절편 값인 0.5207이 logit 값이 된다.
> # 이 때, prob 를 구하면 아래와 같다
> p <- ilogit(coef(fit.ex6.3.centered)["(Intercept)"])
> p
(Intercept)
0.6273001
> # 아니면
> a <- fit.ex6.3.centered$coefficients["(Intercept)"]
> exp(a)/(1+exp(a))
(Intercept)
0.6273001
>
> # 독립변인인 나이에 대한 해석
> # b coefficient
> # 17살일 때를 기준으로 한살씩 증가할 때마다의
> # 마리화나 경험/비경험의 Odds ratio는 -0.2835
> # 이를 수치화하면 (odds ratio는 exp(b)이다)
> exp(coef(fit.ex6.3.centered)["calc_agefirst"])
calc_agefirst
0.7531198
>
> # 이는 17이후에 한살씩 알콜처음 경험을 늦추면
> # 마리화나 경험율 대 미경험 odds ratio가
> # 0.247 낮아진다고 할 수 있다 (0.7531 증가는)
> # 1-0.7531로 보는 것
>
> # 그리고 이에 대한 CI를 보면 아래와 같고
> confint(fit.ex6.3.centered)[2, , drop = F]
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
calc_agefirst -0.3375819 -0.232863
> # 이를 승비로 (odds ratio) 보면
> exp(confint(fit.ex6.3.centered)[2, , drop = F])
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
calc_agefirst 0.7134935 0.7922621
>
> #################################
> #################################
> # 1년이 아니라 3년일 경우
> fit.ex6.3.centered
Call: glm(formula = mj_lifetime ~ calc_agefirst, family = binomial,
data = nsduh)
Coefficients:
(Intercept) calc_agefirst
0.5207 -0.2835
Degrees of Freedom: 842 Total (i.e. Null); 841 Residual
(결측으로 인하여 157개의 관측치가 삭제되었습니다.)
Null Deviance: 1141
Residual Deviance: 968.4 AIC: 972.4
> coef(fit.ex6.3.centered)["calc_agefirst"]
calc_agefirst
-0.283531
> coef(fit.ex6.3.centered)["calc_agefirst"]*3
calc_agefirst
-0.850593
> exp(coef(fit.ex6.3.centered)["calc_agefirst"]*3)
calc_agefirst
0.4271616
>
> # CI 의 경우 아래와 같고
> confint(fit.ex6.3.centered)[2,]*3
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
-1.012746 -0.698589
> # 이에 해당하는 값은
> exp(confint(fit.ex6.3.centered)[2,]*3)
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
0.3632203 0.4972865
>
>
> #################################
> #################################
> ## Multiple regression
> #################################
> fit.ex6.3.adj <- glm(mj_lifetime ~ alc_agefirst +
+ demog_age_cat6 + demog_sex +
+ demog_income,
+ family = binomial, data = nsduh)
> # Regression coefficient table
> round(summary(fit.ex6.3.adj)$coef, 4)
Estimate Std. Error z value
(Intercept) 6.2542 0.5914 10.5759
alc_agefirst -0.2754 0.0276 -9.9922
demog_age_cat626-34 -0.2962 0.3286 -0.9012
demog_age_cat635-49 -0.8043 0.2966 -2.7120
demog_age_cat650-64 -0.6899 0.2985 -2.3109
demog_age_cat665+ -1.2748 0.3043 -4.1893
demog_sexMale -0.0609 0.1618 -0.3763
demog_income$20,000 - $49,999 -0.5309 0.2664 -1.9927
demog_income$50,000 - $74,999 -0.0793 0.3049 -0.2601
demog_income$75,000 or more -0.3612 0.2532 -1.4264
Pr(>|z|)
(Intercept) 0.0000
alc_agefirst 0.0000
demog_age_cat626-34 0.3675
demog_age_cat635-49 0.0067
demog_age_cat650-64 0.0208
demog_age_cat665+ 0.0000
demog_sexMale 0.7067
demog_income$20,000 - $49,999 0.0463
demog_income$50,000 - $74,999 0.7948
demog_income$75,000 or more 0.1538
> coef.values <- coef(fit.ex6.3.adj)
> data.frame(coef.values)
coef.values
(Intercept) 6.25417324
alc_agefirst -0.27541454
demog_age_cat626-34 -0.29618703
demog_age_cat635-49 -0.80427437
demog_age_cat650-64 -0.68990572
demog_age_cat665+ -1.27475385
demog_sexMale -0.06088993
demog_income$20,000 - $49,999 -0.53087558
demog_income$50,000 - $74,999 -0.07930897
demog_income$75,000 or more -0.36119745
>
> logit.values <- exp(coef(fit.ex6.3.adj))
> data.frame(logit.values)
logit.values
(Intercept) 520.1791313
alc_agefirst 0.7592573
demog_age_cat626-34 0.7436483
demog_age_cat635-49 0.4474125
demog_age_cat650-64 0.5016234
demog_age_cat665+ 0.2794998
demog_sexMale 0.9409268
demog_income$20,000 - $49,999 0.5880898
demog_income$50,000 - $74,999 0.9237545
demog_income$75,000 or more 0.6968414
>
> confint(fit.ex6.3.adj)
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
(Intercept) 5.1309585 7.45103372
alc_agefirst -0.3312435 -0.22314999
demog_age_cat626-34 -0.9490643 0.34307731
demog_age_cat635-49 -1.4002915 -0.23429671
demog_age_cat650-64 -1.2893673 -0.11559836
demog_age_cat665+ -1.8854986 -0.68944515
demog_sexMale -0.3790935 0.25566208
demog_income$20,000 - $49,999 -1.0591496 -0.01305382
demog_income$50,000 - $74,999 -0.6785210 0.51882750
demog_income$75,000 or more -0.8643471 0.13016735
> exp(confint(fit.ex6.3.adj))
프로파일링이 완료되길 기다리는 중입니다...
2.5 % 97.5 %
(Intercept) 169.1792047 1721.6419228
alc_agefirst 0.7180303 0.7999949
demog_age_cat626-34 0.3871031 1.4092777
demog_age_cat635-49 0.2465251 0.7911270
demog_age_cat650-64 0.2754450 0.8908329
demog_age_cat665+ 0.1517534 0.5018544
demog_sexMale 0.6844816 1.2913163
demog_income$20,000 - $49,999 0.3467506 0.9870310
demog_income$50,000 - $74,999 0.5073668 1.6800566
demog_income$75,000 or more 0.4213265 1.1390190
> confint.logit.values <- exp(confint(fit.ex6.3.adj))
프로파일링이 완료되길 기다리는 중입니다...
>
> # logit.values 를 Adjusted Odds Ratio로 하고 출력
> cbind("AdjOR" = logit.values, confint.logit.values)
AdjOR 2.5 %
(Intercept) 520.1791313 169.1792047
alc_agefirst 0.7592573 0.7180303
demog_age_cat626-34 0.7436483 0.3871031
demog_age_cat635-49 0.4474125 0.2465251
demog_age_cat650-64 0.5016234 0.2754450
demog_age_cat665+ 0.2794998 0.1517534
demog_sexMale 0.9409268 0.6844816
demog_income$20,000 - $49,999 0.5880898 0.3467506
demog_income$50,000 - $74,999 0.9237545 0.5073668
demog_income$75,000 or more 0.6968414 0.4213265
97.5 %
(Intercept) 1721.6419228
alc_agefirst 0.7999949
demog_age_cat626-34 1.4092777
demog_age_cat635-49 0.7911270
demog_age_cat650-64 0.8908329
demog_age_cat665+ 0.5018544
demog_sexMale 1.2913163
demog_income$20,000 - $49,999 0.9870310
demog_income$50,000 - $74,999 1.6800566
demog_income$75,000 or more 1.1390190
> # 절편은 coefficient 해석과 (odds ratio) 관계 없으므로 생략
> cbind("AdjOR" = logit.values, confint.logit.values)[-1,]
AdjOR 2.5 % 97.5 %
alc_agefirst 0.7592573 0.7180303 0.7999949
demog_age_cat626-34 0.7436483 0.3871031 1.4092777
demog_age_cat635-49 0.4474125 0.2465251 0.7911270
demog_age_cat650-64 0.5016234 0.2754450 0.8908329
demog_age_cat665+ 0.2794998 0.1517534 0.5018544
demog_sexMale 0.9409268 0.6844816 1.2913163
demog_income$20,000 - $49,999 0.5880898 0.3467506 0.9870310
demog_income$50,000 - $74,999 0.9237545 0.5073668 1.6800566
demog_income$75,000 or more 0.6968414 0.4213265 1.1390190
> # 소수점 3으로 제약
> round(cbind("AdjOR" = logit.values,
+ confint.logit.values)[-1,],3)
AdjOR 2.5 % 97.5 %
alc_agefirst 0.759 0.718 0.800
demog_age_cat626-34 0.744 0.387 1.409
demog_age_cat635-49 0.447 0.247 0.791
demog_age_cat650-64 0.502 0.275 0.891
demog_age_cat665+ 0.279 0.152 0.502
demog_sexMale 0.941 0.684 1.291
demog_income$20,000 - $49,999 0.588 0.347 0.987
demog_income$50,000 - $74,999 0.924 0.507 1.680
demog_income$75,000 or more 0.697 0.421 1.139
>
> # OR은 coefficient 값을 이야기하는 것을 다시 확인
>
> # 또한 wald significant test
> # P-values
> car::Anova(fit.ex6.3.adj, type = 3, test.statistic = "Wald")
Analysis of Deviance Table (Type III tests)
Response: mj_lifetime
Df Chisq Pr(>Chisq)
(Intercept) 1 111.8504 < 2.2e-16 ***
alc_agefirst 1 99.8435 < 2.2e-16 ***
demog_age_cat6 4 23.0107 0.000126 ***
demog_sex 1 0.1416 0.706685
demog_income 3 5.4449 0.141974
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
경고메시지(들):
1: printHypothesis(L, rhs, names(b))에서:
one or more coefficients in the hypothesis include
arithmetic operators in their names;
the printed representation of the hypothesis will be omitted
2: printHypothesis(L, rhs, names(b))에서:
one or more coefficients in the hypothesis include
arithmetic operators in their names;
the printed representation of the hypothesis will be omitted
>