Rcmdr.kr: An R Commander User in Korea

통계 > 적합성 모델 > 서열(Ordinal) 회귀 모델...
Statistics > Fit models > Ordinal regression model...

MASS 패키지에 있는 housing 데이터셋을 활용해보자. 먼저 housing 데이터셋을 활성화 시킨다. https://rcmdr.tistory.com/215

Sat는 거주자의 현 거주환경에 대한 만족도에 관한 변수로서, High > Medium > Low 라는 서열화된 요인을 갖고 있다. Sat를 반응변수로, 나머지를 설명변수의 후보군으로 모델을 구성하는 것이다. Freq 변수는 Sat, Infl, Type, Cont 등으로 구별되는 집단 구성원의 숫자를 뜻한다. 아래와 같이 모형을 만든다.

OrdRegModel.1 <- polr(Sat ~ Infl + Type + Cont , weights = Freq, method="logistic", 
  data=housing, Hess=TRUE)
summary(OrdRegModel.1)

OrdRegModel.1 이라는 모형이 만들어졌다면, '모델 > 가설 검정 > 분산분석표...' 메뉴 기능을 이용할 수 있다.

MASS::housing()

library(MASS, pos=16)
data(housing, package="MASS")

'도구 > 패키지 적재하기...' 메뉴 기능을 선택하고, MASS 패키지를 찾아서 선택한다.

그리고나서, '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하면 하위 선택 창으로 이동한다. 아래와 같이 MASS 패키지를 선택하고, housing 데이터셋을 찾아 선택한다.

housing 데이터셋이 활성화된다. R Commander 상단의 메뉴에서 < 활성 데이터셋 없음> 이 'housing'로 바뀐다.

summary(housing)
str(housing)

'통계 > 요약 > 활성 데이터셋' 메뉴 기능을 선택하여 housing 데이터셋의 요약 정보를 살펴보자. 아울러 입력창에 str(housing)을 입력하고 <실행하기> 버튼을 누르자.

데이터셋의 내부는 다음과 같다:

Frequency Table from a Copenhagen Housing Conditions Survey

Description

The housing data frame has 72 rows and 5 variables.

Usage

housing

Format

Sat

Satisfaction of householders with their present housing circumstances, (High, Medium or Low, ordered factor).

Infl

Perceived degree of influence householders have on the management of the property (High, Medium, Low).

Type

Type of rental accommodation, (Tower, Atrium, Apartment, Terrace).

Cont

Contact residents are afforded with other residents, (Low, High).

Freq

Frequencies: the numbers of residents in each class.

Source

Madsen, M. (1976) Statistical analysis of multiple contingency tables. Two examples. Scand. J. Statist. 3, 97–106.

Cox, D. R. and Snell, E. J. (1984) Applied Statistics, Principles and Examples. Chapman & Hall.

References

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth edition. Springer.

Examples

options(contrasts = c("contr.treatment", "contr.poly"))

# Surrogate Poisson models
house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family = poisson,
                  data = housing)
## IGNORE_RDIFF_BEGIN
summary(house.glm0, cor = FALSE)
## IGNORE_RDIFF_END

addterm(house.glm0, ~. + Sat:(Infl+Type+Cont), test = "Chisq")

house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont))
summary(house.glm1, cor = FALSE)

1 - pchisq(deviance(house.glm1), house.glm1$df.residual)

dropterm(house.glm1, test = "Chisq")

addterm(house.glm1, ~. + Sat:(Infl+Type+Cont)^2, test  =  "Chisq")

hnames <- lapply(housing[, -5], levels) # omit Freq
newData <- expand.grid(hnames)
newData$Sat <- ordered(newData$Sat)
house.pm <- predict(house.glm1, newData,
                    type = "response")  # poisson means
house.pm <- matrix(house.pm, ncol = 3, byrow = TRUE,
                   dimnames = list(NULL, hnames[[1]]))
house.pr <- house.pm/drop(house.pm %*% rep(1, 3))
cbind(expand.grid(hnames[-1]), round(house.pr, 2))

# Iterative proportional scaling
loglm(Freq ~ Infl*Type*Cont + Sat*(Infl+Type+Cont), data = housing)


