COMMunication
RESearch.NET

This is an old revision of the document!

---

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