로지스틱 회귀분석에 대해 질문 올립니다.

SOL 2017.09.15 13:10 Views : 13917

안녕하세요. 교수님.

로지스틱 회귀분석 중 교수님께 여쭙고 싶은 내용이 있어서 질문을 올립니다.

Result of Multiple Logistic Regression
	Estimate	Std. Error	z value	Pr(>\|z\|)	OR	lcl	ucl
(Intercept)	-3.1239	0.3416	-9.15	0.0000	0.04	0.02	0.08
Age	0.0547	0.0112	4.88	0.0000	1.06	1.03	1.08
AT_PN	-1.3352	0.7450	-1.79	0.0731	0.26	0.05	0.98
glm( Recalcitrant ~Age+AT_PN ,data=df,family=binomial)

여기서 종속 변수는 0 또는 1로 코딩되어있으며,

Age 는 연속형 변수, AT_PN 은 0 또는 1로 코딩되어 있는 명목형 변수입니다.

표에서 95% CI 는 0.05-0.98 로 구간이 1 을 포함하지 않고 있는데 Pr(>|z|) value 는 0.073 으로 0.05 보다 큰 것 같습니다.

이 변수 뿐만 아니라 몇몇 변수에서 동일한 결과를 보이고 있는데,

혹시 이 경우에는 결과를 어떻게 해석을 해야 하는지 궁금합니다.

Pr(>|z|) cut-off 를 0.05 가 아닌 다른 수치로 봐야하는 것인지.. 혹은 lcl-ucl 이 95% CI 가 아닌 것인지 궁금합니다.

추가적으로 한가지만 더 답변해주시면 정말 감사드리겠습니다.

위 분석은 "2차 치료까지 필요할 정도로 질환이 recalcitrant 하였는지 여부"를 종속 변수로 둔 것인데,

AT_PN 이라는 finding 여부의 영향을 age 를 covariate 로 하여 adjusting 하고자 합니다.

그렇다면 위 방법처럼 로지스틱 회귀 분석에 AT_PN 및 age 두 변수만 놓고 돌리면 틀린 것은 없는 것인지요?

현재 Stepwise regression 옵션은 풀고 multiple regression 투입 p value 을 1 로 설정하여 두 변수에 의한 회귀 결과가 강제로 나오도록 하여 분석하고 있습니다.

항상 감사드립니다.

I want to

List

Comment 12

스튜던트

2017.09.15 23:10

1.AT_PN 의 회귀계수(-1.3352)와 표준오차(0.7450) 이 맞게 구해졌다는 가정하에, 오즈비의 95% 신뢰구간을 구해보면

exp(-1.3352-1.96*0.7450), exp(-1.3352+1.96*0.7450) -> (0.06,1.13) 으로 나옵니다. 오즈비가 1을 포함하므로

유의한 설명변수가 아니라며, 이는 p-value가 0.05보다 큰 것과 일치 합니다.

아마도 올리주신 표에 오즈비의 신뢰구간이 잘못 기재되었거나, 잘못 계산된것으로 보입니다

2. 네. 위의 회귀식을 해석하면, '나이가 일정할때, AP-PN이 1일때 recalcitrant가 일어날 오즈는 AP-PN이 0일때 비하여 0.26배다'가 됩니다

Comment
cardiomoon

2017.09.15 23:38

웹R에서 OR 계산하는 부분을 다시 검토해보고 글 올리겠습니다.

Comment
천감기

2017.09.15 23:42

1. 번에 대해서는 원자료로 한번 분석해보고 싶어지네요..일반적인 결과가 아니므로 착오 등이 있지는 않은지 궁금합니다.

2. 에 대해서 일반 위의논지는 맞습니다 age를 보정한 APPN의 여향을 보려고 하는.....

그런데 다른 질문이 생기죠?

왜. 성별 다른 변수는 아닌 나이만 분석했느냐?

나이만이 중요하냐...

그러므로 일반적으로 여러 변수를 투입해서 단변수분석부터 하면서 변수들을 제거해가는 방식으로 하죠.

cutㅡoff를 1로했다는 것은 stepwise 등 하지 않고 그냥 enter했다는 뜻인것 같네요. 그렇게 하기도하고 특히 변수가 작을 때 일 때 조금더 유리하죠.

Comment

로지스틱 회귀분석에서 회귀계수의 신뢰구간