# multinomial model
library(nnet)
(house.mult<- multinom(Sat ~ Infl + Type + Cont, weights = Freq,
                       data = housing))
house.mult2 <- multinom(Sat ~ Infl*Type*Cont, weights = Freq,
                        data = housing)
anova(house.mult, house.mult2)

house.pm <- predict(house.mult, expand.grid(hnames[-1]), type = "probs")
cbind(expand.grid(hnames[-1]), round(house.pm, 2))

# proportional odds model
house.cpr <- apply(house.pr, 1, cumsum)
logit <- function(x) log(x/(1-x))
house.ld <- logit(house.cpr[2, ]) - logit(house.cpr[1, ])
(ratio <- sort(drop(house.ld)))
mean(ratio)

(house.plr <- polr(Sat ~ Infl + Type + Cont,
                   data = housing, weights = Freq))

house.pr1 <- predict(house.plr, expand.grid(hnames[-1]), type = "probs")
cbind(expand.grid(hnames[-1]), round(house.pr1, 2))

Fr <- matrix(housing$Freq, ncol  =  3, byrow = TRUE)
2*sum(Fr*log(house.pr/house.pr1))

house.plr2 <- stepAIC(house.plr, ~.^2)
house.plr2$anova

MASS 패키지에는 birthwt라는 데이터셋이 포함되어 있다. birthwt 데이터셋을 활용하여 bwt라는 2차 데이터셋이 만들어진다. 

bwt <- with(birthwt, {
race <- factor(race, labels = c("white", "black", "other"))
ptd <- factor(ptl > 0)
ftv <- factor(ftv)
levels(ftv)[-(1:2)] <- "2+"
data.frame(low = factor(low), age, lwt, race, smoke = (smoke > 0),
           ptd, ht = (ht > 0), ui = (ui > 0), ftv)
})

bwt <- with(birthwt, {
race <- factor(race, labels = c("white", "black", "other"))
ptd <- factor(ptl > 0)
ftv <- factor(ftv)
levels(ftv)[-(1:2)] <- "2+"
data.frame(low = factor(low), age, lwt, race, smoke = (smoke > 0),
           ptd, ht = (ht > 0), ui = (ui > 0), ftv)
})

bwt 데이터셋은 분석 모형을 만드는데 간혹 예제로 사용되는데, birthwt에서 bwt가 만들어지는 과정이 R Commander 기본 사용자에게는 다소 어렵게 느껴질수 있겠다는 판단이다. 데이터셋 자체에 대한 이해의 어려움 때문에 분석 모형의 구성과 해석으로 나아가지 못하는 경우가 있어, bwt 데이터셋 설명을 하고자 한다.

bwt 데이터셋은 저체중아 출생의 원인을 찾고자 하는 문제의식을 담고 있다. low 변수는 출생당시 몸무게가 2.5kg 미만 여부를 담고 있으며, 반응변수가 된다. 나머지 변수들은 저체중아 출산에 영향을 끼치는가 여부인 설명변수들의 후보군이 되겠다.

options(contrasts = c("contr.treatment", "contr.poly"))
GLM.1 <- glm(low ~ ., binomial, bwt)

통계 > 적합성 모델 > 선형 혼합 모델...
Statistics > Fit models > Linear mixed model...

데이터셋을 활성화시키면, '통계 > 적합성 모델 > 선형 혼합 모델...' 메뉴 기능을 사용할 수 있다. lme4 패키지의 sleepstudy 데이터셋을 이용하여 연습해보자.

sleepstudy 데이터셋을 활성화 시키자. 먼저 lme4 패키지를 호출해야 한다. 그래야 포함된 데이터셋 목록을 확인할 수 있기 때문이다. '도구 > 패키지 적재하기...' 메뉴 기능을 통하여 lme4를 적재한다. 그리고 '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 lme4 패키지에 포함된 데이터셋들 중에서 sleepstudy를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'sleepstudy'로 바뀐다.

https://rcmdr.tistory.com/212

sleepstudy 데이터셋을 이용하여 LMM.1 모형을 만든다. '통계 > 적합성 모델 > 선형 혼합 모델...' 메뉴 기능을 이용할 수 있다.

LMM.1 <- lmer(Reaction ~ Days + (Days | Subject), data=sleepstudy, REML=TRUE)

