R Language
선형 모델 (회귀)

통사론

• lm (수식, 데이터, 부분 집합, 가중치, na.action, 메서드 = "qr", 모델 = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, .. .)

매개 변수

mtcars 데이터 집합의 선형 회귀

내장 된 mtcars 데이터 프레임 에는 무게, 연료 효율 (갤런 당 마일), 속도 등 32 개 차량에 대한 정보가 들어 있습니다 (데이터 세트에 대해 자세히 알아 보려면 help(mtcars) ).

연료 효율 ( mpg )과 중량 ( wt ) 간의 관계에 관심이 있다면 다음 변수를 사용하여 이러한 변수를 계획 할 수 있습니다.

plot(mpg ~ wt, data = mtcars, col=2)

플롯은 (선형) 관계를 보여줍니다!. 선형 모델의 계수를 결정하기 위해 선형 회귀를 수행하려면, lm 함수를 사용합니다.

fit <- lm(mpg ~ wt, data = mtcars)

여기서 ~ 는 "설명 됨"을 의미하므로 공식 mpg ~ wt 는 mpg ~ wt 로 설명 된 mpg를 예측한다는 것을 의미합니다. 출력을 보는 가장 유용한 방법은 다음과 같습니다.

summary(fit)

출력 결과는 다음과 같습니다.

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

다음 정보를 제공합니다.

• mpg의 최적 예측을 제안하는 각 계수 ( wt 및 y 절편)의 추정 기울기는 37.2851 + (-5.3445) * wt
• 절편과 가중치가 우연히 발생하지 않았 음을 나타내는 각 계수의 p- 값
• R ^ 2와 조정 된 R ^ 2와 같은 적합도에 대한 전반적인 추정치는 모델에 의해 mpg 의 변화가 얼마나 많은지를 보여줍니다

예측 된 mpg 를 표시하는 첫 번째 플롯에 선을 추가 할 수 있습니다.

abline(fit,col=3,lwd=2)

그 플롯에 방정식을 추가하는 것도 가능합니다. 먼저 coef 로 coef . 그런 다음 paste0 을 사용하여 적절한 변수와 +/- 사용하여 계수를 축소하여 방정식을 작성합니다. 마지막으로 mtext 사용하여 플롯에 추가합니다.

bs <- round(coef(fit), 3) 
lmlab <- paste0("mpg = ", bs[1],
             ifelse(sign(bs[2])==1, " + ", " - "), abs(bs[2]), " wt ")
mtext(lmlab, 3, line=-2) 

결과는 다음과 같습니다.

회귀 분석 (기본)

계속해서 mtcars 예제에서, 여기에 게시에 적합 할 수있는 선형 회귀의 플롯을 생성하는 간단한 방법이 있습니다.

먼저 선형 모델을 맞추고

fit <- lm(mpg ~ wt, data = mtcars)

그런 다음 관심있는 두 변수를 플롯하고 정의 도메인 내에 회귀선을 추가합니다.

plot(mtcars$wt,mtcars$mpg,pch=18, xlab = 'wt',ylab = 'mpg')
lines(c(min(mtcars$wt),max(mtcars$wt)),
as.numeric(predict(fit, data.frame(wt=c(min(mtcars$wt),max(mtcars$wt))))))

거의 다 왔어! 마지막 단계는 줄거리, 회귀 방정식, rsquare 및 상관 계수를 더하는 것입니다. 이것은 vector 함수를 사용하여 수행됩니다.

rp = vector('expression',3)
rp[1] = substitute(expression(italic(y) == MYOTHERVALUE3 + MYOTHERVALUE4 %*% x), 
          list(MYOTHERVALUE3 = format(fit$coefficients[1], digits = 2),
                        MYOTHERVALUE4 = format(fit$coefficients[2], digits = 2)))[2]
rp[2] = substitute(expression(italic(R)^2 == MYVALUE), 
             list(MYVALUE = format(summary(fit)$adj.r.squared,dig=3)))[2]