로지스틱 회귀분석에서 신뢰구간을 구하는 것은 다음에 따릅니다. 원 자료가 없으므로 survival 패키지의 colon 데이터를 예로 들겠습니다.

require(survival)

Loading required package: survival

colon<-na.omit(colon)
result=glm(status~rx+sex+age+obstruct+perfor,family=binomial,data=colon)
summary(result)


Call:
glm(formula = status ~ rx + sex + age + obstruct + perfor, family = binomial, 
    data = colon)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4485  -1.1949  -0.9553   1.1368   1.4223  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.346152   0.262470   1.319    0.187    
rxLev       -0.091277   0.116251  -0.785    0.432    
rxLev+5FU   -0.619551   0.118099  -5.246 1.55e-07 ***
sex         -0.102240   0.096405  -1.061    0.289    
age         -0.002325   0.004052  -0.574    0.566    
obstruct     0.197600   0.122820   1.609    0.108    
perfor       0.367541   0.284735   1.291    0.197    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2461.7  on 1775  degrees of freedom
Residual deviance: 2423.3  on 1769  degrees of freedom
AIC: 2437.3

Number of Fisher Scoring iterations: 4

오즈비 및 신뢰구간

오즈비 및 신뢰구간은 다음과 같이 구할수 있습니다.

# Odds ratio
exp(coef(result))

(Intercept)       rxLev   rxLev+5FU         sex         age    obstruct 
  1.4136174   0.9127645   0.5381858   0.9028127   0.9976777   1.2184750 
     perfor 
  1.4441789

# confidence Interval 
exp(confint(result))

Waiting for profiling to be done...

                2.5 %    97.5 %
(Intercept) 0.8455328 2.3671160
rxLev       0.7266781 1.1463109
rxLev+5FU   0.4266480 0.6779222
sex         0.7472524 1.0904893
age         0.9897779 1.0056315
obstruct    0.9580619 1.5510076
perfor      0.8311164 2.5534969

웹 R에서 odds ratio의 계산

웹R에서 odds raio의 계산은 위의 공식을 따릅니다. 다음과 같은 함수를 만들어 계산합니다.

require(moonBook)

Loading required package: moonBook

extractOR

function (x, digits = 2) 
{
    suppressMessages(a <- confint(x))
    result = data.frame(exp(coef(x)), exp(a))
    result = round(result, digits)
    result = cbind(result, round(summary(x)$coefficient[, 4], 
        4))
    colnames(result) = c("OR", "lcl", "ucl", "p")
    result
}
<environment: namespace:moonBook>

extractOR(result)

              OR  lcl  ucl      p
(Intercept) 1.41 0.85 2.37 0.1872
rxLev       0.91 0.73 1.15 0.4323
rxLev+5FU   0.54 0.43 0.68 0.0000
sex         0.90 0.75 1.09 0.2889
age         1.00 0.99 1.01 0.5661
obstruct    1.22 0.96 1.55 0.1076
perfor      1.44 0.83 2.55 0.1968

회귀계수로부터 manual 로 계산하기(1)

회귀모형 중 마지막 설명변수인 perfor에 대해 손으로 계산해보겠습니다.

(res=0.3675+c(-1,1)*1.96*0.2847)

[1] -0.190512  0.925512

exp(res)

[1] 0.8265358 2.5231598

손계산이 정확한가?

여기서 회귀계수인 0.3675에서 표준에러의 1.96배를 빼준 것인데 1.96이라는 숫자는 모형의 선형가정이 맞다고 볼때 자료의 갯수가 많아지면 t분포에 따르고 1.96에 수렴하기 때문입니다. sample size가 작을 경우 1.96보다 커집니다. 정확하게는 장도에 따라 달라지는데 이 경우 자유도가 1769이므로 다음과 같이 계산됩니다.

qt(0.975,1769)

[1] 1.961306

즉 1.96보다는 약간 커집니다.

또한 회귀계수를 소수점 4자리에서 반올림하여 계산한 것과 싵제 값은 차이가 있습니다.

회귀계수로부터 manual 로 계산하기(2)

따라서 회귀계수의 계산을 반올림 없이 하면 다음과 같이 됩니다.

res2=summary(result)$coef[7,1]+c(-1,1)*qt(0.975,1769)*summary(result)$coef[7,2]
res2

[1] -0.1909117  0.9259936

exp(res2)

[1] 0.8262055 2.5243751

역시 confint함수로 계산한 것과는 0.03 정도 차이가 나게 됩니다.

