COMMunication
RESearch.NET

This is an old revision of the document!

---

j <- log10(10)
j
10^j
k <- log2(10)
k
2^k

set.seed(101)
n <- 200
p.cancer <- 0.08
p.mutant <- 0.39

c <- runif(n, 0, 1)
canc <- ifelse(c>=p.cancer, "nocancer", "cancer")
c <- runif(n, 0, 1)
gene <- ifelse(c>=p.mutant, "norm", "mutated")

da <- data.frame(gene, canc)
da
tab <- table(da)
tab

odd.canc.mutated <- 10/119
odd.canc.norm <- 23/198
odd.ratio <- odd.canc.mutated/odd.canc.norm
odd.ratio
log(odd.ratio)

##########
n.mut <- 23+117
n.norm <- 6+210

p.cancer.mut <- 23/(23+117)
p.cancer.norm <- 6/(6+210)
set.seed(101)
c <- runif(n.mut, 0, 1)
canc.mut <- ifelse(c>=p.cancer.mut, "nocancer", "cancer")
c <- runif(n.norm, 0, 1)
canc.norm <- ifelse(c>=p.cancer.norm, "nocancer", "cancer")
table(canc.mut)
table(canc.norm)


load("nsduh2019_adult_sub_rmph.RData")

# Shorter name
nsduh <- nsduh_adult_sub
tab <- table(nsduh$demog_sex, nsduh$mj_lifetime)
tab


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 
# p = e^x/(1+e^x)
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)
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"))

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)

# 절편 해석 
# demog_sexFemale 일 때, 
exp(-0.13504)/(1 + exp(-0.13504)) # 위의 절편에 대한 p 값 계산 
PF.yes
# p = e^x/1+e^x
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))
p <- ilogit(adash) 
p 
# 이 값은 PM.yes 값과 같다
PM.yes

# coefficient 아웃풋을 이용하여 b값은 아래와 같이 얻는다
b <- summary(fit.ex6.2)$coefficient[2,1]
b
exp(1)^b

# 위의 b값에 대한 CI을 구하기 위해서 confint 펑션을 사용한다
confint(fit.ex6.2)
# b값에 대한 것만 알고 싶으므로 
confint(fit.ex6.2)[2, , drop=F]
# 그리고 이 값들의 실제 odds ratio값을 보려면 
exp(confint(fit.ex6.2)[2,])


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