rp[3] = substitute(expression(Pearson-R == MYOTHERVALUE2), 
             list(MYOTHERVALUE2 = format(cor(mtcars$wt,mtcars$mpg), digits = 2)))[2]

legend("topright", legend = rp, bty = 'n')

벡터 함수를 적용하여 RMSE와 같은 다른 매개 변수를 추가 할 수 있습니다. 10 개의 요소로 된 전설이 필요하다고 상상해보십시오. 벡터 정의는 다음과 같습니다.

rp = vector('expression',10)

r[1] .... ~ r[10] 을 정의해야합니다 r[10]

출력은 다음과 같습니다.

가중치

때때로 우리는 모델이 다른 데이터보다 데이터 포인트 나 예제에 더 많은 가중치를 주길 원합니다. 이는 모델을 학습하는 동안 입력 데이터의 가중치를 지정하여 가능합니다. 예제에 대해 비 균일 가중치를 사용할 수있는 시나리오에는 일반적으로 두 가지 종류가 있습니다.

• 분석 가중치 (Analytic Weights) : 서로 다른 관측치의 정밀도를 반영합니다. 예를 들어, 각 관찰이 지리적 영역의 평균 결과 인 데이터를 분석하는 경우 분석 가중치는 추정 된 분산의 역수에 비례합니다. 관측 수를 기준으로 비례 가중치를 제공하여 데이터의 평균을 처리 할 때 유용합니다. 출처

• 표본 가중치 (Inverse Probability Weights - IPW) : 데이터가 수집 된 것과 다른 모집단으로 표준화 된 통계를 계산하기위한 통계 기법. 서로 다른 샘플링 모집단과 목표 추론 (대상 모집단)을 가진 연구 설계는 적용 분야에서 일반적입니다. 누락 된 값이있는 데이터를 다룰 때 유용합니다. 출처

lm() 함수는 분석 가중치를 수행합니다. 샘플링 가중치의 경우 survey 패키지를 사용하여 설문 조사 설계 객체를 작성하고 svyglm() 실행합니다. 기본적으로 survey 패키지는 샘플링 가중치를 사용합니다. (참고 : lm() 및 svyglm() 과 family gaussian() 은 모두 가중치가 최소 제곱을 최소화하여 계수를 계산하기 때문에 동일한 표준 지점 추정치를 산출합니다. 표준 오차를 계산하는 방법은 다릅니다.

테스트 데이터

data <- structure(list(lexptot = c(9.1595012302023, 9.86330744180814, 
8.92372556833205, 8.58202430280175, 10.1133857229336), progvillm = c(1L, 
1L, 1L, 1L, 0L), sexhead = c(1L, 1L, 0L, 1L, 1L), agehead = c(79L, 
43L, 52L, 48L, 35L), weight = c(1.04273509979248, 1.01139605045319, 
1.01139605045319, 1.01139605045319, 0.76305216550827)), .Names = c("lexptot", 
"progvillm", "sexhead", "agehead", "weight"), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -5L))

분석 가중치

lm.analytic <- lm(lexptot ~ progvillm + sexhead + agehead, 
                            data = data, weight = weight)
summary(lm.analytic)

산출

Call:
lm(formula = lexptot ~ progvillm + sexhead + agehead, data = data, 
    weights = weight)

Weighted Residuals:
         1          2          3          4          5 
 9.249e-02  5.823e-01  0.000e+00 -6.762e-01 -1.527e-16

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.016054   1.744293   5.742    0.110
progvillm   -0.781204   1.344974  -0.581    0.665
sexhead      0.306742   1.040625   0.295    0.818
agehead     -0.005983   0.032024  -0.187    0.882

Residual standard error: 0.8971 on 1 degrees of freedom
Multiple R-squared:  0.467, Adjusted R-squared:  -1.132 
F-statistic: 0.2921 on 3 and 1 DF,  p-value: 0.8386

샘플링 가중치 (IPW)

library(survey)
data$X <- 1:nrow(data)             # Create unique id
        
