Differences

This shows you the differences between two versions of the page.

--- logistic_regression:examples_r [2024/12/10 14:34] – hkimscil
+++ logistic_regression:examples_r [2024/12/12 10:27] (current) – hkimscil
@@ Line 1: / Line 1: @@
+====== R code ======
 <code>
 # log cacluation
@@ Line 266: / Line 268: @@
 # 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
+    0      0
+    0      1
+    0      0
+    0      0
+    0      0
+    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
+mutant nocancer
+mutant   cancer
+mutant nocancer
+mutant nocancer
+mutant nocancer
+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
+    534
+> table(nsduh$mj_lifetime)
+ No Yes
+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)
+프로파일링이 완료되길 기다리는 중입니다...
+.5 %     97.5 %
+(Intercept)   -0.3054835 0.03475989
+demog_sexMale  0.1185985 0.61808060
+> # b값에 대한 것만 알고 싶으므로
+> confint(fit.ex6.2)[2, , drop=F]
+프로파일링이 완료되길 기다리는 중입니다...
+.5 %    97.5 %
+demog_sexMale 0.1185985 0.6180806
+> # 그리고 이 값들의 실제 odds ratio값을 보려면
+> exp(confint(fit.ex6.2)[2,])
+프로파일링이 완료되길 기다리는 중입니다...
+.5 %   97.5 %
+.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)
+.9952302
+> or <- exp(coef(fit.ex6.3)["alc_agefirst"])
+> 1-or
+alc_agefirst
+.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
+.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)
+.6273001
+> # 아니면
+> a <- fit.ex6.3.centered$coefficients["(Intercept)"]
+> exp(a)/(1+exp(a))
+(Intercept)
+.6273001
+>
+> # 독립변인인 나이에 대한 해석
+> # b coefficient
+> # 17살일 때를 기준으로 한살씩 증가할 때마다의
+> # 마리화나 경험/비경험의 Odds ratio는  -0.2835
+> # 이를 수치화하면 (odds ratio는 exp(b)이다)
+> exp(coef(fit.ex6.3.centered)["calc_agefirst"])
+calc_agefirst
+.7531198
+>
+> # 이는 17이후에 한살씩 알콜처음 경험을 늦추면
+> # 마리화나 경험율 대 미경험 odds ratio가
+> # 0.247 낮아진다고 할 수 있다 (0.7531 증가는)
+> # 1-0.7531로 보는 것
+>
+> # 그리고 이에 대한 CI를 보면 아래와 같고
+> confint(fit.ex6.3.centered)[2, , drop = F]
+프로파일링이 완료되길 기다리는 중입니다...
+.5 %    97.5 %
+calc_agefirst -0.3375819 -0.232863
+> # 이를 승비로 (odds ratio) 보면
+> exp(confint(fit.ex6.3.centered)[2, , drop = F])
+프로파일링이 완료되길 기다리는 중입니다...
+.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
+.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
+.4271616
+>
+> # CI 의 경우 아래와 같고
+> confint(fit.ex6.3.centered)[2,]*3
+프로파일링이 완료되길 기다리는 중입니다...
+.5 %    97.5 %
+-1.012746 -0.698589
+> # 이에 해당하는 값은
+> exp(confint(fit.ex6.3.centered)[2,]*3)
+프로파일링이 완료되길 기다리는 중입니다...
+.5 %    97.5 %
+.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)
+프로파일링이 완료되길 기다리는 중입니다...
+.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))
+프로파일링이 완료되길 기다리는 중입니다...
+.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
+.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
+경고메시지(들):
+: 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
+: 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>