REML(Restricted maximum likelihood, 제한적 최대우도) 대신 ML(Maximum likelihood, 최대우도) 인자를 사용해보자.

LMM.2 <- lmer(Reaction ~ Days + (Days | Subject), data=sleepstudy, REML=FALSE)

?lmer   # lmer 함수 도움말 보기

선형 혼합 모형(Linear mixed model)을 만들었다면, '모델 > 가설 검정 > 분산분석표...' 메뉴 기능을 이용할 수 있다. LMM.1을 만들었으니 가능하다.

Anova(LMM.1, type="II")

lme4::sleepstudy()

data(sleepstudy, package="lme4")

'데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하면 하위 선택 창으로 이동한다. 아래와 같이 lme4 패키지를 선택하고, sleepstudy 데이터셋을 찾아 선택한다.

sleepstudy 데이터셋이 활성화된다. R Commander 상단의 메뉴에서 < 활성 데이터셋 없음> 이 'sleepstudy'로 바뀐다.

summary(sleepstudy)
str(sleepstudy)

'통계 > 요약 > 활성 데이터셋' 메뉴 기능을 통해서 sleepstudy 데이터의 요약 정보를 살펴보자. str() 함수를 이용하여 sleepstudy 데이터셋의 내부 구조를 살펴보자.

데이터셋의 내부는 다음과 같다:

Reaction times in a sleep deprivation study

Description

The average reaction time per day for subjects in a sleep deprivation study. On day 0 the subjects had their normal amount of sleep. Starting that night they were restricted to 3 hours of sleep per night. The observations represent the average reaction time on a series of tests given each day to each subject.

Format

A data frame with 180 observations on the following 3 variables.

Reaction

Average reaction time (ms)

Days

Number of days of sleep deprivation

Subject

Subject number on which the observation was made.

Details

These data are from the study described in Belenky et al. (2003), for the sleep-deprived group and for the first 10 days of the study, up to the recovery period.

References

Gregory Belenky, Nancy J. Wesensten, David R. Thorne, Maria L. Thomas, Helen C. Sing, Daniel P. Redmond, Michael B. Russo and Thomas J. Balkin (2003) Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: a sleep dose-response study. Journal of Sleep Research 12, 1–12.

Examples

str(sleepstudy)
require(lattice)
xyplot(Reaction ~ Days | Subject, sleepstudy, type = c("g","p","r"),
       index = function(x,y) coef(lm(y ~ x))[1],
       xlab = "Days of sleep deprivation",
       ylab = "Average reaction time (ms)", aspect = "xy")
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy))

통계 > 평균 > 일-요인 반복-측정 ANOVA/ANCOVA...

Statistics > Means > One-factor repeated-measures ANOVA/ANCOVA...

데이터셋을 활성화시켰다면, '통계 > 평균 > 일-요인 반복-측정 ANOVA/ANCOVA' 메뉴 기능을 사용할 수 있다. carData 패키지의  OBrienKaiser 데이터셋을 이용하여 연습해보자.

먼저, OBrienKaiser 데이터셋을 활성화 시키자. '데이터 > 패키지에 있는 데이터 >

첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 carData 패키지에 포함된 데이터셋들 중에서 OBrienKasier를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'OBrienKasier'로 바뀐다.

https://rcmdr.tistory.com/95

'통계 > 평균 > 일-요인 반복-측정 ANOVA/ANCOVA' 메뉴 기능을 선택하자. 아래와 같은 화면이 등장한다.

메뉴 창의 세부 항목을 아래와 같이 설정한다.

data(OBrienKaiser, package="carData")
Anova(lm(cbind(pre.1, post.1, fup.1) ~ gender *treatment, data=OBrienKaiser), 
  idata=data.frame(Phase=factor(c('pretest', 'posttest', 'followup'))), idesign = ~Phase, 
  test.statistic="Pillai")

'선택기능' 메뉴 창의 중간에 있는 '반복-측정 평균들'에서 '요인 반응 평균 그리기를 선택하자. 그래픽 창이 등장하며 세개의 단계 'pretest-posttest-followup'의 각 1 level에 해당하는 평균 값 그래프를 생산한다.