# Build survey design object with unique id, ipw, and data.frame
des1 <- svydesign(id = ~X,  weights = ~weight, data = data)
        
# Run glm with survey design object
prog.lm <- svyglm(lexptot ~ progvillm + sexhead + agehead, design=des1)

산출

Call:
svyglm(formula = lexptot ~ progvillm + sexhead + agehead, design = des1)

Survey design:
svydesign(id = ~X, weights = ~weight, data = data)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 10.016054   0.183942  54.452   0.0117 *
progvillm   -0.781204   0.640372  -1.220   0.4371  
sexhead      0.306742   0.397089   0.772   0.5813  
agehead     -0.005983   0.014747  -0.406   0.7546  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2078647)

Number of Fisher Scoring iterations: 2

다항식 회귀 분석을 통한 비선형 성 검사

때로는 선형 회귀로 작업 할 때 데이터의 비선형 성을 검사해야합니다. 이를 수행하는 한 가지 방법은 다항식 모델을 적용하고 선형 모델보다 데이터에 더 잘 맞는지 확인하는 것입니다. 변수 관계가 본질적으로 다항식이라고 믿어지기 때문에 이차 또는 고차 모델에 적합하다는 것을 나타내는 이론적 인 이유와 같은 다른 이유가 있습니다.

mtcars 데이터 세트에 2 차 모델을 적용 해 mtcars . 선형 모델의 경우 mtcars 데이터 집합에서 선형 회귀를 참조하십시오.

먼저 변수 mpg (Miles / gallon), disp (Displacement (cu.in.)) 및 wt (Weight (1000 lbs))의 산점도를 만듭니다. mpg 와 disp 사이의 관계는 비선형으로 나타난다.

plot(mtcars[,c("mpg","disp","wt")])

선형 적합은 disp 가 중요하지 않음을 나타냅니다.

fit0 = lm(mpg ~ wt+disp, mtcars)
summary(fit0)

# Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 34.96055    2.16454  16.151 4.91e-16 ***
#wt          -3.35082    1.16413  -2.878  0.00743 ** 
#disp        -0.01773    0.00919  -1.929  0.06362 .  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.917 on 29 degrees of freedom
#Multiple R-squared:  0.7809,    Adjusted R-squared:  0.7658

그런 다음 2 차 모델의 결과를 얻기 위해 I(disp^2) . R^2 때 새로운 모델이 더 잘 나타나고 모든 변수가 중요합니다.

fit1 = lm(mpg ~ wt+disp+I(disp^2), mtcars)
summary(fit1)

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 41.4019837  2.4266906  17.061  2.5e-16 ***
#wt          -3.4179165  0.9545642  -3.581 0.001278 ** 
#disp        -0.0823950  0.0182460  -4.516 0.000104 ***
#I(disp^2)    0.0001277  0.0000328   3.892 0.000561 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.391 on 28 degrees of freedom
#Multiple R-squared:  0.8578,    Adjusted R-squared:  0.8426 

우리가 3 개의 변수를 가지기 때문에, 적합 모델은 다음과 같이 표현되는 표면이다.

mpg = 41.4020-3.4179*wt-0.0824*disp+0.0001277*disp^2

다항식 회귀를 지정하는 또 다른 방법은 매개 변수 raw=TRUE 인 poly 를 사용하는 것입니다. 그렇지 않으면 직교 다항식 이 고려됩니다 (자세한 내용은 help(ploy) 참조). 우리는 다음과 같은 결과를 얻는다.

summary(lm(mpg ~ wt+poly(disp, 2, raw=TRUE),mtcars))

마지막으로, 추정 된 표면의 플롯을 보여줄 필요가 있다면? 음, R 3D 플롯을 만드는 데는 여러 가지 옵션이 있습니다. 여기에 우리가 사용 Fit3d 에서 p3d 패키지로 제공된다.

library(p3d)
Init3d(family="serif", cex = 1)
Plot3d(mpg ~ disp+wt, mtcars)
Axes3d()
Fit3d(fit1)

