Rcmdr.kr: An R Commander User in Korea

모델 > 수치적 진단 > 비선형성에 대한 RESET 검정...

Models > Numerical diagnostics > RESET test for nonlinearity...

데이터셋 하나를 만들어보자. x, y1, y2라는 변수를 가진 rawData라는 데이터셋이다.

x <- c(1:30)
y1 <- 1 + x + x^2 + rnorm(30)
y2 <- 1 + x + rnorm(30)
rawData <- as.data.frame(cbind(x, y1, y2))

rawData 데이터셋을 활성화시키자. <활성 데이터셋 없음>을 누르면 활성데이터셋 후보에 rawData가 있을 것이다. 두개의 단순 회귀 모형을 만든다.

RegModel.1 <- lm(y1 ~ x, data=rawData)
summary(RegModel.1)
RegModel.2 <- lm(y2 ~ x, data=rawData)
summary(RegModel.2)

'모델 > 수치적 진단 > 비선형성에 대한 RESET 검정...'을 RegModel.1, RegModel.2에 각각 시행한다. RESET 검정의 하위 메뉴 창이 등장한다. 기본 설정을 그대로 사용하고, 예(OK) 버튼을 누른다.

RESET은 Regression Equation Specification Error Test의 약자이다.

RegModel.1과 RegModel.2의 분포적 특징을 산점도를 활용하여 먼저 시각적으로 살펴보자.

첫째 그래프인 plot(y1, x)의  점들의 분포는 선형이라는 보다는 이차방정식적(quadratic)으로 보인다. 그렇다면, lm(y1 ~ x)라는 선형 회귀 모형을 lm(y1 ~ x + I(x^2))라는 이차방정식적 회귀 모형으로 변화시켜 만들어볼 수 있지 않을까?

RegModel.1_1 <- lm(y1 ~ x + I(x^2), data=rawData)
summary(RegModel.1_1)

RegModel.1_1이 RegModel.1보다 모형의 설명력이 높다고 할 수 있다.

모델 > 수치적 진단 > 변환 반응하기...

Models > Numerical diagnostics > Response transformation...

어느 데이터셋이 활성화되고, 그 데이터셋을 활용한 분석 모형이 만들어지면, '모델 > 수치적 진단 > 변환 반응하기...' 메뉴 창이 활성화된다.

carData 패키지의 Prestige 데이터셋을 이용하여 연습해보자. 먼저 Prestige 데이터셋을 활성화시킨다. '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 통하여 carData 패키지에서 Prestige 데이터셋을 찾아서 선택한다. 다음, Prestige 데이터셋을 활용하여 선형 모델을 만든다. '통계 > 적합성 모델 > 선형 모델...' 메뉴 기능을 선택하여 모형을 만든다.

data(Prestige, package="carData")
LinearModel.1 <- lm(prestige ~ education + income, data=Prestige)
summary(LinearModel.1)

'모델 > 수치적 진단> 변환 반응하기...'메뉴 기능을 선택하면, 다음과 같은 하위 메뉴창이 등장한다. 기본 설정을 이용하자. 예(OK) 버튼을 누른다.

선형 모델을 구성하는 어느 설명변수 또는 설명변수들의 정규성을 향한 변환 기법은 '통계 > 요약 > 정규성을 향해서 변환시키기...' 메뉴 기능을 통해서 연습할 수 있다. 

https://rcmdr.tistory.com/198

한편, 위에 있는 선형 모델, LinearModel.1의 반응변수의 변환 값은 1.2017(Est Power)이다. 이 값을 이용하여 모형을 보다 최적화하는 방법은 아래와 같다.

LM.1_pT <- powerTransform(LinearModel.1, family="bcPower")
LinearModel.1_pT <- lm(bcPower(prestige, LM.1_pT$roundlam) ~ education + income, data=Prestige)
summary(LinearModel.1_pT)   #LinearModel.1의 변환 모형 요약 정보
summary(LinearModel.1)      #LinearModel.1의 요약 정보

bcPower(Prestige$prestige, LM.1_pT$roundlam)  # 변환된 prestige 변수 사례 값
Prestige$prestige                             # 원래의 prestige 변수 사례 값

통계 > 요약 > 정규성을 향해서 변환시키기...

Statistics > Summaries > Transform toward normaltiy...

어느 데이터셋을 활성화시키면, '통계 > 요약 > 정규성을 향해서 변환시키기...' 메뉴 기능이 아래와 같이 활성화될 것이다.

carData 패키지에 포함된 Prestige 데이터셋을 이용하여 연습해보자.