샘플 수가 적은 경우(자유도가 적은 경우) confint()함수로 계산한 것과 손으로 반올림해서 계산한 것은 차이가 날수 있을 것 같습니다. 이경우 자유도가 얼마죠?

자료를 올려봐주시면 검토하는데 도움이 되겠습니다. 스튜던트님 천감기님 같이 검토 부탁드립니다.

Comment

SOL

2017.09.16 14:04

파일 업로드 하였습니다. 일부 민감한 정보를 제외한 raw data 를 업로드 하였습니다.

Recalcitrant ~ (Age + ) Type, Recalcitrant ~ (Age + ) AT_PN 에서 동일한 이슈가 발생하였습니다.

많은 선생님들께 감사드립니다.
- Coding_for_upload.csv [File Size:12.6KB/Download:63]
Comment
스튜던트

2017.09.17 00:38

SOL님이 올려주신 95% 신뢰구간을 가지고 역으로 추정해보니,

AT_PN 모델에서는 표준오차에 1.96을 곱한게 아니라 1.76을 곱하였네요..

원래자료는 자유도가 약 270개 정도 되기 때문에

n수가 부족해서 나타나는 현상은 아닌 것 같고,

아마 문교수님이 말씀하신 것 처럼 비선형 회귀분석이기 때문에 발생하는

현상 같아요.

이건 좀 더 찾아봐야 겠습니다. 알게되면 다시 글 올리겠습니다

Comment

일단 올려주신 자료를 R에서 로지스틱회귀분석을 해 보았습니다.

data=read.csv("Downloads/Coding_for_upload.csv",stringsAsFactors = FALSE)
fit=glm(Recalcitrant~Age+AT_PN,family=binomial,data=data)
summary(fit)