품질 평가

회귀 모델을 작성한 후 결과를 확인하고 모델이 적절하고 현재의 데이터와 잘 작동하는지 여부를 결정하는 것이 중요합니다. 잔차 플롯 및 기타 진단 플롯을 검사하여이를 수행 할 수 있습니다.

# fit the model
fit <- lm(mpg ~ wt, data = mtcars)
# 
par(mfrow=c(2,1))
# plot model object
plot(fit, which =1:2)

이 플롯은 모델을 빌드하는 동안 만들어진 두 가지 가정을 검사합니다.

1. 예측 변수 (이 경우 mpg )의 예상 값은 예측 변수의 선형 조합 (이 경우 wt )에 의해 주어집니다. 우리는이 추정치가 편향 적이기를 기대합니다. 따라서 나머지는 예측 변수의 모든 값에 대해 평균을 중심으로 배치해야합니다. 이 경우 잔차는 양 끝에서 양수이고 중간에서 음수가되는 경향이있어 변수간에 비선형 관계가 있음을 알 수 있습니다.
2. 실제 예측 변수는 일반적으로 추정치를 중심으로 분산됩니다. 따라서 나머지는 정상적으로 분산되어야합니다. 정상적으로 분산 된 데이터의 경우, 정상 QQ 플롯의 점은 대각선에 있거나 대각선에 가깝습니다. 끝 부분에는 약간의 비뚤어 짐이 있습니다.

'예측'기능 사용

일단 모델이 만들어지면 새로운 데이터로 테스트하기위한 주요 기능이 predict 됩니다. 이 예제에서는 mtcars 내장 데이터 집합을 사용하여 변위에 대해 갤런 당 마일을 회귀합니다.

my_mdl <- lm(mpg ~ disp, data=mtcars)
my_mdl

Call:
lm(formula = mpg ~ disp, data = mtcars)

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

내가 변위가있는 새로운 데이터 소스를 가지고 있다면 갤런 당 추정 마일을 볼 수 있습니다.

set.seed(1234)
newdata <- sample(mtcars$disp, 5)
newdata
[1] 258.0  71.1  75.7 145.0 400.0

newdf <- data.frame(disp=newdata)
predict(my_mdl, newdf)
       1        2        3        4        5 
18.96635 26.66946 26.47987 23.62366 13.11381

이 프로세스의 가장 중요한 부분은 원본 데이터와 동일한 열 이름을 사용하여 새 데이터 프레임을 만드는 것입니다. 이 경우 원래 데이터에는 disp 라는 레이블이 지정된 열이 있었으므로 같은 이름의 새 데이터를 호출해야했습니다.

주의

몇 가지 일반적인 함정을 살펴 보겠습니다.

1. 새 객체에 data.frame을 사용하지 않습니다.

   predict(my_mdl, newdata)
   Error in eval(predvars, data, env) : 
      numeric 'envir' arg not of length one
   

2. 새 데이터 프레임에서 같은 이름을 사용하지 않음 :

   newdf2 <- data.frame(newdata)
   predict(my_mdl, newdf2)
   Error in eval(expr, envir, enclos) : object 'disp' not found
   

정확성

예측의 정확성을 확인하려면 새 데이터의 실제 y 값이 필요합니다. 이 예에서 newdf 는 'mpg'및 'disp'에 대한 열이 필요합니다.

newdf <- data.frame(mpg=mtcars$mpg[1:10], disp=mtcars$disp[1:10])
#     mpg  disp
# 1  21.0 160.0
# 2  21.0 160.0
# 3  22.8 108.0
# 4  21.4 258.0
# 5  18.7 360.0
# 6  18.1 225.0
# 7  14.3 360.0
# 8  24.4 146.7
# 9  22.8 140.8
# 10 19.2 167.6

p <- predict(my_mdl, newdf)

#root mean square error
sqrt(mean((p - newdf$mpg)^2, na.rm=TRUE))
[1] 2.325148