#title 회귀분석 [[TableOfContents]] ==== 회귀분석 기초 ==== * '회귀(regression)' 개념은 19세기 말 영국의 생물통계학자 골튼(F.Galton)에 의해 처음 이용됨. * 평균으로의 회귀. * 변수 사이의 함수적 관계를 나타내는 수학적 회귀방정식을 구하고 독립변수를 특정한 값에 따른 종속변수의 값을 예측하는 기법 * 회귀분석은 서로 영향을 주고 받으면서 변화하는 인과관계를 갖는 두 변수 사이의 관계를 분석 * 산포도를 그려보면 회귀분석을 할 것인지 말 것인지 대략 결정할 수 있다. (양/음의 상관관계) * 현상의 정량적 이해나 Y의 예측을 위하여 유효한 회귀 모델을 구축하는 것이 회귀분석의 주목적. * 잔차(residual)분석: 모델검증 * 회귀진단: 특이한 이상치나 모수의 추정에 큰 영향을 주는 관측치 검출 * 목적변수가 양적변수일때 적용 가능 * 설명변수는 양적이거나 질적이거나 상관없으나, 질적변수는 데이터의 값이 0 또는 1만 취하는 더미변수를 도입하여 질적인 데이터를 변환해야 한다. * 종속변수(y)는 정규분포라 가정 * 회귀식과 표준편차 * y = a + bx ± 1σ : 68.0%가 이 구간내에 들어온다. * y = a + bx ± 2σ : 95.5%가 이 구간내에 들어온다. * y = a + bx ± 3σ : 99.7%가 이 구간내에 들어온다. ==== 인과 관계 성립 조건 ==== * 종속 변수와 독립 변수간에 상관 관계가 있어야 한다. (독립 변수가 변하면 종속 변수도 변해야 함) * 독립 변수가 종속 변수보다 먼저 존재해야 한다. (시간의 우선성) * 종속 변수와 독립 변수간의 상관 관계가 제3의 변수로 설명될 여지가 없어야 한다. (제3의 변수가 종속 변수를 설명해서는 안 된다) ==== 변수변환(Box-Cox변환) ==== 사용법은 다음과 같다. {{{ library("TeachingDemos") y <- rlnorm(500, 3, 2) par(mfrow=c(2,2)) qqnorm(y) qqnorm(bct(y,1/2)) qqnorm(bct(y,0)) hist(bct(y,0)) }}} attachment:회귀분석/lm02.jpg {{{ library("car") y <- rlnorm(500, 3, 2) hist(box.cox(y,0)) }}} ==== 단순회귀분석 ==== 근속년수와 연봉의 회귀분석 {{{ #x: 근속년수 x <- c(2,3,4,6,8,9,10,12,13,15,17,19,20,23,25,26,27,29,30,33) #y: 연봉 y <- c(1500,1640,1950,2700,2580,3200,3750,2960,4200,3960,4310,3990,4010,4560,4800,5300,5050,5450,5950,5650) result <- lm(y~x) summary(result) }}} 결과 {{{ Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -470.0 -290.8 -168.0 293.9 789.4 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1708.117 173.489 9.846 1.13e-08 *** x 130.960 9.089 14.409 2.52e-11 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 386.6 on 18 degrees of freedom Multiple R-squared: 0.9202, Adjusted R-squared: 0.9158 F-statistic: 207.6 on 1 and 18 DF, p-value: 2.522e-11 }}} 결과해석 * 회귀식: y = 130.960 * x + 1708.117 * 자유도: 18 * 결정계수(coefficient of determination): 0.9202(조정된 결정계수: 0.9158)는 이 모델이 92%를 설명함을 뜻함. * Pr(>|t|) : F검정의 결과 (p-value) = 2.52e-11 * 귀무가설: 회귀식에 의미가 없다. * 대립가설: 회귀식에 의미가 있다. * p-value = 2.522e-11이므로 유의수준 0.05에서 대립가설 채택할 근거 있다. 즉, 회귀식에 의미가 있다. * 잔차분석 * 더빈 왓슨비(durbin-watson ration; 오차항간의 계열상관 유무조사) * 바트렛의 검정(bartlett's test; 오차항의 분산의 균일성 조사) * 잔차: Residual standard error: 386.6 * 근속년수가 50년인 사람의 연봉은? 130.960 * 50 + 1708.117 = 8256.117 만원이다. * 회귀식에 대입해서 잔차/표준편차 = 표준잔차, -2.5 <= 표준잔차 <= 2.5 면 정상치, 범위를 벗어나면 이상치일 가능성이 높다. ==== 중회귀분석 ==== 공던지기 결과과 악력, 신장, 체중의 중회귀분석 데이터 {{{ #x1:악력, x2:신장, x3:체중, y1:공던지기 x1 <- c(28,46,39,25,34,29,38,23,42,27,35,39,38,32,25) x2 <- c(146,169,160,156,161,168,154,153,160,152,155,154,157,162,142) x3 <- c(34,57,48,38,47,50,54,40,62,39,46,54,57,53,32) y1 <- c(22,36,24,22,27,29,26,23,31,24,23,27,31,25,23) data <- data.frame(y1,x1,x2,x3) }}} 상관분석 {{{ #시각적인 관찰(산포도)이 필요하면 pairs(data)하면 된다. cov.wt(data, cor=T) > cov.wt(data, cor=T) $cov y1 x1 x2 x3 y1 16.31429 21.07143 20.01429 28.91429 x1 21.07143 48.66667 26.71429 53.64286 x2 20.01429 26.71429 52.25714 45.60000 x3 28.91429 53.64286 45.60000 82.54286 $center y1 x1 x2 x3 26.20000 33.33333 156.60000 47.40000 $n.obs [1] 15 $cor y1 x1 x2 x3 y1 1.0000000 0.7478150 0.6854618 0.7879319 x1 0.7478150 1.0000000 0.5297304 0.8463624 x2 0.6854618 0.5297304 1.0000000 0.6943082 x3 0.7879319 0.8463624 0.6943082 1.0000000 > }}} 결과해석 * x1(악력)과 x3(체중)의 상관계수가 가장 높다. 0.8463624 * x1(악력)과 x2(신장)의 상관계수가 가장 낮다. 0.5297304 * 공던지기(y1)와 상관계수가 높은 순서는 체중(x3), 악력(x1), 신장(x2)의 순이다. pairs의 결과 attachment:회귀분석/lm01.jpg?width=60% 회귀분석 {{{ summary(lm(y1~x1+x2+x3)) > summary(lm(y1~x1+x2+x3)) Call: lm(formula = y1 ~ x1 + x2 + x3) Residuals: Min 1Q Median 3Q Max -3.9976 -1.3822 0.3611 1.1549 3.9291 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -13.2173 17.6038 -0.751 0.469 x1 0.2014 0.1842 1.093 0.298 x2 0.1710 0.1316 1.300 0.220 x3 0.1249 0.1667 0.750 0.469 Residual standard error: 2.532 on 11 degrees of freedom Multiple R-squared: 0.6913, Adjusted R-squared: 0.6072 F-statistic: 8.213 on 3 and 11 DF, p-value: 0.003774 > }}} 결과해석 * 결정계수: 0.6913 (조정: 0.6072), 회귀식은 약 69%를 설명함. * p-value: 0.003774으로 유의수준 0.05에서 대립가설을 채택. 즉, 회귀식은 유의함. * 회귀식: y1 = -13.2173 + x1 * 0.2014 + x2 * 0.1710 + x3 * 0.1249 * Pr(>|t|)가 모두 0.05를 넘는 수준이다. 별로 안 좋다. * 목적변수에 대한 각 설명변수의 영향력(t) * 악력(x1) = 1.093 * 신장(x2) = 1.300 * 체중(x3) = 0.750 * 신장이 공던지기(y1)에 가장 많은 영향력을 끼치고 있다. * 편회귀계수의 유의성 판정(F) * 악력(x1) = 0.298 * 신장(x2) = 0.220 * 체중(x3) = 0.469 * 모두 0.05를 넘는다. 다중공선성을 의심해야 한다. 교호작용이 있는지도 봐야 한다. ==== 변수선택 ==== * 변수를 잘못 선택하면? * 추정의 정밀도가 나빠진다. * y의 예측치는 편의를 갖게되고, 오차분산의 추정치는 과대평가하게 된다. * 다중공선성(multicollinearity)의 문제 발생 * 변수선택방법 * 전진선택(forward selection) * 후진소거(backward elimination) * 둘다(전진선택, 후진소거) * 이는 [로지스틱 회귀분석]의 [http://databaser.net/moniwiki/wiki.php/%EB%A1%9C%EC%A7%80%EC%8A%A4%ED%8B%B1%ED%9A%8C%EA%B7%80%EB%B6%84%EC%84%9D#s-4 여기]에 사용법이 기술되어 있다. 이렇게 해도 된다. {{{ library("MASS") fit <- lm(y1~x1+x2+x3,data=mydata) step <- stepAIC(fit, direction="both") step$anova # display results }}} ==== 더미변수 ==== * 만약 성별, 요일과 같은 변수가 회귀식에 포함되어야 하면 더미변수를 추가해야 함. * p종류의 범주가 있으면 p-1개의 더미변수가 필요함. * 성별은 2가지라 0, 1로 그냥 세팅하면 되지만 요일은 다음과 같이 해야 함. * y = x + x1 + x2 + x3 + x4 + x5 + x6 {{{ x1 x2 x3 x4 x5 x6 월 0 0 0 0 0 0 화 1 0 0 0 0 0 수 0 1 0 0 0 0 목 0 0 1 0 0 0 금 0 0 0 1 0 0 토 0 0 0 0 1 0 일 0 0 0 0 0 1 }}} ==== 교호작용 ==== 가끔씩 두 개 이상의 변수 짬뽕으로 결과에 영향을 끼치는 경우가 있다. 교호작용이 없는 모델 {{{ m1 <- lm(formula=Height~Girth + Volume, data=trees) summary(m1) }}} {{{ > summary(m1) Call: lm(formula = Height ~ Girth + Volume, data = trees) Residuals: Min 1Q Median 3Q Max -9.7855 -3.3649 0.5683 2.3747 11.6910 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 83.2958 9.0866 9.167 6.33e-10 *** Girth -1.8615 1.1567 -1.609 0.1188 Volume 0.5756 0.2208 2.607 0.0145 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 5.056 on 28 degrees of freedom Multiple R-squared: 0.4123, Adjusted R-squared: 0.3703 F-statistic: 9.82 on 2 and 28 DF, p-value: 0.0005868 }}} 교호작용이 있는 모델 {{{ m2 <- lm(formula=Height~Girth + Volume + Girth:Volume, data=trees) summary(m2) }}} {{{ > summary(m2) Call: lm(formula = Height ~ Girth + Volume + Girth:Volume, data = trees) Residuals: Min 1Q Median 3Q Max -6.7781 -3.5574 -0.1512 2.3631 10.5879 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 75.40148 8.49147 8.880 1.7e-09 *** Girth -2.29632 1.03601 -2.217 0.035270 * Volume 1.86095 0.47932 3.882 0.000604 *** Girth:Volume -0.05608 0.01909 -2.938 0.006689 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.482 on 27 degrees of freedom Multiple R-squared: 0.5546, Adjusted R-squared: 0.5051 F-statistic: 11.21 on 3 and 27 DF, p-value: 5.898e-05 }}} 두 변수의 교호작용은 곱하기(Girth*trees)다. 아래의 A, B의 결과는 같다. A {{{ fitted(m2) }}} B {{{ coef(m2)[1] + coef(m2)[2]*trees$Girth + coef(m2)[3]*trees$Volume + coef(m2)[4]*trees$Girth*trees$Volume }}} ==== 다중공선성(multicollinearity) ==== * 설명변수 사이에 선형관계가 두 개 이상 있을 때 다중공선성이 있다고 말한다. * 상관분석(R함수: cov.wt())해서 상관계수가 1에 가까운 설명변수를 버린다. * 다중공선성 진단 * 설명변수 x1, x2, x3, x4가 있을 경우 * x1 ~ x2 + x3 + x4 -> 상관계수 R1 -> 허용도(tolerance) = 1 - R1^2 * x2 ~ x1 + x3 + x4 -> 상관계수 R2 -> 허용도(tolerance) = 1 - R2^2 * x3 ~ x1 + x2 + x4 -> 상관계수 R3 -> 허용도(tolerance) = 1 - R3^2 * x4 ~ x1 + x2 + x3 -> 상관계수 R4 -> 허용도(tolerance) = 1 - R4^2 * (1 - R,,j,,^2) = 허용도가 작을 때는 R,,j,,가 1에 가까우므로 x,,j,,는 다른 설명변수의 선형결합으로 표현할 수 있다. * 분산팽창계수(VIF: variance inflation factor) = 1 / (1 - R,,i,,^2)가 10보다 크면 다중공선성이 존재한다고 경험적으로 판단함. * VIP = 1 / (1 - R,,j,,^2) = 1 / (1 - 0.7432^2) = 2.233869 분산팽창계수(VIF)는 Design 라이브러리의 vif()함수를 사용한다. 다음은 그 결과이다. {{{ #install.packages("Design",dependencies=T) library("Design") result <- lm(y1~x1+x2+x3) vif(result) > library("Design") 요구된 패키지 Hmisc를 로드중입니다 다음의 패키지를 첨가합니다: 'Hmisc' The following object(s) are masked from package:base : format.pval, round.POSIXt, trunc.POSIXt, units 요구된 패키지 survival를 로드중입니다 요구된 패키지 splines를 로드중입니다 다음의 패키지를 첨가합니다: 'survival' The following object(s) are masked from package:Hmisc : untangle.specials Design library by Frank E Harrell Jr Type library(help='Design'), ?Overview, or ?Design.Overview') to see overall documentation. 다음의 패키지를 첨가합니다: 'Design' The following object(s) are masked from package:survival : Surv, survfit > result <- lm(y1~x1+x2+x3) > vif(result) x1 x2 x3 3.607546 1.975832 5.010688 > }}} 결과해석 * x1의 VIF: 3.607546 * x2의 VIF: 1.975832 * x3의 VIF: 5.010688 * VIF가 10이 넘지 않으므로 다중공선성은 없는 것으로 판단된다. ==== 회귀진단 ==== 잔차는.. * 정규분포를 따라야하고, * 분산이 일정하고, * 특별한 추세를 보이지 않아야 한다. 즉, 잔차는 정규성, 등분산성, 독립성을 만족해야 한다. {{{ par(mfrow=c(2,2)) plot(lm(formula = y1 ~ x1 + x2 + x3 + x1:x2)) }}} attachment:회귀분석/regression_diagnosis.png * 좌상: 예측값, 잔차 --> 0으로부터 고르게 분포하고 있는가? * 우상: 표준화된 잔차의 Q-Q Plot --> 정규성을 띄는가? * 좌하: 예측값, 스튜던트 잔차(표준화된 잔차)의 제곱근 --> 직선에 가까워야 한다.. 분산이 일정하지 않음을 의심해야 함. * 우하: 레버러지와 스튜던트 잔차 plot --> rs = (k + 1 / n) * 2; k = predictor의 숫자, n = sample_size; rs보다 크가면 outlier '''정규성''' {{{ > out <- lm(formula = y1 ~ x1 + x2 + x3 + x1:x2) > shapiro.test(resid(out)) Shapiro-Wilk normality test data: resid(out) W = 0.9669, p-value = 0.8099 > }}} 결과에서 볼 수 있듯이 p-value = 0.8099로 '정규분포와 차이가 없다'라는 귀무가설을 조낸 지지한다. 즉, 잔차는 정규분포를 따른다. '''등분산성''' 브로슈-파간 테스트(Breusch- Pagan Test)로 등분산인지 알아본다. {{{ > library("lmtest") > bptest(out) studentized Breusch-Pagan test data: out BP = 5.6201, df = 4, p-value = 0.2294 > }}} p-value가 0.2294로 유의수준 0.05에서 귀무가설 지지한다. 즉, 등분산이다. 만약 이분산(heterogenous variance)일 경우 가중최소제곱(weighted least squares; nls()사용)의 사용. 변수변환을 통해 교정할 수도 있다. {{{ > out <- lm(formula = y1 ~ x1 + x2 + x3 + x1:x2, weights=I(x1^(2))) > summary(out) Call: lm(formula = y1 ~ x1 + x2 + x3 + x1:x2, weights = I(x1^(2))) Residuals: Min 1Q Median 3Q Max -106.194 -34.272 4.358 30.456 85.645 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 206.98239 62.41049 3.316 0.00779 ** x1 -6.74461 1.84209 -3.661 0.00438 ** x2 -1.25865 0.40492 -3.108 0.01109 * x3 0.44112 0.13692 3.222 0.00915 ** x1:x2 0.04199 0.01105 3.800 0.00349 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 61.32 on 10 degrees of freedom Multiple R-squared: 0.8834, Adjusted R-squared: 0.8368 F-statistic: 18.94 on 4 and 10 DF, p-value: 0.0001167 }}} '''독립성''' {{{ > library("lmtest") > dwtest(out) Durbin-Watson test data: out DW = 2.3524, p-value = 0.7083 alternative hypothesis: true autocorrelation is greater than 0 > }}} * d통계량 * 0 : 양의 자기상관, p-value = 1 * 2 : 독립, p-value = 0 * 4 : 음의 자기상관, p-value = -1 * 대략적으로 DW값이 1보다 작거나 3보다 크면 자기상관이 확실히 있다고 판달할 수 있으며, 1.5~2.5 사이에 있을 경우 독립이라고 판단 * 대략 독립. 참고: * [https://www.r-bloggers.com/visualising-residuals/ Visualising Residuals] ==== 두 회귀식의 비교 ==== {{{ data(airquality) lm1<-lm(Ozone~.,airquality) # full model lm2<-lm(Ozone~Solar.R+Wind +Month+Day,airquality) # reduced model anova(lm2,lm1) }}} ==== 다변량 회귀분석 ==== 여기에서는 그냥 lm()함수의 다변량 회귀분석의 사용방법만 보자. {{{ y2 <- x1 #x1을 종속변수라 가정하고... summary(lm(cbind(y1,y2)~x2+x3)) > y2 <- x1 > summary(lm(cbind(y1,y2)~x2+x3)) Response y1 : Call: lm(formula = y1 ~ x2 + x3) Residuals: Min 1Q Median 3Q Max -3.5060 -1.5862 0.2502 0.8539 5.3777 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -9.8745 17.4765 -0.565 0.5825 x2 0.1493 0.1311 1.139 0.2770 x3 0.2678 0.1043 2.567 0.0247 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.552 on 12 degrees of freedom Multiple R-squared: 0.6578, Adjusted R-squared: 0.6008 F-statistic: 11.53 on 2 and 12 DF, p-value: 0.001605 Response y2 : Call: lm(formula = y2 ~ x2 + x3) Residuals: Min 1Q Median 3Q Max -5.4716 -1.9151 -0.2964 1.9562 7.1935 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 16.5997 27.1675 0.611 0.552589 x2 -0.1079 0.2038 -0.529 0.606186 x3 0.7095 0.1622 4.375 0.000904 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.967 on 12 degrees of freedom Multiple R-squared: 0.7228, Adjusted R-squared: 0.6766 F-statistic: 15.65 on 2 and 12 DF, p-value: 0.0004537 > }}} ==== 결과의 해석 ==== 뭘 또 이렇게까지나 해야 하는지?? 안 해도 되는거 같으다. {{{ > summary(model) Call: lm(formula = 통계 ~ 수학 + 결석, data = result) Residuals: Min 1Q Median 3Q Max -5.348 -2.274 -1.276 2.954 5.673 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 53.6832 14.1808 3.786 0.00431 ** 수학 0.6073 0.1984 3.062 0.01353 * 결석 -1.9346 0.9144 -2.116 0.06348 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.721 on 9 degrees of freedom Multiple R-squared: 0.8289, Adjusted R-squared: 0.7909 F-statistic: 21.8 on 2 and 9 DF, p-value: 0.0003544 }}} * 회귀직선이 y=통계, x1=수학, x2=결석 이라면, '''''y''''' = 53.6832 + 0.6073'''''x1''''' + -1.9346'''''x2''''' 으로 추정된다. * -1.9346'''''x2'''''으로 마이너스(-)값을 가진다는 것은 결석을 할 수록 통계성적이 나빠진다는 것을 의미한다. * 'β1=0'이라는 귀무가설에 대한 T 검정통계량은 3.062이며, 이것은 P값은 0.01353으로 유의수준 0.05에서는 귀무가설을 기각하게 된다. 즉, 유의미한 값이다. * 결정계수 R^^2^^은 아래의 anova(분산분석)의 결과에서 R^^2^^ = (541.69 + 61.97)/(541.69 + 61.97 + 124.59) = 0.8289186405이고, 조정된 R^^2^^값은 0.7909 이다. __즉, 통계성적은 수학성적과 결석횟수 2개의 변수에 의해 약 83%가 설명되어진다.__ * F 검정통계량은 21.8이며, F분포표의 α = 0.05, v1=2, v2=9의 값은 4.26이다. 즉, 21.8 > 4.26 이므로 귀무가설(β1=β2=0)은 기각한다. 즉, 유의미한 값이다. p-value가 나오므로 F분포표를 볼 필요는 없다. 유의수준 0.05 > 0.0003544 이므로 귀무가설을 기각하면 된다. * p-value가 작으면 대립가설을 지지하고, 커지면 귀무가설을 지지한다. 0.0003544는 귀무가설이 맞다고 가정했을 경우 얻어진 검정통계량보다 더 극단적인 결과가 나올 확률이다. 즉, 통계적인 수치에 근거하여 의사결정을 했을 때에 오류(1종 오류)가 나올 확률. * 유의수준 0.05라는 것은 95%만 평균에 가까운 정상으로 생각하고, 평균가 멀리 떨어진 5%는 비정상으로 생각한다는 것이다. {{{ > anova(model) #자유도, 제곱합, 제곱평균 F값, Analysis of Variance Table Response: 통계 Df Sum Sq Mean Sq F value Pr(>F) 수학 1 541.69 541.69 39.1297 0.0001487 *** 결석 1 61.97 61.97 4.4762 0.0634803 . Residuals 9 124.59 13.84 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 }}} ||요인||제곱합||자유도||제곱평균||F비|| ||회귀||603.66||11||301.83||21.802950|| ||잔차||124.59||9||13.84|| || ||계||728.25||11|| || || * F 검정통계량은 21.8이며, F분포표의 α = 0.05, v1=2, v2=9의 값은 4.26이다. 즉, 21.8 > 4.26 이므로 귀무가설(β1=β2=0)은 기각한다. 즉, 유의미한 값이다. * 결정계수 R^^2^^은 아래의 R^^2^^ = 603.66/728.25 = 0.82891864 이다. * 자유도란 자유로운 정도. 랜덤한 정도. 수식에 의해 자유가 제한되지 않는 정도인데 보통 표본의 개수를 의미한다. '''참고: p-value 뽑아내기''' {{{ fit <- lm(Volume ~ Girth + Height, data=trees) f <- summary(fit)$fstatistic pf(f[1],f[2],f[3],lower.tail=F) }}} '''참고: 95% 신뢰구간''' {{{ confint(model, level=0.95) 2.5 % 97.5 % (Intercept) 502.5676 524.66261 poly(q, 3)1 1919.2739 2232.52494 poly(q, 3)2 -264.6292 48.62188 poly(q, 3)3 707.3999 1020.65097 }}} ==== 표준화된 회귀계수(beta) ==== * B : 비표준화 회귀계수 * beta : 표준화 회귀계수 {{{ model <- lm(Girth ~ Height + Volume, data=trees) summary(model) }}} 결과 {{{ > summary(model) Call: lm(formula = Girth ~ Height + Volume, data = trees) Residuals: Min 1Q Median 3Q Max -1.34288 -0.56696 -0.08628 0.80283 1.11642 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.81637 1.97320 5.482 7.45e-06 *** Height -0.04548 0.02826 -1.609 0.119 Volume 0.19518 0.01096 17.816 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7904 on 28 degrees of freedom Multiple R-squared: 0.9408, Adjusted R-squared: 0.9366 F-statistic: 222.5 on 2 and 28 DF, p-value: < 2.2e-16 }}} 각각의 계수(B)는 나오는데, 어떤 변수가 더 큰 영향을 끼치는지 알 수 없다. 이런 경우 beta를 보면되는데, lm결과에 나오지 않는다. 다음과 같이 하면 된다. {{{ #install.packages("QuantPsyc") library("QuantPsyc") lm.beta(model) }}} 결과 {{{ > lm.beta(model) Height Volume -0.09235166 1.02236871 }}} Volume이 거의 모든 영향을 끼쳤다. 물론 beta를 보기 전에 Height 변수는 모델에서 제거 되었을 것이다. ==== predict ==== {{{ y <- seq(1: nrow(x1)) dau <- x1$dau fit <- lm(dau ~ poly(y, 3)) summary(fit) plot(dau~y) lines(y, predict(fit, data.frame(y=y)), col="purple") newdata <- data.frame(y = seq(139,145)) predict(fit, newdata) out <- lm(dau~y + I(y^2) + I(y^3)) summary(out) p <- 145 3.131e+06 - p*2.637e+04 + (p^2) * 2.438e+02 - (p^3)*1.138e+00 }}} ==== 마우스로 콕콕 찝어서 이상치 제거 후 회귀분석하기 ==== {{{ plot(y ~ x, data=df) notin <- identify(df$x, df$y, labels=row.names(df)) fit <- lm(y ~ x - 1, data=df[!rownames(df) %in% notin,]) summary(fit) plot(y ~ x, data=df[!rownames(df) %in% notin,]) abline(fit) }}} ==== coef의 신뢰구간 ==== {{{ se <- sqrt(diag(vcov(bu.model)))[2] ci.upr <- coef(bu.model)[2] + 1.96*se ci.lwr <- coef(bu.model)[2] - 1.96*se }}} ==== confidence, prediction level ==== * interval = "confidence"는 평균적인 것, 즉, 오차항(intercept)을 고려치 않아(0으로 가정) 범위가 prediction에 비해 좁다. * interval = "prediction"은 single value에 대한 것, 즉, 오차항(intercept)가 고려해야 함 {{{ # Make up some data x <- seq(0, 1, length.out = 20) y <- x^2 + rnorm(20) # fit a quadratic model <- lm(y ~ I(x^2)) fitted <- predict(model, interval = "confidence", level=.95) # plot the data and the fitted line plot(x, y) lines(x, fitted[, "fit"]) # now the confidence bands lines(x, fitted[, "lwr"], lty = "dotted") lines(x, fitted[, "upr"], lty = "dotted") }}} 참고 * http://www.sthda.com/english/articles/40-regression-analysis/166-predict-in-r-model-predictions-and-confidence-intervals/ * https://thebook.io/006723/ch08/02/03/ * https://www.real-statistics.com/regression/confidence-and-prediction-intervals/ ==== R-squared ==== R^2의 의미 * 회귀식이 설명하는 변동량 * 실제값과 예측값의 상관관계(R^2는 피어슨 상관계수의 제곱) 회귀식이 설명하는 변동량 {{{ d <- data.frame(y=(1:10)^2,x=1:10) model <- lm(y~x,data=d) d$prediction <- predict(model,newdata=d) #R-squared 1-sum((d$prediction-d$y)^2)/sum((mean(d$y)- d$y)^2) }}} 독립변수가 많아지면 R^2가 커지는 경향이 있다. 그래서 조정된 R^2를 사용한다. 여러 독립변수를 사용한 회귀식간의 비교는 조정된 R^2로 비교해야 한다. 조정된 R^2은 다음과 같이 구한다. {{{ m <- lm(Girth~Height, trees) summary(m) }}} {{{ > summary(m) Call: lm(formula = Girth ~ Height, data = trees) Residuals: Min 1Q Median 3Q Max -4.2386 -1.9205 -0.0714 2.7450 4.5384 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -6.18839 5.96020 -1.038 0.30772 Height 0.25575 0.07816 3.272 0.00276 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.728 on 29 degrees of freedom Multiple R-squared: 0.2697, Adjusted R-squared: 0.2445 F-statistic: 10.71 on 1 and 29 DF, p-value: 0.002758 }}} ==== Adjusted R-squared ==== {{{ n <- length(trees$Girth) #표본수 p <- 1 #독립변수수 R2 <- cor(trees$Girth, fitted(m))^2 #R-squared 1 - ((n-1)*(1-R2)/(n-p-1)) #Adjusted R-squared }}} ==== as.lm.nls ==== {{{ as.lm.nls <- function(object, ...) { if (!inherits(object, "nls")) { w <- paste("expected object of class nls but got object of class:", paste(class(object), collapse = " ")) warning(w) } gradient <- object$m$gradient() if (is.null(colnames(gradient))) { colnames(gradient) <- names(object$m$getPars()) } response.name <- if (length(formula(object)) == 2) "0" else as.character(formula(object)[[2]]) lhs <- object$m$lhs() L <- data.frame(lhs, gradient) names(L)[1] <- response.name fo <- sprintf("%s ~ %s - 1", response.name, paste(colnames(gradient), collapse = "+")) fo <- as.formula(fo, env = as.proto.list(L)) do.call("lm", list(fo, offset = substitute(fitted(object)))) } }}} {{{ predict(as.lm.nls(fit), data.frame(x=xx), interval="prediction", level=.99) }}} ==== 비모수 방법 ==== {{{ x <- c(3,4,5,6,7) y <- c(4.8, 6, 7.2, 7.8, 8.2) fit1 <- lm(y ~ x) summary(fit1) plot(x,y) abline(fit1) library("zyp") fit2 <- zyp.sen(y ~ x) fit2 abline(coef=coef(fit2), col="red") }}} attachment:회귀분석/np_regression.png --빨간선이 비모수 방법 ==== 오차 계산 ==== 이거 은근히 귀찮은데, 역시 만들어 놨군. {{{ y <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14) x <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69) model <- lm(y ~ x) library("DMwR") regr.eval(y, fitted(model)) plot(x,y) abline(model) }}} {{{ > regr.eval(y, fitted(model)) mae mse rmse mape 0.41276110 0.24190200 0.49183534 0.08399855 }}} ==== t-sql로 구현한 단순회귀분석 ==== {{{ select (count(*)*sum(x*y)-(sum(x)*sum(y)))/(count(*)*sum(x*x)-(sum(x)*sum(x))) slope , avg(y)-((count(*)*sum(x*y))-(sum(x)*sum(y)))/((count(*)*sum(power(x,2)))-power(sum(x),2))*avg(x) intercept from tbl_xy }}} {{{ create view v_test as select 98 x, 1 y union all select 90, 2 union all select 90, 3 union all select 65, 4 union all select 66, 5 union all select 67, 6 union all select 13, 7 union all select 44, 8 union all select 27, 9 go declare @x nvarchar(255) , @y nvarchar(255) , @tablenm nvarchar(255) , @where nvarchar(1000) , @a nvarchar(255) , @b nvarchar(255) , @sql nvarchar(max); set @x = 'x'; set @y = 'y'; set @tablenm = 'select * from v_test'; set @where = ''; set @sql = ' declare @cnt int , @a decimal(38,5) , @b decimal(38,5) , @c decimal(38,5) , @d decimal(38,5) , @xavg decimal(38,5) , @yavg decimal(38,5) , @a1 decimal(38,5) , @b1 decimal(38,5); select @cnt = count(*) , @a = sum(convert(decimal(38,5), ' + @x + ')) , @b = sum(convert(decimal(38,5), ' + @y + ')) , @c = sum(convert(decimal(38,5), power(' + @x + ', 2))) , @d = sum(convert(decimal(38,5), ' + @x + '*' + @y + ')) , @xavg = avg(convert(decimal(38,5), ' + @x + ')) , @yavg = avg(convert(decimal(38,5), ' + @y + ')) from (' + @tablenm + ') as t; set @a1 = ((@b*@c) - (@a*@d)) / (@cnt*@c - power(@a, 2)); set @b1 = ((@cnt*@d) - (@a*@b)) / (@cnt*@c - power(@a, 2)); select 회귀식 , round(r2, 4) [R-squared] , round(1-(1-r2)*(@cnt-1)/(@cnt-2),4) [Adjusted R-squared] , round(corr, 4) 상관계수 from ( select ''y = ''+ convert(nvarchar, @a1) + '' + '' + convert(nvarchar, @b1) + ''x'' 회귀식 , 1 - var(' + @y + '-(@a1 + @b1*' + @x + ')) / var(' + @y + ') r2 , ((1.00/(@cnt - 1)*sum((' + @x + '-@xavg) * (' + @y + '-@yavg))) / sqrt(var(' + @x + ') * var(' + @y + '))) corr from (' + @tablenm + ') as t where 1=1 ' + @where + ' ) t select row_number() over(order by ' + @x + ') seq , ' + @x + ' x , ' + @y + ' y , power(' + @x + ', 2) x2 , ' + @x + '*' + @y + ' xy , @a1 + @b1*' + @x + ' estimatedvalue_linear --, @a1 * log(@a) + @b1 *' + @x + ' estimatedvalue_log , ' + @y + '-(@a1 + @b1*' + @x + ') diffvalue from (' + @tablenm + ') as t where 1=1 ' + @where + ' order by 2; '; print @sql; exec(@sql); }}} attachment:회귀분석/lm01.gif?width=20% attachment:회귀분석/lm02.gif?width=25% --그림 출처:http://www.sqlservercentral.com/articles/SQLCLR/71942/ ==== Two-Stage Least Squares ==== y ~ x 모델이 있다. (이 모델은 오차가 좀 있다) 그런데 x는 k에 좌지우지 된다. y랑은 상관없고. 그러면.. * x' ~ k 와 같이 x를 추정(x')하고, * y ~ x' 와 같이 추정된 x인 x'를 모델에 넣는다. 하면 오차가 줄어들게 된다. 이렇게 2단계로 regression하는 것을 Two-Stage Least Squares(2SLS)라고 한다. 여기서, k는 도구 변수(instrument variable)라고 한다. ==== 참고자료 ==== * [https://cran.r-project.org/web/packages/sjPlot/vignettes/plot_interactions.html Plotting Interaction Effects of Regression Models] * [http://en.wikibooks.org/wiki/R_Programming/Linear_Models R Programming/Linear Models] * [Linear Regression]