repeatedMeasuresPlot(OBrienKaiser, within=c("pre.1", "post.1", "fup.1"), within.names="Phase", 
  within.levels=list(Phase=c("pretest", "posttest", "followup")), print.tables=FALSE, plot.means=TRUE, 
  response.name="score")

만약 '개인간 요인 (없거나 하나 이상 선택)'에서 gender 와 treatment 두개를 모두 선택하고 그래프를 출력하면 아래와 같이 gender와 treatment변로 평균 그림이 생산된다.

foreign 패키지 안에는 SPSS 파일이 예제로 포함되어 있다.

R Commander에서 SPSS 파일 불러오기 연습은 아래를 참조할 수 있다:

https://rcmdr.tistory.com/30

모델 > 그래프 > 기초 진단 그림

Models > Graphs > Basic diagnostic plots

데이터셋을 활성화시키고, 분석 모형을 만들었다면, '모델 > 그래프 > 기초 진단 그림...' 메뉴 기능을 사용할 수 있다. carData 패키지의 Prestige 데이터셋을 이용하여 연습해보자.

먼저, Prestige 데이터셋을 활성화 시키자. '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 carData 패키지에 포함된 데이터셋들 중에서 Prestige를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'Prestige'로 바뀐다.

Prestige 데이터셋을 이용하여 LinearModel.1 모형을 만든다. '통계 > 적합성 모델 > 선형 모델...' 메뉴 기능을 이용할 수 있다. 

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

LinearModel.1 모형이 활성화된 이후, '모델 > 그래프 > 기초 진단 그림...' 메뉴 기능을 선택하면 아래와 같은 그래픽 장치 창이 등장한다.

oldpar <- par(oma=c(0,0,3,0), mfrow=c(2,2))
plot(LinearModel.1)
par(oldpar)

모델 > 그래프 > Predictor 효과 그림...
Models > Graphs > Predictor effect plots...

데이터셋을 활성화시키고, 분석 모형을 만들었다면, '모델 > 그래프 > Predictor 효과 그림...' 메뉴 기능을 사용할 수 있다. carData 패키지의 Prestige 데이터셋을 이용하여 연습해보자.

먼저, Prestige 데이터셋을 활성화 시키자. '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 carData 패키지에 포함된 데이터셋들 중에서 Prestige를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'Prestige'로 바뀐다.

Prestige 데이터셋을 이용하여 LinearModel.1 모형을 만든다. '통계 > 적합성 모델 > 선형 모델...' 메뉴 기능을 이용할 수 있다.

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

LinearModel.1 모형이 활성화된 이후, '모델 > 그래프 > Predictor 효과 그림...' 메뉴 기능을 선택하면 아래와 같은 하위 화면이 등장한다. 예측변수 두개를 모두 선택한다. 메뉴 창 오른쪽 끝 아래에 있는 예(OK) 버튼을 누른다.

plot(predictorEffects(LinearModel.1))

?predictorEffects   #  Predictor 효과 그림 도움말보기

모델 > 그래프 > 영향력 색인 그림...

Models > Graphs > Influence index plot...

데이터셋을 활성화시키고, 분석 모형을 만들었다면, '모델 > 그래프 > 영향력 색인 그림...' 메뉴 기능을 사용할 수 있다. carData 패키지의 Prestige 데이터셋을 이용하여 연습해보자.

먼저, Prestige 데이터셋을 활성화 시키자. '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 carData 패키지에 포함된 데이터셋들 중에서 Prestige를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'Prestige'로 바뀐다.

Prestige 데이터셋을 이용하여 LinearModel.1 모형을 만든다. '통계 > 적합성 모델 > 선형 모델...' 메뉴 기능을 이용할 수 있다. 

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

LinearModel.1 모형이 활성화된 이후, '모델 > 그래프 > 영향력 색인 그림...' 메뉴 기능을 선택하면 아래와 같은 하위 화면이 등장한다. 메뉴 창 오른쪽 끝 아래에 있는 예(OK) 버튼을 누른다.

influenceIndexPlot(LinearModel.1, id=list(method="y", n=2), vars=c("Cook", "Studentized", "Bonf", 
  "hat"))

?influenceIndexPlot   # 영향력 색인 그림의 도움말 보기