부록 A: 샘플 세션
부록 A: 샘플 세션 (A Sample Session)
아래 세션은 R 환경의 여러 기능을 직접 실행해 보면서 익숙해지도록 준비된 예제예요. 처음에는 낯설고 헷갈리는 기능이 많겠지만, 이 헷갈림은 곧 사라질 거예요. 하나씩 명령을 따라 치다 보면 R이 어떻게 생겼는지 손으로 느낄 수 있을 거예요.
출처: R 공식 매뉴얼
본문
먼저 여러분 플랫폼에 맞게 R을 실행해 볼게요 (실행 방법은 Invoking R을 참고하세요). R 프로그램이 시작되면 배너가 나타나요. (여기서 R 코드를 보여줄 때는 왼쪽의 프롬프트를 생략할 거예요 — 혼동을 피하기 위해서요.)
도움말과 첫 객체 만들기
help.start() 명령으로 R 온라인 도움말의 HTML 인터페이스를 시작해 볼게요. 여러분 컴퓨터의 웹 브라우저가 열리면서 도움말 창이 뜰 거예요. 마우스로 이 기능을 잠깐 살펴보고, 도움말 창은 아이콘으로 최소화한 뒤 다음 단계로 넘어가요.
help.start()
이제 첫 계산 객체를 만들어 볼게요. 정규분포에서 난수를 두 벡터로 뽑아 x, y 좌표로 삼는 거예요.
x <- rnorm(50)
y <- rnorm(x)
위 코드는 표준정규분포에서 난수 50개를 뽑아 x에 넣고, 다시 x와 같은 길이의 난수를 y에 넣어요. 이렇게 생성된 의사난수(pseudo-random) 정규 벡터 두 개가 x, y 좌표가 됩니다.
plot(x, y)
plot()은 점들을 평면에 그려 줘요. 그래픽 창이 자동으로 열리면서 산점도가 나타날 거예요.
ls()
ls()로 현재 R 워크스페이스에 어떤 객체들이 있는지 확인할 수 있어요. 방금 만든 x와 y가 보일 거예요.
rm(x, y)
rm()은 더 이상 필요 없는 객체를 지우는 함수예요. 정리(clean up)라고 생각하면 돼요. x, y를 삭제했어요.
데이터 프레임과 단순 선형 회귀
이번에는 데이터를 만들고 회귀 분석까지 해 볼게요. 먼저 1부터 20까지의 정수 벡터를 만들게요.
x <- 1:20
이제 x = (1, 2, ..., 20)이 됐어요. 이 값들에 가중치를 줄 표준편차 벡터를 만들어 볼게요.
w <- 1 + sqrt(x)/2
w는 표준편차를 담은 '가중치' 벡터예요. 이제 x와 y 두 열로 된 데이터 프레임을 만들고 확인해 볼게요.
dummy <- data.frame(x=x, y= x + rnorm(x)*w)
dummy
data.frame()으로 두 열(x, y)을 가진 데이터 프레임을 만들고, dummy라고 부르기로 했어요. 그 아래 dummy를 입력해 내용을 눈으로 확인해요.
이제 이 데이터로 단순 선형 회귀를 적합시켜 보고 그 분석 결과를 살펴볼게요.
fm <- lm(y ~ x, data=dummy)
summary(fm)
lm()은 선형 모델(linear model)을 적합시키는 함수예요. 수식 y ~ x에서 물결표(~) 오른쪽의 x가 설명 변수, 왼쪽의 y가 반응 변수예요. 즉 y가 x에 의존한다고 모델링하는 거죠. summary()로 분석 결과를 요약해 볼 수 있어요.
우리는 실제 표준편차를 알고 있으니, 가중 회귀(weighted regression)도 시도해 볼 수 있어요.
fm1 <- lm(y ~ x, data=dummy, weight=1/w^2)
summary(fm1)
표준편차를 이미 알고 있으니 그 역수의 제곱을 가중치로 써서 회귀를 다시 적합시킨 거예요.
데이터 프레임 다루기와 비모수 회귀
attach()를 쓰면 데이터 프레임의 열을 변수처럼 바로 쓸 수 있어요.
attach(dummy)
attach(dummy)로 데이터 프레임의 열들을 변수로 바로 꺼내 쓸 수 있게 됐어요. 이제 비모수(nonparametric) 국소 회귀 함수를 만들어 볼게요.
lrf <- lowess(x, y)
lowess()로 분포에 대한 가정이 없는(비모수) 국소 회귀 함수를 만들었어요. 이제 그래프로 결과를 비교해서 볼게요.
plot(x, y)
lines(x, lrf$y)
plot()으로 표준 산점도를 그리고, lines()로 국소 회귀 결과(lrf$y)를 선으로 더해 줬어요.
이제 회귀선들을 한 화면에 겹쳐 그려 볼게요. 좌표는 abline()으로 그릴 수 있어요.
abline(0, 1, lty=3)
abline(coef(fm))
abline(coef(fm1), col = "red")
abline(0, 1, lty=3)은 절편 0, 기울기 1인 참 회귀선(점선)이에요.abline(coef(fm))은 가중치를 쓰지 않은 회귀선이에요.abline(coef(fm1), col = "red")은 가중 회귀선(빨간색)이에요.
이제 검색 경로에서 데이터 프레임을 빼낼게요.
detach()
detach()로 데이터 프레임을 검색 경로에서 제거했어요. 다음으로 회귀 진단 플롯을 그려 볼게요.
plot(fitted(fm), resid(fm),
xlab="Fitted values",
ylab="Residuals",
main="Residuals vs Fitted")
적합값(fitted) 대 잔차(residual)를 그린 표준 회귀 진단 플롯이에요. **이분산성(heteroscedasticity)**이 있는지 확인하는 용도인데, 여기서 보이나요? 분산이 고르게 퍼져 보이면 문제없는 모양이에요.
qqnorm(resid(fm), main="Residuals Rankit Plot")
qqnorm()은 정규 점수 플롯으로, 잔차의 **왜도(skewness), 첨도(kurtosis), 이상치(outliers)**를 확인해 주는 그래프예요. (여기 데이터에서는 크게 유용하진 않아요.)
rm(fm, fm1, lrf, x, dummy)
분석이 끝났으니 다시 정리해 줄게요.
Michelson 데이터 읽기와 분산 분석
다음 절에서는 Michelson이 빛의 속도를 측정한 고전 실험 데이터를 다뤄 볼게요. 이 데이터는 morley 객체에 이미 들어 있지만, read.table() 함수를 설명하기 위해 파일에서 직접 읽어 볼게요.
filepath <- system.file("data", "morley.tab" , package="datasets")
filepath
system.file()로 데이터 파일의 경로를 얻었어요. (선택사항) file.show(filepath)로 파일 내용을 눈으로 확인해 볼 수도 있어요.
mm <- read.table(filepath)
mm
Michelson 데이터를 데이터 프레임으로 읽어 mm에 저장하고 확인했어요. 이 데이터에는 5개의 실험(Expt 열)이 있고, 각각 20번의 실행(Run 열)이 있어요. sl은 적절히 부호화된 빛의 속도 측정값이에요.
Expt와 Run을 요인(factor)으로 바꿔 볼게요.
mm$Expt <- factor(mm$Expt)
mm$Run <- factor(mm$Run)
실험 번호와 실행 번호를 숫자가 아니라 범주형 **요인(factor)**으로 바꾼 거예요. 이제 데이터 프레임을 검색 경로에 올리고, 실험들을 상자 그림으로 비교해 볼게요.
attach(mm)
plot(Expt, Speed, main="Speed of Light Data", xlab="Experiment No.")
attach(mm)로 데이터 프레임을 기본 위치(2번)에 올렸고, plot()으로 다섯 실험의 단순 상자 그림(boxplot)을 비교했어요.
이제 무작위 블록(randomized block) 설계로 분산 분석을 해 볼게요.
fm <- aov(Speed ~ Run + Expt, data=mm)
summary(fm)
aov()로 분산 분석을 수행했어요. Run(실행)과 Expt(실험)를 요인으로 두고, Speed를 반응 변수로 분석한 거예요.
이제 Run을 뺀 부분 모형을 적합시켜 공식적인 분산 분석으로 비교해 볼게요.
fm0 <- update(fm, . ~ . - Run)
anova(fm0, fm)
update()로 Run 항을 제외한 부분 모형(fm0)을 만들고, anova()로 두 모형을 공식적으로 비교했어요.
detach()
rm(fm, fm0)
다음으로 넘어가기 전에 다시 정리해 줄게요.
등고선과 이미지 플롯
이번에는 좀 더 그래픽적인 기능인 등고선(contour)과 이미지(image) 플롯을 살펴볼게요.
x <- seq(-pi, pi, len=50)
y <- x
x는 구간 [-pi, pi]에 등간격으로 50개 들어 있는 벡터예요. y도 같게 만들었어요.
f <- outer(x, y, function(x, y) cos(y)/(1 + x^2))
outer()로 함수 cos(y)/(1 + x^2)의 값을 행과 열로 담은 정사각 행렬 f를 만들었어요. 행은 x, 열은 y로 인덱스됩니다.
oldpar <- par(no.readonly = TRUE)
par(pty="s")
그래픽 파라미터를 저장하고(oldpar), 플롯 영역을 "정사각형(square)"으로 설정했어요. 나중에 원래대로 되돌리기 위해 저장해 둔 거예요.
contour(x, y, f)
contour(x, y, f, nlevels=15, add=TRUE)
f의 등고선 지도를 그리고, 더 자세하게 볼 수 있도록 등고선을 추가로 더 그렸어요.
fa <- (f-t(f))/2
fa는 f의 "비대칭 부분(asymmetric part)"이에요. (t()는 전치(transpose) 함수예요.)
contour(x, y, fa, nlevels=15)
fa의 등고선 플롯을 그렸어요.
par(oldpar)
그래픽 파라미터를 원래대로 복원했어요.
image(x, y, f)
image(x, y, fa)
고밀도 이미지 플롯을 만들었어요. 필요하면 이 화면의 하드카피(하드카피 인쇄)를 뽑을 수도 있어요.
objects(); rm(x, y, f, fa)
objects()로 현재 목록을 확인한 뒤, 다음 단계로 넘어가기 전에 정리했어요.
복소수 연산
R은 복소수 연산도 지원해요.
th <- seq(-pi, pi, len=100)
z <- exp(1i*th)
1i는 복소수 i를 나타내는 표현이에요. exp(1i*th)로 단위원 위의 복소수들을 만들어 z에 담았어요.
par(pty="s")
plot(z, type="l")
복소수 인자를 플롯하면 허수부를 실수부에 대해 그리게 돼요. 원 그래프가 나오는 게 정상이에요.
w <- rnorm(100) + rnorm(100)*1i
단위원 안에서 점을 샘플링하고 싶다면, 한 방법은 실수부·허수부가 각각 표준정규인 복소수를 뽑는 거예요.
w <- ifelse(Mod(w) > 1, 1/w, w)
Mod()는 복소수의 크기(절댓값)예요. 단위원 바깥에 있는 점은 역수(1/w)로 사상해 안쪽으로 옮겨 넣었어요.
plot(w, xlim=c(-1,1), ylim=c(-1,1), pch="+",xlab="x", ylab="y")
lines(z)
모든 점이 단위원 안에 들어왔지만, 분포가 균등하지는 않아요. 두 번째 방법은 균등분포를 쓰는 거예요.
w <- sqrt(runif(100))*exp(2*pi*runif(100)*1i)
plot(w, xlim=c(-1,1), ylim=c(-1,1), pch="+", xlab="x", ylab="y")
lines(z)
두 번째 방법은 균등분포(runif)를 사용해요. 이번에는 점들이 원반 위에 훨씬 고르게 퍼져 보일 거예요.
rm(th, w, z)
다시 정리해 줄게요.
세션 종료
q()
q()는 R 프로그램을 종료하는 함수예요. R 워크스페이스를 저장할지 물어보는데, 지금처럼 탐구적인 세션에서는 저장하지 않는 편이 좋아요. 이렇게 명령 하나씩 따라 하면서 R의 핵심 객체, 데이터 프레임, 회귀·분산분석, 그래픽, 복소수 연산까지 한 바퀴 돌아봤어요.
더 알아보기
- 도움말 시스템:
help.start()로 시작하는 온라인 도움말과help(함수명)형태의 개별 도움말 - 데이터 프레임:
data.frame(),attach()/detach(),read.table()로 데이터를 다루는 기본기 - 모델링:
lm()(선형 회귀),aov()(분산 분석),update()·anova()를 이용한 모형 비교 - 그래픽:
plot(),lines(),abline(),contour(),image()등 기본 플로팅 함수들 - 이후 챕터에서는 이번에 나온 함수들을 더 깊이, 그리고 실제 데이터 분석에 어떤 의미가 있는지 자세히 다뤄요.