응용통계학 (6주차) 4월 14일
모형 진단
- 응용통계학 : 모형 진단
- Checking error assumptions
- Finding unusual observations
- Multicollinearity (다중공선성)
- 다중공선성의 해결
- Robust regression
- M-estimation
선형모형에서 추정량의 특정 성질이나 이후 추론과정들은 모형의 가정에 의존하게 된다. 즉, 가정에 대한 검토과정이 정당성을 획득하기 위해 필요한 절차이다. 모형 가정에 대한 이슈는 다음과 같이 요약될 수 있다.
- $y = X\beta + \epsilon$
- 오차항 가정
- 모형가정
- 오차항 : 정규성, 등분산성, 독립성
- stochastic part; Stochastic refers to the property of being well described by a random probability distribution.
- 모형 : $ E(y)=Xβ $ : 선형성
- deterministric part = symmetric part
- 특이 관측치 : 간혹 소수의 관측치들이 모형에 잘 맞지 않아 모형 전체에 영향을 미칠 수 있다.
- outlier 있나?
- 어떻게 처리해서 모델을 봐 볼까?
오차항을 관측 불가하므로 잔차에 기반하여 가정을 검토한다.
$$e=y−X\hat{β}=(I−H)y.$$
오차항은 서로 독립이더라도 잔차들은 그렇지 않음을 유의하자.
$e_i$ vs $\hat{y}_i$(x축 vs y축) 산점도(residual plot)가 가장 기본이 되는 진단 방법이다. 만약 가정이 모두 만족된다면 수직축 기준으로 대칭적이고 일정한 정도를 보이는 변동만이 나타날 것이다.
- Note : 가로축에 설명변수를 넣은 산점도를 이용하기도 한다.
- Simple linear에서 $\hat{y} = \hat{\beta}_0 + \hat{\beta}x$ 이니까 $\epsilon$ vs $\hat{y}$ 나 $\epsilon$ vs $x$ 나 거의 같다.
library(faraway)
data(savings,package="faraway")
str(savings)
-
sr: savings rate - personal saving divided by disposable income -
pop15: percent population under age of 15 -
pop75: percent population over age of 75 -
dpi: per-capita disposable income in dollars -
ddpi: percent growth rate of dpi
lmod <- lm(sr ~ pop15+pop75+dpi+ddpi,savings)
plot(fitted(lmod),residuals(lmod),xlab="Fitted",ylab="Residuals")
abline(h=0)
Q. What do you see in this plot?
- fitted value 가장 끝(오른)쪽 데이터 갈수록 작아지는 거 아닌가? 의심, 하지만 뒤로 갈수록 데이터가 없으니 모른다..$\rightarrow $ 주관적
library(car)
ncvTest(lmod)
등분산 검정 패키지,
- $H_0$ : 등분산, $H_1$ : 이분산
- p값을 보니 유의수준 5%에서 기각 못하니, 귀무가설 기각한다는 증거가 없다.
- 따라서 등분산이다.
문제가 발견되었을 때 몇가지 해결책
- 이분산성과 함께 비선형 추세가 보이는 경우 설명변수에 대한 적절한 변환이 필요할 수 있다.
- log 등
- 하지만 해석이 달라지지 않을까?
- 이분산성만이 순수하게 관측되는 경우 가중최소제곱법 활용가능.
- 최소제곱법에 비해 추정량 성질이 좋아지고, 모형 가정도 만족한다.
- 하지만 실제 활용이 어렵다.
- $\sum \frac{1}{w}(y_i - X\beta)^2$ 이런 식으로..하지만 등분산 가정이 깨진다.
- $\frac{1}{w} \rightarrow var(y_i) = \sigma^2 w_i$
- 반응변수에 대한 변환을 통해 등분산을 만족하도록 할 수 있다.
- $h(y)$ 분산 안정화 변환 $\rightarrow$ 쓰는 분야만 쓰는 경향.. 너무 전통적인 방법이다. 변환하면 해석이 달라져 모델을 세운 의미가 없을 수 있다.
gala: Species diversity on the Galapagos Islands
-
Species: the number of plant species found on the island -
Endemics: the number of endemic species -
Area: the area of the island (km$^2$) -
Elevation: the highest elevation of the island (m) -
Nearest: the distance from the nearest island (km) -
Scruz: the distance from Santa Cruz island (km) -
Adjacent: the area of the adjacent island (square km)
data(gala, package="faraway")
par(mfrow=c(1,2))
lmod <- lm(Species ~ Area + Elevation + Scruz + Nearest + Adjacent, gala)
plot(fitted(lmod),residuals(lmod),xlab="Fitted",ylab="Residuals")
abline(h=0)
lmod <- lm(sqrt(Species) ~ Area + Elevation + Scruz + Nearest + Adjacent, gala)
plot(fitted(lmod),residuals(lmod),xlab="Fitted",ylab="Residuals")
abline(h=0)
루트변환에 의한 분산 안정화
Q. 분산안정화가 꼭 필요한가? 가중최소제곱법을 고려함으로써 얻게 되는 이득은 무엇인가?
- (왼쪽 그림) 데이터 그대로 보니 비선형 trend가 있는 것 같다.
- 루트변환을 통해 분산 안정화를 시도했다.
- (오른쪽 그림) trend가 없어졌으니 등분산 만족하나?
- 모형이 달라지니 해석도 달라지지 않을까?
검정절차와 구간추정을 위해서는 정규성 가정에 대한 검토가 어느 정도 필요하다 (대표본의 경우 크게 중요하지는 않지만). Q-Q plot에 의해 정규성을 간단히 살펴보거나 잔차에 대한 정규성 검정을 통해 판단한다. Q-Q plot은 관측된 표본분위수가 정규성 가정 하에서의 이상적인 분위수와 일치하는지를 그림으로 나타낸 것이다. 즉, 직선 상에 분포해야 정규성에서 벗어나지 않는 것으로 볼 수 있다.
qqnorm, qqline 에서 qq의 의미 : quantile + theotical quantile + norm or line
qqplot
- 표준 100개 뽑았을때 순서대로 나열하면 가운데 값(50번째)는 0에 가까워야 해,
- 아니라면?
- 정규분포로 보기 어렵다고 판단.
- x축은 이론적 분위수, y축은 관측 분위수
par(mfrow=c(1,2))
lmod <- lm(sr ~ pop15+pop75+dpi+ddpi,savings)
qqnorm(residuals(lmod),ylab="Residuals",main="")
qqline(residuals(lmod))
hist(residuals(lmod),xlab="Residuals",main="")
Q. Do you think it is truely normal?
- (왼쪽 그림) 정규분포의 이상 직선에 가까운가?
- (오른쪽 그림) 히스토그램 그려보니 정규분포 모양에 비해 약간 찌그려져 있는 것 같다
- 모두 주관적...
몇 가지 분포로부터 표본을 생성하여 Q-Q plot의 모양을 살펴보자.
par(mfrow=c(1,3))
n <- 50
for(i in 1:3) {x <- rnorm(n); qqnorm(x); qqline(x)}
- light tail
- 실제 정규분포를 따라도, qqplot이 안 맞을 수도 있다.
- 심하게 벗어나지 않으면 정규분포 아니라고 안 해도 뭐..
par(mfrow=c(1,3))
for(i in 1:3) {x <- exp(rnorm(n)); qqnorm(x); qqline(x)}
- $log X \sim N$
- $X \sim exp(N)$
par(mfrow=c(1,3))
for(i in 1:3) {x <- rcauchy(n); qqnorm(x); qqline(x)}
- $-\infty \sim \infty$
- rcauchy - 전형적인 정규분포는 벗어나..
- 정규분포와 흡사하지만 굉장히 다르다.
- heavy tail을 가진 cauchy 분포
- 오차가 +보다 크게, -보다 작게 극단적으로 나타나네?
정규분포와 cauchy 분포 비교
n=100
x = rnorm(n)
y = rcauchy(n)
par(mfrow=c(1,2))
hist(x)
hist(y)
sort(y)[1:10]
sort(y)[90:100]
par(mfrow=c(1,3))
for(i in 1:3) {x <- runif(n); qqnorm(x); qqline(x)}
- 상하한이 정해져 있다.
- light tail 오차 발생에 함께 존재
- 0에서 1까지
Q-Q plot에서 나타나는 문제는 보통 heavy-tail 분포에서 기인하는 경우가 많으나 때로는 몇 개의 이상치로 인한 것일 수도 있다. 후자의 경우에는 해당 관측치의 제외여부를 검토할 필요가 있다.
- heavy tail 분포
-
- $\epsilon$ 이 heavy tail인가?
-
- outlier가 존재하나?
-
- outlier로 볼때, 해석에 따라 다르게 처리할 수 있다.
- 이상치를 포함할건지, 말건지..
정규성 가정이 위배되더라도 최소제곱추정량은 여전히 BLUE(best linear unbiased estimator)이다. 하지만 최적의 추정량은 아닐 수 있으며 robust 추정량이 더 나은 성질을 나타내기도 한다. 또한 검정이나 구간추정의 부정확성을 초래할 수 있다. 하지만 반복해서 언급하였듯이 대표본의 경우 크게 문제가 되지 않는 경우가 많다.
정규성에 위배가 나타날 경우 몇 가지 해결책은 다음과 같다.
- Short-tail(=light-tail)의 경우 파급이 크지 않아 크게 신경쓰지 않아도 무방하다.
- 비대칭 분포의 경우 적절한 변환이 필요할 수 있다.
- 해석에 영향이 가니까 해석이 중요하면 조심해야 한다.
- Heavy-tail의 경우 bootstrap이나 permutation test와 같은 방식으로 추론하거나 robust method를 이용한다.
- bootstrap이나 permutation test는 최소제곱 분포를 사용하지만 robust method를 이용하면 추정방식 자체를 바꾸겠다는 것.
다음과 같은 정규성 검정도 가능하다.
shapiro.test(residuals(lmod))
$H_0$ : 정규성 가정 만족
$H_1$ : 정규성 가정 만족하지 않음
p값이 유의수준 5%보다 크기 때문에 정규성 가정을 만족한다고 할 수 있다.
독립성이 위배되는 것은 그 경우의 수가 매우 많아 확인이 매우 어렵다. 만약, 시간 순으로 관측된 자료라면 시간 순으로 일종의 상관성이 존재할 가능성이 있다.
다음은 북반구의 온도와 그에 대한 8가지 proxy에 의한 선형모형 적합결과이다.
A data frame with 1001 observations on the following 10 variables.
-
nhtemp: Northern hemisphere average temperature (C) provided by the UK Met Office (known as HadCRUT2) -
wusa: Tree ring proxy information from the Western USA. -
jasper: Tree ring proxy information from Canada. -
westgreen: Ice core proxy information from west Greenland -
chesapeake: Sea shell proxy information from Chesapeake Bay -
tornetrask: Tree ring proxy information from Sweden -
urals: Tree ring proxy information from the Urals -
mongolia: Tree ring proxy information from Mongolia -
tasman: Tree ring proxy information from Tasmania -
year: Year 1000-2000AD
각각 동일연도의 관측값들임, 따라서 시간순으로 상관성이 있을 가능성이 존재한다.
data(globwarm,package="faraway")
lmod <- lm(nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals + mongolia + tasman, globwarm)
plot(residuals(lmod) ~ year, na.omit(globwarm), ylab="Residuals")
abline(h=0)
잔차가 같은 부호가 연속되는 경향이 있음을 알 수 있다. 이런 현상은 오차가 양의 상관성을 띌 때 전형적으로 나타난다.
$$Corr(\epsilon_i,\epsilon_{i+1})>0$$
n <- length(residuals(lmod))
plot(tail(residuals(lmod),n-1) ~ head(residuals(lmod),n-1), xlab= expression(hat(epsilon)[i]),ylab=expression(hat(epsilon)[i+1]))
abline(h=0,v=0,col=grey(0.75))
그림만 보면 전형적인 양의 상관관계에 있음, 독립성 깨진 것.
Durbin-Watson test 잔차가 AR(1) 프로세스인지를 검정한다.
$$\epsilon_{t+1} = \rho \epsilon _t + v_t$$ $$H_0:ρ=0 \text{ vs } H_1:ρ≠0$$
- $\rho>0$ 이면 양의 자기 상관이 있다.
- $\rho<0$ 이면 음의 자기 상관이 있다.
- $0 < \rho < 1$
- 잔차 $v_t \sim^{iid} N(0,\sigma^2)$
- $H_0$ : 독립성에 문제 없다 정도, 확신은 아냐..
- but 독립적으로 관찰된 것들이면 독립성에 문제가 없다.
- $H_1$ : 독립성 깨진다.
require(lmtest)
dwtest(nhtemp ~ wusa + jasper + westgreen + chesapeake + tornetrask + urals + mongolia + tasman, data=globwarm)
더빈왓슨 검정의 p값이 유의수준 5%보다 작아서 귀무가설을 기각 시키고 독립성이 성립 안 한다는 결론.
- bootstrap 하거나
- AR(1) 모형에 포함시켜서 가정하거나
- outliers
- influential observations
- leverage points
위 개념들은 명확하게 구분되지 않는다.
$H$의 대각원소 $h_i$를 leverage라 한다. $var(\hat{\epsilon}_i) = \sigma^2( 1 − h_i)$ 로부터 $h_i$가 크면 해당관측치쪽으로 적합된 직선이 당겨지는 효과가 있게 된다. $h_i$는 $X$-space 에서 극단적인 쪽에 위치할 때 발생함이 알려져 있다.
- $0 \le h \le 1$
- $h$ 가 크면 $var(\epsilon)$이 작아지고, $\epsilon$이 0에 가까워진다.
$\sum _i h_i = K$ 이므로 $\frac{K}{N}$의 2배~4배를 넘는 관측치에 대해서 주의가 필요함이 알려져 있다.
lmod <- lm(sr ~ pop15 + pop75 + dpi + ddpi, savings)
hatv <- hatvalues(lmod)
head(hatv)
sum(hatv)
- 변수 5개
- 데이터 50개
- $\frac{5}{50} = 0.1 = h_i =$leverage
이상치는 데이터의 특성에서 다소 떨어진 관측치를 의미한다. 이상치는 전체 적합 결과에 영향을 미칠 수도 그렇지 않을 수도 있다.
cars_outliers1 <- data.frame(speed=c(20), dist=c(220))
cars_outliers2 <- data.frame(speed=c(30), dist=c(130))
cars_outliers3 <- data.frame(speed=c(35), dist=c(30))
# 기존 데이터에 이상치 데이터를 추가한 데이터 프레임
new_cars1 <- rbind(cars, cars_outliers1)
new_cars2 <- rbind(cars, cars_outliers2)
new_cars3 <- rbind(cars, cars_outliers3)
# plot 창 변환
par(mfrow = c(1, 3))
# 시각화
plot(new_cars1$speed, new_cars1$dist, xlim=c(0, 28), ylim=c(0, 230),xlab="x", ylab="y", pch="*", col="red", cex=2)
abline(lm(dist ~ speed, data=new_cars1), col="blue", lwd=3, lty=2)
abline(lm(dist ~ speed, data=cars), col="black", lwd=3, lty=2)
plot(new_cars2$speed, new_cars2$dist, xlim=c(0, 28), ylim=c(0, 230),xlab="x", ylab="y", pch="*", col="red", cex=2)
abline(lm(dist ~ speed, data=new_cars2), col="blue", lwd=3, lty=2)
abline(lm(dist ~ speed, data=cars), col="black", lwd=3, lty=2)
plot(new_cars3$speed, new_cars3$dist, xlim=c(0, 28), ylim=c(0, 230),xlab="x", ylab="y", pch="*", col="red", cex=2)
abline(lm(dist ~ speed, data=new_cars3), col="blue", lwd=3, lty=2)
abline(lm(dist ~ speed, data=cars), col="black", lwd=3, lty=2)
- 위로 올라온 그래프
- 조금 변한 그래프
- 밑으로 내려온 그래프
소수의 이상치는 큰 자료에서는 큰 문제가 되지 않을 수 있다.
이상치 발견의 경우
- 데이터 입력 오류를 살펴보아라.
- 배경을 살펴보아라. 때로는 이상치의 발견이 의미있을 수도 있다.
- 이상치가 데이터의 특성을 반영하지 않는다고 생각하여 제거하고자 하는 경우 어떤 것을 제거하는지 기록하고 보고한다.
- 이상치가 자연스러운 데이터 생성과정에서 발생한 것이고 그 수가 적지 않다면 모형 적합과정에서 포함하되 robust method 등을 이용하는 것이 좋다.
Ex :
stardataset
data(star, package="faraway")
str(star)
별 밝기/온도
plot(star$temp,star$light,xlab="log(Temperature)", ylab="log(Light Intensity)")
lmod <- lm(light ~ temp, star)
abline(lmod)
4개의 이상치!
range(rstudent(lmod)) # `rstudent` function computs the studentized residuals distributed as t(n-p)
lmodi <- lm(light ~ temp, data=star, subset=(temp>3.6))
plot(star$temp,star$light,xlab="log(Temperature)", ylab="log(Light Intensity)")
abline(lmod)
abline(lmodi,lty=2)
Big changes!
영향점은 제거여부가 모형의 적합에 큰 변화를 초래하는 것을 말한다. 이상치나 leverage가 영향점이 될 수 있다.
Cook statistics : changes in β with and without ith observation
$$D_i = \frac{1}{K} \epsilon^{2}_{i} \frac{h_i}{1−h_i}$$
$h_i$가 1에 가까우면 $\beta$에 대한 변화가 커진다.
- leverage이지만 영향점일수도, 아닐수도 있다.
$\epsilon$이 커질수록 $\beta$에 대한 변화가 커진다.
$h_i$가 1에 가깝게 크더라도 $\epsilon$이 굉장히 작으면 영향이 없을 것이다.
lmod <- lm(sr ~ pop15+pop75+dpi+ddpi,savings)
cook <- cooks.distance(lmod)
plot(cook)
lmodi <- lm(sr ~ pop15+pop75+dpi+ddpi,savings,subset=(cook < max(cook)))
sumary(lmodi) # 50번째쯤에 존재하는 영향치를 제거한 모델
sumary(lmod) # 원래 모델
R에서는 기본적인 오차항에 대한 가정 검토를 위한 그림을 제공한다.
plot(lmod)
1번 그림
- outlier 있지만 등분산에 크게 문제는 없어보인다.
2번 그림
- outlier 의심값들이 존재하지만, 직선에 거의 근접한 모습이다.
3번 그림
- pass..
4번 그림
- Libya 가 가장 큰 통계량을 가지고 있지만 아주 심한 영향은 아니다.
- cook 통계량은 1이 넘어가면 문제 발생!
잔차는 선형성을 제외한 다른 체계적인 요소가 남아있는지를 판단할 수 있는 근거가 된다.
만약 특정 예측변수(설명변수)가 다른 변수들의 선형결합으로 표현되면 $X^\top X$의 역행렬이 존재하지 않아 최소제곱추정량이 정의되지 않는다.
- 설명변수 간 상관성이 크게(강하게) 존재한다면?
완전히 표현되지 않더라도 $X^\top X$가 singular에 가까울 수 있다. 이 경우 추정량은 정의되지만 분산이 매우 커지는 결과를 초래한다. 보통 변수들간에 상관성이 존재할 때 이런 일이 발행하지만 변수의 개수가 많아지면 거의 필연적으로 이런 현상이 나타남이 알려져 있다.
- $Var(\hat{\beta}) = \sigma^2(X^{-1} X)^{-1}$
- 여기서 $(X^{-1} X)^{-1}$는 존재하지 않아도 문제고, 작아도 문제다 $\to$ 분산이 커지게 되니까!
분산의 증가는 부정확한 추정을 의미할 뿐 아니라 계수에 대한 검정의 검정력을 떨어뜨리게 된다.
이러한 다중공선성은 $X^\top X$의 고유치(eigenvalue)로부터 확인이 가능하다. 고유치가 0인 것이 있으면 exact collinearity를 나타내고 0에 가까운 것이 있으면 다중공선성이 어느 정도 존재하는 것으로 볼 수 있다. 다음과 같은 condition number로 기준을 삼기도 하는데 10 혹은 30이 넘으면 주의를 요하는 것으로 판단한다.
- 10 혹은 30의 범위의 기준은 절대적인 것이 아니다.
- $A$: 고유치 $>= 0$
- $A^{-1}$ 존재 $\to \lambda_i >0$
$$\kappa = \sqrt{ \frac{\lambda_{max}} {\lambda_{min}} } $$
$\lambda_{min}$이 0에 가까우면 $\kappa$가 커지겠지?
다음 식을 살펴보자.
$$var(\hat{\beta}_j) = \sigma^2 \big( \frac{1}{1-R^{2}_{j}} \big) \frac{1}{\sum_i (x_{ij} - \bar{x}_j)^2}$$
여기서 $R^{2}_{j}$는 $j$번째 설명변수를 반응변수로 두고 나머지 설명변수로 적합한 모형의 결정계수이다. 여기서 $\frac{1}{(1−R^{2}_{j})}$를 분산팽창계수라 한다 (VIF; Variance Inflation Factor). (rull of thumb : VIF > 10 → multicollinearity)
- $R^{2}_{j} = \frac{SSR}{SST}$
- 성명변수 $X_j~X_1, \dots X_{j+1}$
- 설명력이 크다면 $X_j$를 대체할 수도 있다. $X_j$ 없이 잘 설명되니까!
- 만약 분산팽창계수가 10을 넘으면 빼고 모형을 설정해보는 것이 좋은 시도!
만약 모든 설명변수들을 직교하도록 만들 수 있다면 공선성의 측면에서는 바람직할 것이다. 그래서 실험설계를 할 때 이러한 요소를 고려하기도 한다.
- 직교; $corr = 0 \to R^{2}_{j} = 0$
The dataset contains the following variables
-
Age:Age in years -
Weight: Weight in lbs -
HtShoes:Height in shoes in cm -
Ht:Height bare foot in cm -
Seated:Seated height in cm -
Arm:lower arm length in cm -
ThighThigh length in cm -
Leg: Lower leg length in cm -
hipcenter: horizontal distance of the midpoint of the hips from a fixed location in the car in mm
data(seatpos, package="faraway")
lmod <- lm(hipcenter ~ ., seatpos)
summary(lmod)
- 모든 설명변수가 유의수준 5%를 넘는다.
- 하지만 모델 자체의 p값은 5%보다 낮았다.
설명변수들간의 상관행렬을 살펴보자.
round(cor(seatpos[,-9]),2)
Seated 변수는 HtShoes와 Ht 와 상관관계가 높았다. 각각 0.93
$X^{\top}X$의 고유치 또한 살펴보자.
x <- model.matrix(lmod)[,-1]
e <- eigen(t(x) %*% x)
e$val
- X에 절편 있으면 여기서는 지워도 상관이 없을 것.
-
e$val은 max(N)임
sqrt(e$val[1]/e$val)
library(regclass)
VIF(lmod)
vif(lmod)
다중공선성의 존재가 의심된다.
해결하는 간단한 방법 중 하나는 일부 변수들을 제거하는 것이다.
lmod2 <- lm(hipcenter ~ Age + Weight + Ht, seatpos)
summary(lmod2)
- 수정된 설명력이 0.6258로, 위 모형보다는 효율적이다.
결정계수 측면에서 큰 차이가 없다.
변수를 제거하지 않고 개선하는 방법 중 Ridge나 Lasso와 같은 축소추정법을 이용할 수도 있다.
- 분산 크기를 줄여야 겠다.
오차항이 정규분포가 아니거나 이상치가 다수 존재하는 경우 최소제곱추정량에 대한 대안이 필요할 수 있다.
- 분산 안정화
- 최소제곱법: 모든 값의 제곱합
- 로버스트: 중요도 다르게 가중치 주기
다음과 같은 식을 최소로 하는 계수를 찾는 방법
$$sum^{n}_{i=1} \rho (y_i - x^{\top}_{i} \beta)$$
- $= \sum_{i=1} (y_i - X_i \beta)^2$
Some possible choices:
- $\rho(x)=x^2$ : least squares
- ex) $\bar{X} = argmin_{\theta} \sum (x_i - \theta)^2$
- 더 강건(robust)한 $L_2$
- $\rho(x)=|x|$ : least absolute deviation regression or $L_1$ regression
- ex) $m = argmin_{\theta} \sum^{n}_{i=1} |x_i - \theta|^2$
- 차이 커지면 오차가 기하급수적으로 증가
- 극단값 많이 있으면 설명력이 낮아짐
- $\rho (x)=\frac{x^2}{2}⋅I(|x|≤c)+(c|X|−\frac{c^2}{2})⋅I(|X|>c)$ : Huber's method and this is a compromise between $L_1$ and least squares.
- 절대값에 증거한 측정이 강건(robust)한 계산
data(gala, package="faraway")
lsmod <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,gala)
summary(lsmod)
다음은 Huber의 방법이다.
- 밑에서 rlm의 r은 robust 이다.
library(MASS)
rlmod <- rlm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,gala) # Huber's method
summary(rlmod)
다음은 적합에서 사용된 가중치를 작은 순으로 출력하여 본 것이다.
wts <- rlmod$w
names(wts) <- row.names(gala)
head(sort(wts),10)
LAD 방법도 다음과 같이 적합 가능하다.
library(quantreg)
l1mod <- rq(Species ~Area+Elevation+Nearest+Scruz+Adjacent, data=gala)
summary(l1mod)
계수추정치에 변화가 다소 있으며 adjacent 가 유의미하지 않음을 시사한다.