통계 > 적합성 모델 > 일반화 선형 혼합 모델... Statistics > Fit models > Generalized linear mixed model...
'도구 > 패키지 적재하기...' 메뉴 기능을 이용하여 lme4 패키지를 찾아서 적재하자. lme4 패키지에는 일반화 선형 혼합 모델을 만들고 분석하는데 필요한 glmer()와 예제 데이터셋 cbpp가 포함되어 있다.
'데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 통하여 lme4 패키지에 있는 cbpp 데이터셋을 찾아서 선택하자. 그러면 R Commander의 상단에 있는 <활성 데이터셋 없음>이 'cbpp'로 활성화될 것이다. https://rcmdr.tistory.com/240
require(lme4)
data(cbpp, package="lme4")
GLMM.1 <- glmer(incidence / size ~ period + (1 | herd ), family=binomial(logit), data=cbpp,
weights=size)
summary(GLMM.1)
exp(coef(GLMM.1)) # Exponentiated coefficients ("odds ratios")
Anova(GLMM.1) # period 변수의 영향 여부 검정
요인형 변수를 두개 이상 가지고 있는 데이터셋이 활성화되어 있다면, '통계 > 비율 > 이-표본 비율 검정..' 메뉴 기능을 이용할 수 있다. carData 패키지에 있는 Chile 데이터셋을 활용해서 연습해보자. 먼저, '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 통하여 Chile 데이터셋을 활성화시키자. R Commander 상단에 'Chile'라는 데이터셋이 활성화되었는지 확인하자.
Statistics > Means > Two-factor repeated-measures ANOVA/ANCOVa...
데이터셋을 활성화시켰다면, '통계 > 평균 > 일-요인 반복-측정 ANOVA/ANCOVA' 메뉴 기능을 사용할 수 있다. carData 패키지의 OBrienKaiser 데이터셋을 이용하여 연습해보자.
먼저, OBrienKaiser 데이터셋을 활성화 시키자. '데이터 > 패키지에 있는 데이터 >
첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 carData 패키지에 포함된 데이터셋들 중에서 OBrienKasier를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'OBrienKasier'로 바뀐다.
Statistics > Contingency tables > Enter and analyze two-way table...
'통계 > 분할표 > 이원표 입력 및 분석하기...' 메뉴 기능을 선택하면 하위 창이 등장한다. '변수 이름', '행과 열의 수', '사례 수' 등을 입력할 수 있다. 아래의 내용은 'chisq.test' 함수 도움말 문서에 나오는 사례를 입력한다. 아래의 입력 스크립트를 참조 할 수 있다.
?chisq.test # stats 패키지의 chisq.test 도움말 보기
## From Agresti(2007) p.39
M <- as.table(rbind(c(762, 327, 468), c(484, 239, 477)))
dimnames(M) <- list(gender = c("F", "M"),
party = c("Democrat","Independent", "Republican"))
(Xsq <- chisq.test(M)) # Prints test summary
Xsq$observed # observed counts (same as M)
Xsq$expected # expected counts under the null
Xsq$residuals # Pearson residuals
Xsq$stdres # standardized residuals
## Effect of simulating p-values
x <- matrix(c(12, 5, 7, 7), ncol = 2)
chisq.test(x)$p.value # 0.4233
chisq.test(x, simulate.p.value = TRUE, B = 10000)$p.value
# around 0.29!
## Testing for population probabilities
## Case A. Tabulated data
x <- c(A = 20, B = 15, C = 25)
chisq.test(x)
chisq.test(as.table(x)) # the same
x <- c(89,37,30,28,2)
p <- c(40,20,20,15,5)
try(
chisq.test(x, p = p) # gives an error
)
chisq.test(x, p = p, rescale.p = TRUE)
# works
p <- c(0.40,0.20,0.20,0.19,0.01)
# Expected count in category 5
# is 1.86 < 5 ==> chi square approx.
chisq.test(x, p = p) # maybe doubtful, but is ok!
chisq.test(x, p = p, simulate.p.value = TRUE)
## Case B. Raw data
x <- trunc(5 * runif(100))
chisq.test(table(x)) # NOT 'chisq.test(x)'!
Sat는 거주자의 현 거주환경에 대한 만족도에 관한 변수로서, High > Medium > Low 라는 서열화된 요인을 갖고 있다. Sat를 반응변수로, 나머지를 설명변수의 후보군으로 모델을 구성하는 것이다. Freq 변수는 Sat, Infl, Type, Cont 등으로 구별되는 집단 구성원의 숫자를 뜻한다. 아래와 같이 모형을 만든다.
통계 > 적합성 모델 > 선형 혼합 모델... Statistics > Fit models > Linear mixed model...
데이터셋을 활성화시키면, '통계 > 적합성 모델 > 선형 혼합 모델...' 메뉴 기능을 사용할 수 있다. lme4 패키지의 sleepstudy 데이터셋을 이용하여 연습해보자.
sleepstudy 데이터셋을 활성화 시키자. 먼저 lme4 패키지를 호출해야 한다. 그래야 포함된 데이터셋 목록을 확인할 수 있기 때문이다. '도구 > 패키지 적재하기...' 메뉴 기능을 통하여 lme4를 적재한다. 그리고 '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 lme4 패키지에 포함된 데이터셋들 중에서 sleepstudy를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'sleepstudy'로 바뀐다.
Statistics > Means > One-factor repeated-measures ANOVA/ANCOVA...
데이터셋을 활성화시켰다면, '통계 > 평균 > 일-요인 반복-측정 ANOVA/ANCOVA' 메뉴 기능을 사용할 수 있다. carData 패키지의 OBrienKaiser 데이터셋을 이용하여 연습해보자.
먼저, OBrienKaiser 데이터셋을 활성화 시키자. '데이터 > 패키지에 있는 데이터 >
첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택하고, 다음 화면에서 carData 패키지에 포함된 데이터셋들 중에서 OBrienKasier를 찾아서 선택한다. 그러면, R Commander 상단의 <활성 데이터셋 없음> 버튼이 'OBrienKasier'로 바뀐다.
어느 데이터셋을 활성화시키면, '통계 > 요약 > 정규성을 향해서 변환시키기...' 메뉴 기능이 아래와 같이 활성화될 것이다.
carData 패키지에 포함된 Prestige 데이터셋을 이용하여 연습해보자.
'데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...' 메뉴 기능을 선택한후 carData 패키지에 있는 Prestige 데이터셋을 찾아서 선택한다. 그러면, R Commander 상단에 있는 <활성 데이터셋 없음> 이라는 버튼이 'Prestige'로 바뀌면서 '통계 > 요약 > 정규성을 향해서 변환하기...' 메뉴 기능을 이용할 수 있게 된다. 이 메뉴를 선택한다. 그러면, 아래와 같은 화면이 등장한다. Prestige 데이터셋에서 정규성을 향해서 변환시킬 변수를 선택하는 메뉴 창인데, income 변수를 선택해보자.
car 패키지에 포함된 powerTransform() 함수를 사용하여 출력창에 아래와 같은 결과를 생산한다. 출력내용에서 'Est Power'값으로 제시하는 0.1793을 확인한다. Prestige 패키지에 있는 income 변수에 0.1793 승(제곱)을 하면, 그 사례 값들은 정규 분포적 특성으로 변환될 수 있다는 의미이다.
'통계 > 요약 > 정규성 검정...' 메뉴 기능이 있다. 이 기능을 사용하여, Prestige 데이터셋의 income 변수와 정규성을 향해서 변환된 income 변수의 새로운 변수, 편의상 variable_Income,의 정규성 검정을 시도해보자.
먼저, 변환된 변수를 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)
R Commander를 설치하는 과정에서 의존패키지인 sem이 함께 설치된다. 안내가 나오면서 추가 설치를 하겠는가 물어보기도 한다. 위의 화면처럼, sem 패키지에 포함된 함수를 사용하는 <확인적 요인 분석...> 기능은 처음에는 비활성화되어 있다.
만약 Rcmdr 패키지가 호출될 때, sem 패키지가 자동으로 호출된다면, '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...'에 sem 패키지가 carData, sandwich 처럼 메뉴창안의 패키지 목록에 포함되어 있어야 할 것이다. 하지만, 위의 화면에선 보이지 않는다. sem 패키지를 추가로 호출해주어야 한다는 뜻이다.
library(sem) #설치한 sem 패키지를 호출하기
sem 패키지가 호출되면, '데이터 > 패키지에 있는 데이터 > 첨부된 패키지에서 데이터셋 읽기...'기능에서 sem 패키지의 데이터셋을 선택할 수 있게 된다. 아래 화면을 살펴보라.
HS.data를 선택하자. HS.data 데이터셋이 활성화되면, 처음에 미활성화된 <확인적 요인 분석...> 메뉴가 활성화된다. <확인적 요인 분석...> 메뉴를 선택하면, 두개의 하위 창을 보게된다. 변수들을 선택하여 요인으로 묶는 <데이터> 창과, 연산을 통하여 획득하고자하는 통계지수(index) 목록의 <선택기능> 창이다.
위에 보이는 스크립트는 무엇을 나타내는가? 설명을 추가한다. 1. 선택된 HS.data는 여러개의 변수들을 포함하고 있다. 2. 최초의 연구목적에서 개념적으로 구성된 몇 몇 요인들이 있고, 이 요인들을 세부적으로 구성하는 것이 하위 변수들이다. 3. 변수들 몇 개씩을 묶어서 요인으로 재구성하는것, 연구적 의도에서 보면, 요인을 구성한다고 가정하여 세분화된 변수들의 사례적 값들이 실제로 요인을 구성하는지를 확인하는 작업이 <확인적 요인 분석>이다. 4. HS.data에 포함된 여러개의 변수들을 spatial, verbal, memory, math라는 네개의 개념화된 요인으로 변수들을 재그룹화 한것이다. 'spatial: cubes, flags, paper, visual' 'verbal: general, paragrap, sentence, wordc, wordm' 'memory: figurer, figurew, numberf, numberr, object, wordr' 'math: arithmet, deduct, numeric, problemr, series' 5. spatial, verbal, memory, math 라는 요인의 이름은 최초의 설계에서 등장하는 개념적 요인을 뜻하는 것이다. 편의상으로 factor.1, factor.2, factor.3, factor.4 등으로 이름을 붙여도 무방하다. 6. fit.indices라는 옵션에 AIC, BIC 두개의 통계지수가 포함되어 있는데, 이것은 <선택기능> 창의 기본설정이며, 원하는 지수를 추가로 선택할 수 있다. CFI, RMSEA 등이 선택될 수 있다.