'데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택한후 carData 패키지에 있는 Prestige 데이터셋을 찾아서 선택한다. 그러면, R Commander 상단에 있는 <활성 데이터셋 없음> 이라는 버튼이 'Prestige'로 바뀌면서 '통계 > 요약 > 정규성을 향해서 변환하기...' 메뉴 기능을 이용할 수 있게 된다. 이 메뉴를 선택한다. 그러면, 아래와 같은 화면이 등장한다.  Prestige 데이터셋에서 정규성을 향해서 변환시킬 변수를 선택하는 메뉴 창인데, income 변수를 선택해보자.

car 패키지에 포함된 powerTransform() 함수를 사용하여 출력창에 아래와 같은 결과를 생산한다. 출력내용에서 'Est Power'값으로 제시하는 0.1793을 확인한다. Prestige 패키지에 있는 income 변수에 0.1793 승(제곱)을 하면, 그 사례 값들은 정규 분포적 특성으로 변환될 수 있다는 의미이다.

'통계 > 요약 > 정규성 검정...' 메뉴 기능이 있다. 이 기능을 사용하여, Prestige 데이터셋의 income 변수와 정규성을 향해서 변환된 income 변수의 새로운 변수, 편의상 variable_Income,의 정규성 검정을 시도해보자.

https://rcmdr.kr/88

먼저, 변환된 변수를 Prestige 데이터셋에 추가해야 한다. '데이터 > 활성 데이터셋에 있는 변수 관리하기 > 새로운 변수 계산하기...' 메뉴 기능을 통하여 income 변수가 변환된 variable_Income 변수를 만들어보자. '계산 표현식'에 income^0.1793이 입력된 것을 확인한다. 여기서 0.1793은 powerTransform() 함수를 통하여 앞서 추출한 income 변수의 승(제곱) 값이다.

Prestige 데이터셋에 variable_Income 이라는 새로운 변수가 추가된다. 이제 '통계 > 요약 > 정규성 검정...' 메뉴 기능을 선택하고, 차례로 두번의 정규성 검정을 진행한다. 첫번째는 income 변수로, 두번째는 variable_Income변수를 선택한다. 두개의 정규성 검정 결과를 비교하기 위해서다. 아래 화면은 두번째 검정, variable_Income 변수를 선택하는 메뉴 창이다.

두번의 정규성 검정에서 제시하는 유의 확률값을 살펴보자. variable_Income 변수의 정규성 검정 결과는 p-value=0.07573으로 variable _income 변수는 정규분포적 특징을 갖는다고 할 수 있다.

?powerTransform  # car 패키지의 powerTransform 도움말 보기

# Box Cox Method, univariate
summary(p1 <- powerTransform(cycles ~ len + amp + load, Wool))
# fit linear model with transformed response:
coef(p1, round=TRUE)
summary(m1 <- lm(bcPower(cycles, p1$roundlam) ~ len + amp + load, Wool))

# Multivariate Box Cox uses Highway1 data
summary(powerTransform(cbind(len, adt, trks, sigs1) ~ 1, Highway1))

# Multivariate transformation to normality within levels of 'htype'
summary(a3 <- powerTransform(cbind(len, adt, trks, sigs1) ~ htype, Highway1))

# test lambda = (0 0 0 -1)
testTransform(a3, c(0, 0, 0, -1))

# save the rounded transformed values, plot them with a separate
# color for each highway type
transformedY <- bcPower(with(Highway1, cbind(len, adt, trks, sigs1)),
                        coef(a3, round=TRUE))
## Not run: scatterplotMatrix( ~ transformedY|htype, Highway1) 

# With negative responses, use the bcnPower family
m2 <- lm(I1L1 ~ pool, LoBD)
summary(p2 <- powerTransform(m2, family="bcnPower"))
testTransform(p2, .5)
summary(powerTransform(update(m2, cbind(LoBD$I1L2, LoBD$I1L1) ~ .), family="bcnPower"))

## Not run:  
  # takes a few seconds:
  # multivariate bcnPower, with 8 responses
  summary(powerTransform(update(m2, as.matrix(LoBD[, -1]) ~ .), family="bcnPower"))
  # multivariate bcnPower, fit with one iteration using starting values as estimates
  summary(powerTransform(update(m2, as.matrix(LoBD[, -1]) ~ .), family="bcnPower", itmax=1))

## End(Not run)

# mixed effects model
## Not run: 
  # uses the lme4 package
  data <- reshape(LoBD[1:20, ], varying=names(LoBD)[-1], direction="long", v.names="y")
  names(data) <- c("pool", "assay", "y", "id")
  data$assay <- factor(data$assay)
  require(lme4)
  m2 <- lmer(y ~ pool + (1|assay), data)
  summary(l2 <- powerTransform(m2, family="bcnPower", verbose=TRUE))

## End(Not run)