Call:
glm(formula = Recalcitrant ~ Age + AT_PN, family = binomial, 
    data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6841  -0.4841  -0.3443  -0.3014   2.7659  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.12392    0.34158  -9.145  < 2e-16 ***
Age          0.05467    0.01120   4.882 1.05e-06 ***
AT_PN       -1.33523    0.74500  -1.792   0.0731 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 199.26  on 280  degrees of freedom
Residual deviance: 173.00  on 278  degrees of freedom
AIC: 179

Number of Fisher Scoring iterations: 5

회귀계수의 신뢰구간

confint(fit)

Waiting for profiling to be done...

                  2.5 %      97.5 %
(Intercept) -3.84284465 -2.49633775
Age          0.03342183  0.07768483
AT_PN       -2.99254309 -0.01690978

오즈비

confint(fit)

Waiting for profiling to be done...

                  2.5 %      97.5 %
(Intercept) -3.84284465 -2.49633775
Age          0.03342183  0.07768483
AT_PN       -2.99254309 -0.01690978

exp(confint(fit))

Waiting for profiling to be done...

                 2.5 %     97.5 %
(Intercept) 0.02143255 0.08238617
Age         1.03398662 1.08078198
AT_PN       0.05015971 0.98323239

웹 R의 방법

require(moonBook)

Loading required package: moonBook

extractOR(fit)

              OR  lcl  ucl      p
(Intercept) 0.04 0.02 0.08 0.0000
Age         1.06 1.03 1.08 0.0000
AT_PN       0.26 0.05 0.98 0.0731

이상 웹R의 방법이 문제가 없는 것은 확인하였습니다.

여기서 문제는 회귀계수의 신뢰구간이 손으로 계산한것과 다르다는 것입니다. 왜 그럴까요? AT_PN의 경우 회귀계수 -1.33523 이고 표준오차가 0.74500 이므로 다음과 같습니다.
fit$coef[3]+c(-1,1)*qt(0.975,df=278)*summary(fit)$coef[3,2]
[1] -2.801793  0.131329
이 값은 다음과 차이가 많이 납니다. 
confint(fit)
Waiting for profiling to be done...
                  2.5 %      97.5 %
(Intercept) -3.84284465 -2.49633775
Age          0.03342183  0.07768483
AT_PN       -2.99254309 -0.01690978

이 문제는 조금 더 찾아보고 다시 글 올리겠습니다.

Comment

스튜던트

2017.09.17 19:01

답을 찾은 것 같네요. 저도 교수님 처럼 붙여넣기 하고 싶은데

수식이있다 보니 깨져서 나옵니다..

http://rpubs.com/KOSUNGWOO/308829

여기에 가서 보시거나, 첨부 파일을 읽어보시면 됩니다

parameter.html
- parameter.html [File Size:699.4KB/Download:530]
Comment
cardiomoon

2017.09.17 20:56

스튜던트님 잘 보았습니다. 그러면 가능도비로 p값을 구하려면 어떻게 해야 하나요?

Comment
스튜던트

2017.09.17 22:31

model2<-glm(Recalcitrant~Age,data = sov,family = binomial)

# model1과 비교하기위해 AT_PN변수를 뺀 model2를 적합하였음

library(lmtest)
lrtest(model2,model1)

이렇게 구하거나,

anova(model2,model1,test = "LRT")

이렇게 하면 AT_PN 계수에 대한 가능도비 검정의 p-value를 구할 수 있습니다.

만족할만한 대답이 되었는지는 모르겠습니다만 ㅠㅠ

Comment

                Variables in the Equation
                                95.0% C.I.for EXP(B)
       B    S.E.    Wald    df    Sig.    Exp(B)    Lower    Upper
Step 1a Age    .055    .011    23.829    1    .000    1.056    1.033    1.080
   AT_PN    -1.335    .745    3.212    1    .073    .263    .061    1.133
Constant    -3.124    .342    83.636    1    .000    .044
a. Variable(s) entered on step 1: Age, AT_PN.

요건 SPSS에서 구한 것입니다.

일단 R에서 구한 것과는 다른 값을 제시하고 있고

엑셀과 SPSS는 같은 값을 보여 주네요.

Comment

SOL

2017.09.20 11:26

제가 모르는 별천지의 이야기라서 이해를 하지는 못하였지만 많이 찾아봐주시고 도움 주신 선생님들께 감사드립니다.

Comment

No.	Subject	Author	Date
1214	안녕하십니까 교수님 PSM 결과에 대해 궁금한 것이 있어 여쭤봅니다 [2]	김도훈	2023.06.29
1213	교수님 지난 번 문의 드렸던 PSM 재문의 드립니다. [3]	긍정적인마음가짐	2023.06.27
1212	안녕하세요 ^^ odds ratio 와 관련된 문의드리고자 합니다 ^^ [2]	성열	2023.06.26
1211	mytable() 함수에서 na처리와 pvalue 계산 문의드립니다. [1]	hjk34	2023.06.23
1210	교수님 안녕하세요 ^^ web-R table 작성관련 문의드립니다. [2]	mkw7654	2023.06.21
1209	Methods for continuous variables과 P value corrected for multiple comparisons 방법 궁금합니다. [1]	호호셀	2023.06.19
1208	ordinal logistic regression analysis 문의 [1]	버블	2023.06.18
1207	e-mail 변경 문의 [1]	kimdo	2023.06.16
1206	안녕하세요 number at risk 관련 문의 드립니다 [1]	cebukim	2023.06.13
1205	Propensity matching score 에 대해 문의드립니다. [1]	saga	2023.06.08
1204	콕스비례위험모델의 적용에 관련해서 질문 드립니다. [2]	kng1032	2023.06.07
1203	기초통계 계산 [2]	leesav	2023.06.07
1202	정회원 후원하였습니다. [2]	toast	2023.06.04
1201	logistic regression analysis 가 에러가 뜹니다. [1]	Ikarus86	2023.06.02
1200	PMS caliper 에 대해 여쭤봅니다. [2]	pjhaaa	2023.05.30
1199	ggGam에 대하여 [2]	leesav	2023.05.22
1198	R 에서 PFS 의 confidence interval 을 구하고 싶을 때 어떻게 하면 될까요? [1]	ddeeny	2023.05.21
1197	메타분석 관련 문의드립니다	화이팅	2023.05.12
1196	교수님, 새로운 앱 개발 문의드립니다! [4]	swpapa	2023.05.08
1195	교수님, 안녕하세요? 전용 앱 관련 문의드립니다! [1]	swpapa	2023.05.07

Write Tag

First Page 1 2 3 4 5 6 7 8 9 10 Last Page

	coeff b	s.e.	Wald	p-value	exp(b)	lower	upper
Intercept	-3.12392	0.341588	83.63607	5.95E-20	0.043984
AT_PN	-1.33523	0.745038	3.211856	0.073107	0.263097	0.061086	1.133167
Age	0.054675	0.011201	23.82864	1.05E-06	1.056197	1.033264	1.07964