COMMunication
RESearch.NET

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== R code ======

<code>
# 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")
</code>
====== output ======

<code>
> # 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
> 
</code>