R - 선형회귀와 분산분석
R - 선형회귀와 분산분석
선형회귀
1) 모형이 통계적으로 유의미한가?
- F 통계량/p값 확인
2) 계수들이 유의미한가?
- Summary에서 해당 계수의 t 통계량과 p-값, 신뢰구간 확인
3) 모형이 유용한가?
- R^2 혹은 수정계수
4) 모형이 데이터를 잘 적합하고 있는가?
- resid( )
5) 데이터가 선형회귀가 전제하고 있는 가정을 만족시키는가?
2. 필요 커맨드
> anova( )
> coef( )
> deviance( ): 잔차제곱합
> effects( ): 직교효과들로 이루어진 벡터
> fitted( ): 적합된 y값으로 이루어진 벡터
> resid( ): 모형의 잔차
> vcov( ): 주 매개변수들의 분산-공분산 행렬
3. 절편이 없는 선형회귀 실시
먼저 신뢰구간에 0이 포함되는지 확인
> confint(lm (y ~ x) )
확인 후
> lm( y ~ x + 0)
4. 상호작용 항 / 그냥 곱셈 항
> lm ( y ~ u*v)
> lm ( y ~ u:v)
* 다음 3가지 식은 동일한 것이다.
y ~ u*v = y ~ u+v+u:v = y ~ (u+v)^2
* 표현식 쓰기
> lm ( y~ u + I(u^2) )
* 다항식 회귀분석
> lm ( y ~ poly(x,3, raw = TRUE) )
5. 모형 설계 - 둘다 이용할 것
1) 후진제거법 Backward Elimination
> reduced.model <- step(full.model, direction = "backward")
2) 전진선택법 Forward Selection
> min.model <- lm(y~1)
> fwd.model <- step(min.model, direction = "forward", scope = ( ~ x1 + x2 ... ),
trace = 0 )
6. 데이터의 부분집합에 대해 회귀분석
> lm (y ~ x, subset = 1:100)
subset = 1:floor (length(x)/2)
subset = ( lab == "NJ" )
7. 선형회귀 진단
1) 회귀 잔차 그래프
> m <- lm ( y ~ x)
> plot(m, which = 1)
- which=1이 있으면 잔차그래프만 선택
안하면 4개의 진단 그래프 생성
2) 4개의 진단 그래프
> plot(m)
> outlierTest(m) # car패키지에 있음
3) 영향력이 큰 관측 판별
> influence.measures(m)
4) 잔차의 자기상관 검정 (더빈-왓슨)
> dwtest(m) # lmtest패키지에 있음
> acf(m) # 모델 잔차의 ACF(자기상관함수)를 그래프로 그림
* 모형이 음의 자기상관을 드러낼 수 있는 종류라면
> dwtest(m, alternative = "two.sided")
8. 예측
1) 새로운 값 예측
> m <- lm (y ~ u+v+w)
> preds <- data.frame(u=3.1, v=4.0, w=5.5)
> predict (m, newdata = preds)
2) 예측 구간 구하기
> predict(m, newdata=preds, interval = "prediction")
- 이러한 예측구간들은 정규성으로부터의 편차에 극도로 민감하다.
따라서 만약 반응변수가 정규분포가 아니라는 판단이 되면,
부트스트랩과 같은 비모수적 방법으로 예측구간을 구해야 한다.
분산분석
[1] 모수 검정
1. 2 Sample t-test / paired 2 Sample t-test
1) 탐색적 분석
hist / boxplot / summary
2) 두 집단이 독립적인지 짝을 이루는지에 따라
t.test(x ~ f) 혹은 t.test(x, y)
인자: conf.level = / alternative = " " / var.equal = / data =
t.test(x1, x2, paired = TRUE)
* prop.test(x, n)
인자 conf.level = / alternative =
2. 일원분산분석 One-way ANOVA
1) oneway.test(x ~ f, data = ): 등분산을 가정하지 않음
- 하려면 var.equal = TRUE
2) aov(x ~ f, data = ): 등분산을 가정함
- bartlett.test(x ~ f, data = ) 로 등분산성 검정
(p값이 커야 오차항은 등분산성을 만족한다)
3) 사후검정 다중비교 (Post-hoc pair-wise multiple regression)
: 설정된 유의수준을 그대로 유지하면서,
가능한 모든 수준평균 간의 차이에 대한 사후검정 실시 (Tukey와 Duncan에 해당)
(1) TukeyHSD: Honestly Significan Difference (Studentized Range Distribution 이용)
> m <- aov(x ~ f)
> TukeyHSD(m)
> plot(TukeyHSE(m))
- 신뢰구간이 0을 포함하고 있으면 차이가 유의미하지 않을 수도 있다는 것을 암시함
(2) Duncan's LSR: Least Significant Range
- Tukey에 비해 차이 있다고 결과 나오는 경우가 더 많음
> duncan.test(model, "group", alpha=0.05, console=TRUE ): agricolae 패키지
> LDuncan(model, "group"): laercio 패키지
#텍스트
y1 <- c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8)
y2 <- c(47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9)
y3 <- c(46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2)
y <- c(y1,y2,y3)
n <- rep(10,3)
group <- rep(1:3, n)
#combining into df
group_df <- data.frame(y, group)
sapply(group_df, class)
group_df <- transform(group_df, group =factor(group))
#ANOVA test
aov_model <- aov(y~group, data = group_df)
summary(aov_model)
require("agricolae")
duncan.test(aov_model, "group", alpha = 0.05, console = TRUE)
require("laercio")
LDuncan(aov_model, "group")
(3) 샤페 검정 (Scheffe Test)
- Duncan과 마찬가지로 보수적인 검정
> scheffe.test(model, "vector treatment", alpha=0.05, console=TRUE,
group=TRUE ): agricolae 패키지
각 그룹의 차이가
(검정 결과에 나오는) Minimum Significant Difference보다 작으면
두 그룹은 차이가 없다고 판단함
#텍스트
require("agricolae")
data("sweetpotato")
str(sweetpotato)
attach(sweetpotato)
boxplot(yield~virus,
main = "Yield of sweetpotato - virus",
xlab = "virus", ylab = "yield")
detach(sweetpotato)
require("doBy")
summaryBy(yield~virus, data=sweetpotato,
FUN = c(mean,sd,min,max))
#ANOVA test
aov_sp <- aov(yield~virus, data=sweetpotato)
summary(aov_sp)
#Sheffe test
compa <- scheffe.test(aov_sp,
"virus", group= TRUE,
console = TRUE,
main = "Yield of sp/virus")
3. 이원분산분석 two-way ANOVA
2개의 요인(factors) 내의 요인 수준(factor levels) 간의 조합(combination)의 각각을
개별 집단/처리 (groups, treatment)로 간주하고 이들 간의 평균을 비교하는 것
비용이나 시간 여건이 허락한다면 분석의 신뢰도를 높이기 위해 반복실험을 하여
관측값을 2개 이상 확보하는 것이 좋음
1) 관측값이 하나로만 이루어져 있을 때
aov(y ~ car_type + insurance)
2) 관측값이 두개 이상일 때
aov(score ~ gender.fac + class.fac + gender.fac:class.fac)
여러개의 범주형 변수를 예측변수를 사용하는 다원분산분석을 수행하고 있을 때,
예측변수들 사이에 있을 수 있는 상호작용을 시각적으로 확인하고 싶다면
> interaction.plot(pred1, pred2, 반응변수)
- 선들이 꺾여져 있거나 하면, 비선형성을 의미함
[2] 비모수 검정
1) 두 집단이 독립적인지 짝을 이루는지에 따라
wilcox.test(x~f, data = df) 혹은 wilcox(x, y, mu = 0)
인자: conf.level = / alternative = " " / conf.int =
wilcox.test(x1, x2, paired =TRUE)
2) Kruskal-Wallis 검정 (로버스트 분산분석)
- 분산분석의 비모수적 버전, 데이터의 정규성에 의존하지 않음
- 집단의 중앙값들에 차이가 있는지 알고 싶음
> kruskal.test(y ~ x)
# 분산분석을 통해 모형 비교 (최적의 모형 찾아가기)
> aov(lm1, lm2)
댓글
댓글 쓰기