5장 랜덤화 블럭설계, 라틴정방설계와 분할법 - 예제

Size: px

Start display at page:

Download "5장 랜덤화 블럭설계, 라틴정방설계와 분할법 - 예제"

회원 옥
2 years ago
Views:

1 5 장랜덤화블럭설계, 라틴정방설계와분할법 - 예제 서울시립대통계학과

2 ii

3 차례 제1장 랜덤화 블럭설계 1 제1절 실험계획: 예제 5.1 -플라스틱 강도 제절 자료의 구성 제3절 시각적 분석 제4절 분산분석 제5절 혼합모형 제장 라틴정방설계 7 제1절 실험 계획: 예제 5. - 로켓 추진체 제절 자료의 구성 제3절 시각적 분석 제4절 분산분석 제5절 라틴정방의 구축 제3장 처리 조합의 블럭 13 제1절 실험 계획: 교재 분할법 I - 예제 화학약품의 생성률 제절 자료의 구성 제3절 시각적 분석 제4절 분산분석 제5절 블럭을 고려하지 않는 경우 제6절 혼합모형 제4장 Split-Plot desgin 3 제1절 실험 계획: 교재 분할법 II - 예제 전자제품 수명 제절 자료의 구성 제3절 시각적 분석 제4절 분산분석 iii

4 iv

5 서문 5장에 이원배치법에 대한 예제, R 프로그램과 결과 해석에 대한 강의자료입니다. 이 노트에 있는 R 프로그램을 실행하려면 다음과 같은 패키지들이 필요하다. library(dplyr) library(tidyr) library(ggplot) library(ggbeeswarm) library(agricolae) library(lme4) library(lmertest) v

6 vi

7 제1장 랜덤화 블럭설계 제1절 실험계획: 예제 5.1 -플라스틱 강도 플라스틱 제품의 강도를 측정하는 것이 실험의 목적이다. 랜덤하게 4일을 택해서 각 일마다 온도를 3개 수준으로 랜덤하게 변화시켜서 제품의 강도(intensity)를 측정하였다. 여기서 온도(temp)는 고정효과(τ )이며 선택된 일(day)는 블럭(ρ)에 따른 효과이다. xij = µ + τi + ρj + eij 제절 자료의 구성 이제 실험자료를 입력하여 데이터프레임으로 만들어 보자 intensity<- c(98.0, 97.7, 96.5, 99.0, 98.0, 97.9, 98.6, 98., 96.9, 97.6, 97.3, 96.7) temp<- factor(rep(c(70, 80, 90), times=4)) day<- as.factor(rep(c(1:4), each=3)) df<- data.frame(intensity=intensity, temp=temp, day=day) df 1 intensity temp day

8 벡터를 범주형 변수로 만들어 줄때 두 함수 as.factor() 와 factor() 모두 사용 가능하다. 제3절 시각적 분석 이제 온도의 수준에 따른 변화를 볼 수 있는 그림을 그려보자. 온도가 올라가면 강도가 떨어지는 경향을 볼 수 있다. df %>% ggplot(aes(x = temp, y = intensity, geom_line(aes(group = day)) + color=day)) + geom_point()

9 plot(intensity ~ temp, data=df) 3

10 이제 실험일에 따른 변동을 살펴보자. 실험일에 따라서 온도의 효과가 변하는 것을 볼 수 있다. 단 실험일과 온도의 상호작용은 크게 나타나지 않는다. 유의할 점은 반복이 없기 때문에 상호작용에 대한 추론은 불가능하다 df %>% ggplot(aes(x = day, y = intensity, geom_line(aes(group = temp)) + 제4절 color=temp)) + geom_point() 분산분석 블럭 효과인 실험일(day)를 고정효과로 놓았을 경우 분산분석표는 다음과 같다. ρj : fixed effect, eij N (0, σe ) model<- aov(intensity ~ temp + day, data=df) summary(model) temp Df Sum Sq Mean Sq F value Pr(>F) ** 4

11 day Residuals * -- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 위의 분산분석표에서 온도의 효과를 검정하는 F-통계량의 값은 이고 p-값은 0.007이다. 따라서 5% 유의수준으로 귀무가설을 기각하며 온도에 따라서 강도는 유의하게 다르다. 일반적으로 블럭효과에 대해서는 검정하지 않지만 그래도 p-값이 로서 매우 작으므로 실험일에 따른 변 동이 크다는 것을 알 수 있다. 이는 실험울 수행하는 날에 따라서 관측값에 변동이 크다는 것이다. 단 상호작용이 그림으로 볼 때 나타나지 않기 때문에 온도의 효과는 적절하게 추정할 수 있다. 제5절 혼합모형 고정효과와 임의효과(변량)가 동시에 모형식에 나타나는 모형을 혼합모형(mixed models)이라고 부른다. 혼합모형을 적합시키는 패키지는 lme4 이며 모형을 적합시키는 함수는 lmer이다. 혼합모형으로 부터 얻은 분산분석표에서 p-값을 보려면 패키지 lmertest를 사용해야 한다. 함수 lmer 에서 고정효과에 대한 모형식은 함수 anova와 같다. 함수 lmer 에서 만약 변수 var 을 임의효과로 고려하려면 (1 var) 으로 쓰면 된다. 다음은 플라스틱 강도 자료 실험에서 블럭 효과인 실험일(day, ρ)를 임의효과로 놓았을 경우 분석결과이다. 즉 ρj N (0, σb ), eij N (0, σe ) fit <- lmer(intensity ~ temp + (1 day), data=df) summary(fit) Linear mixed model fit by REML. t-tests use Satterthwaite's method [ lmermodlmertest] Formula: intensity ~ temp + (1 day) Data: df REML criterion at convergence: 14.6 Scaled residuals: Min 1Q Median Q Max Random effects: 5

12 Groups Name Variance Std.Dev. day (Intercept) Residual Number of obs: 1, groups: day, 4 Fixed effects: Estimate Std. Error df t value Pr(> t ) (Intercept) e-11 *** temp temp *** -- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) temp80 temp temp 의 추정치는 이며 오차항(Residual)의 분산 σe 위의 결과에서 블럭효과(day) 를 나타내는 분산 성분 σb 의 추정치는 이다. 이는 급내상관 계수(ICC)는 로서 매우 크다는 것을 의미한다. ICC = σb + σ = σb E 다음은 플라스틱 강도 자료 실험에서 블럭 효과를 임의효과로 놓았을 경우 분산분석표이다. 함수 lmer 에 의해 생 성된 결과를 함수 anova에 적용하면 고정효과에 대한 분산분석과 F-검정만 보여준다. 앞에서 블럭을 고정효과로 놓았을 때 분산분석의 검정 결과와 같다. anova(fit) Type III Analysis of Variance Table with Satterthwaite's method temp Sum Sq Mean Sq NumDF DenDF F value Pr(>F) ** -- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 6

13 제장 라틴정방설계 제1절 실험 계획: 예제 5. - 로켓 추진체 5가지의 로켓 추진체(A, B, C, D, E)의 성능을 비교하기 위하여 라틴정방계획을 사용한 실험이다. 행블럭: 5개의 연료 (R, ρ) 열블럭: 5명의 기사 (C, γ) 처리: 5가지의 로켓 추진체 (trt, τ ) xijk = µ + ρi + γj + τk + eijk 제절 자료의 구성 예제 5.에 있는 자료를 분석을 위하여 데이터프레임으로 만들어 보자. trt <- c("a", "B", "C", "D", "E", "B", "C", "D", "E", "A", "C", "D", "E", "A", "B", "D", "E", "A", "B", "C", "E", "A", "B", "C", "D" ) trt <- factor(trt) R <- factor(rep(1:5, each=5)) C <- factor(rep(1:5, times=5)) y <- c( -1,-5, -6, -1, -1, -8, -1, 5,, 11, -7, 13, 1,, -4, 7

14 1, 6, 1, -, -3, -3, 5, -5, 4, 6) df<- data.frame(trt, R, C, y) df trt R C y 1 A B C D E B C -1 8 D E 4 10 A C D E A B D E A B C E A B C D 함수 tabs() 는모형식을이용하여다음과같이열과행으로구성된자료를보여줄수있다. xtabs(y~ R + C, data = df) C R

15 제3절 시각적 분석 먼저 로켓 추진체, 즉 처리별로 자료의 분포를 보자. 추진체 B 와 C 가 다른 추진체들 보다 관측값이 작게 나오는 것을 알 수 있다. df %>% ggplot() + aes(x = trt, y = y) + geom_boxplot() 원료(R) 뭉치별로 자료의 분포를 보면 큰 차이는 보이지 않는다. df %>% ggplot() + 9

16 aes(x = R, y = y) + geom_boxplot() 기사 (C) 별로자료의분포를보면약간의차이가보인다. df %>% ggplot() + aes(x = C, y = y) + geom_boxplot() 10

17 제4절 분산분석 이제 라틴정방계획법으로 얻은 자료에 대래 분산분석을 적용해 보자. model<- aov(y ~ trt + R + C, data=df) summary(model) Df Sum Sq Mean Sq F value Pr(>F) trt ** R C * Residuals --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 위의 분산분석표에서 추진체(처리)의 효과를 검정하는 F-통계량의 값은 이고 p-값은 0.005이다. 따라서 5% 유의수준으로 귀무가설을 기각하며 추진체에 따라서 성능이 유의하게 다르다. 11

18 제5절 라틴정방의 구축 교과서 5.3절에서는 라틴정방 계획으로 실험을 하는 경우 처리를 랜덤하게 배정하는 방법을 설명하고 있다. 패키지 agricolae 에 포함된 함수 design.lsd()를 이용하면 다음과 같이 처리를 랜덤하게 배정해준다. mytrt <- factor(c("a", "B", "C", "D", "E")) mytrt [1] A B C D E Levels: A B C D E design.lsd(mytrt)$sketch [,1] [,] [,3] [,4] [,5] [1,] "C" "D" "B" "A" "E" [,] "E" "A" "D" "C" "B" [3,] "D" "E" "C" "B" "A" [4,] "B" "C" "A" "E" "D" [5,] "A" "B" "E" "D" "C" 함수 design.lsd()는 실행할 때마다 랜덤하게 배정하기 때문에 기록을 위해서 랜덤 seed 를 지정하면 나중에도 동일한 계획을 얻을 수 있다. design.lsd(mytrt, seed = 134 )$sketch [,1] [,] [,3] [,4] [,5] [1,] "C" "B" "E" "A" "D" [,] "A" "E" "C" "D" "B" [3,] "B" "A" "D" "E" "C" [4,] "D" "C" "A" "B" "E" [5,] "E" "D" "B" "C" "A" 1

19 제3장 처리 조합의 블럭 제1절 실험 계획: 교재 분할법 I - 예제 화학약품의 생성률 이 실험에서는 화학약품의 생성률에 영향을 미치는 두 요인을 고려한 실험이다. 반응온도(temp, α) 3개의 수준 중간원료 제조회사 (company, β) 3개의 수준 이 실험에서는 9개의 처리를 먼저 랜덤하게 선택하고 선택된 처리 하에서 실험을 번 반복하였다. 따라서 처리의 조합이 블럭효과(block, ρ)로 나타난다. xijk = µ + αi + βj + ρij + e(ijk) 위의 모형식에서 상호작용 효과 (αβ)ij 와 1차 랜덤화에 의한 오차 e1(ij) 는 교락되어 블럭효과 ρij 에 합쳐저서 나타난다. ρij = e1(ij) + (αβ)ij 이러한 경우 블럭효과 ρij 는 임의효과가 된다. ρij N (0, σ1 ), 제절 e(ijk) N (0, σ ) 자료의 구성 이제 실험자료를 입력하여 데이터프레임으로 만들어 보자 13 (3.1)

20 temp<- as.factor(rep(c("a1","a", "A3"), each=, times=3)) company<- as.factor(rep(c("b1", "B", "B3"), each=6)) y <-c( 81.0, 80., 84.1, 83., 85., 86.1, 83.3, 8.7, 86., 85.4, 86.6, 87., 81.3, 81.9, 83., 84., 86.0, 86.4) df<- data.frame(temp, company, y) df temp company y 1 A1 B A1 B A B A B A3 B A3 B A1 B A1 B A B A B A3 B A3 B A1 B A1 B A B A B A3 B A3 B 제 3 절 시각적분석 일단각처리에대한관측값의평균을구해보자. dfsum <- df %>% group_by(temp, company) %>% summarise(mean=mean(y), sd=sd(y)) `summarise()` has grouped output by 'temp'. You can override using the `.groups` argument. 14

21 dfsum # A tibble: 9 x 4 # Groups: temp [3] temp company <fct> <fct> mean sd <dbl> <dbl> 1 A1 B A1 B 83 3 A1 B A B A B A B A3 B A3 B A3 B 이제 처리의 평균값을 가지고 온도에 따른 변화를 살펴보자. 이 경우 제조회사 원료에 대해서는 색깔을 다르게 하여 상호작용 효과도 볼 수 있다. 아래 상호작용 그림을 보면 온도에 따라서 화학약품의 생성률이 크게 변하는 것을 알 수 있다. 유의한 상호작용은 관측되지 않는다. dfsum %>% ggplot(aes(x = temp, y = mean, geom_line(aes(group = company)) + color=company)) + geom_point() 15

22 함수 interaction.plot() 은상호작용그림을평균값을계산하지않고원래자료를이용하여다음과같이그릴 수있다. with(df, interaction.plot(x.factor = temp, trace.factor = company, response = y)) 16

23 위에서 함수 with() 은 이용하고자 하는 변수가 있는 데이터프레임을 지정하는데사용한다. 함수 with()의 첫 번쨰 인자는 앞의 예제와 같이 df 와 같은 데이터 프레임을 지정한다. 두 번째 인자에는 함수를 이용한 명 령문을 넣어준다. 앞의 프로그램에서 함수 interaction.plot() 안에서 사용된 변수들( temp,company,y) 들은 데이터프레임 df에 있는 변수들이다. 이제 제조회사에 따른 변화를 살펴보자. 제조회사에 따른 생성률의 변화는 크지 않다. with(df, interaction.plot(x.factor = company, trace.factor =temp, 17 response = y))

24 제4절 분산분석 이제 분산분석을 하여 처리의 효과에 대한 검정을 해보자. 실험에서 각 처리의 조합을 블럭으로 해주어야 한다. 다음 anova 함수에서 두 처리의 조합을 temp:company 로 표시한다. 사실 temp:company는 두 처리 temp와 company의 상호작용(interaction)을 의미한다. 다음으로 처리의 조합 temp:company 이 임의효과라는 것을 Error(temp:company)와 같이 지정해 준다. model<- aov(y ~ temp + company + Error(temp:company), data=df) Warning in aov(y ~ temp + company + Error(temp:company), data = df): Error() model is singular summary(model) Error: temp:company Df Sum Sq Mean Sq F value temp company Residuals Pr(>F) *** * -- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 18

25 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 위의 분산분석표에서 온도의 효과를 검정하는 F-통계량의 값은 이고 p-값은 이다. 따라서 5% 유의수준으로 귀무가설을 기각하며 온도에 따라서 생성률이 매우 유의하게 다르다. 온도의 효과를 검정하는 F-통계량의 값은 이고 p-값은 이다. 따라서 5% 유의수준으로 귀무가설을 기각하며 원료 제조회사에 따라서도 생성률이 유의하게 다르다. 제5절 블럭을 고려하지 않는 경우 만약에 처리 조합으로 생긴 블럭효과를 고려하지 않으면 어떤 일이 일어날까? 만약 생성률 실험자료를 완전 랜덤화 이원배치법에 의하여 얻은 자료라고 생각한다면 반복이 있으므로 상호작용 효과를 추론할 수 있다. 따라서 상호작용 효과를 고정효과로 놓고 분산분석을 적용할 것이다. ρij = (αβ)ij : fixed effect, e(ijk) N (0, σ ) (3.) 아래 프로그램은 상호작용 효과를 고정효과로 생각한 것이다. model<- aov(y ~ temp + company + temp:company, data=df) summary(model) Df Sum Sq Mean Sq F value Pr(>F) temp e-07 *** company *** temp:company Residuals Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 분산분석의 결과는 위와 같으며 온도와 제조회사에 대한 F-검정 통계량을 보면 임의효과 모형에서 나온 것보다 크다. 이는 F-검정 통계량을 만들 때 분모에 사용된 평균 오차제곱합 M SE 와 자유도가 달라서 나타나는 현상이다. 또한 자유도도 두 모형에서 온도에 대한 F-검정의 차이를 보자. 모형 anova 항 M SA M SE F0 임의효과 모형 [식 (3.1)] Error(temp:company) 고정효과 모형 [식 (3.)] temp:company

26 위의 표에서와 같이 실험계획에 따라서 나누어 주는 평균 오차제곱합 M SE 와 자유도가 다르기 때문에 검정의 결과가 다르게 나타난다. 실험계획에서 통계적 추론을 하는 경우 자료의 구조는 같아도 실험의 방법(랜덤화의 방법)이 다르면 가설 검정의 방법이 다르다. 따라서 실험의 방법에 따른 적절한 통계적 추론 방법을 선택하는 것이 중요하다. 제6절 혼합모형 처리들의 조합을 임의효과로 보는 모형 (3.1) 을 lmer로 적합시키는 프로그램은 다음과 같다. 분산분석 결과는 anova() 에서 임의효과 Error(temp:company)를 사용하는 결과와 동일하다. fit <- lmer(y ~ temp + company + (1 temp:company ), data = df) summary(fit) Linear mixed model fit by REML. t-tests use Satterthwaite's method [ lmermodlmertest] Formula: y ~ temp + company + (1 temp:company) Data: df REML criterion at convergence: 9.4 Scaled residuals: Min 1Q Median 3Q Max Random effects: Groups Name Variance Std.Dev. temp:company (Intercept) Residual Number of obs: 18, groups: temp:company, 9 Fixed effects: (Intercept) Estimate Std. Error df t value Pr(> t ) tempa ** tempa *** companyb ** 0 1.4e-09 ***

27 companyb Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) tempa tempa3 cmpnb tempa tempa companyb companyb anova(fit) Type III Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) temp *** company * --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 1

29 제4장 Split-Plot desgin 제1절 실험 계획: 교재 분할법 II - 예제 전자제품 수명 전자부품의 수명이 온도(580, 600, 60, 640도)와 시간(5, 10, 15분)에 의해 어떤 영향을 받는지에 대한 실험이다. 이 실험은 split-plot 설계를 적용하여 관측값을 얻었다. 온도를 먼저 랜덤하게 선택하고 선택된 온도에서 3개의 가열 시간에 대한 실험을 임의 순서로 진행하였다. 또한 각 실험은 3번 반복 하였다. 온도 (temp, α) : 주구, main plot - 1차 랜덤화 요인 시간 (time, β) : 분할구, split-plot, sub-plot - 차 랜덤화 요인 반복 (rep, r) : 반복 요인 xijk = µ + rk + αi + γik + βj + (αβ)ij + e(ijk) (4.1) 위의 모형식에서 반복과 온도의 상호작용 효과 (αr)ik 와 1차 랜덤화에 의한 오차 e1(ik) 는 교락되어 블럭효과 γik 에 합쳐저서 나타난다. γik = (αr)ik + e1(ik) 제절 자료의 구성 이제 실험자료를 입력하여 데이터프레임으로 만들어 보자 rep<- as.factor(rep(c(1:3), each=1)) temp<- as.factor(rep(c(580, 600, 60, 640), each=3, times=3)) time<- as.factor(rep(c(5, 10, 15), times=1)) 3

30 y <-c(17, 33, 175, 158, 138, 15, 9, 186, 155, 3, 7, 156, 188, 01, 195, 16, 130, 147, 160, 170, 161, 01, 181, 17, 16, 170, 13, 1, 185, 180, 167, 181, 18, 18, 01, 199) df <- data.frame(rep, temp, time, y) 함수 xtab 을 이용하면 반복에 따라서 자료 구조를 쉽게 볼 수 있다. xtabs( y ~time + temp + rep, df),, rep = 1 temp time ,, rep = temp time ,, rep = 3 temp time 제3절 시각적 분석 이제 온도의 수준에 따른 변화를 볼 수 있는 그림을 그려보자. 온도가 증가하면서 수명이 줄어들었다가 다시 늘어나는 현상을 볼 수 있다. 4

31 with(df, interaction.plot(x.factor = temp, trace.factor = time, response = y)) 가열시간의수준에따른변화를볼수있는그림을그려보자. 가열시간이증가하더러도수명이크게변하지않는 것을알수있다. with(df, interaction.plot(x.factor = time, trace.factor = temp, response = y)) 5

32 제4절 분산분석 이제 모형식 (4.1) 에 대한 분산분석을 실시해 보자. 여기서 유의할 점은 모형식 (4.1) 에서 블럭효과 γik 는 임의효과로 생각하며 반복 수준과 온도 수준의 조합이다. 따라서 블럭효과 γik 에 대한 항을 Error(rep:temp)로 사용한다. γik N (0, σ1 ), e(ijk) N (0, σe ) model<- aov(y ~ rep + temp*time + Error(rep:temp), data=df) Warning in aov(y ~ rep + temp * time + Error(rep:temp), data = df): Error() model is singular summary(model) Error: rep:temp Df Sum Sq Mean Sq F value Pr(>F) rep temp Residuals ** -- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: Within Df Sum Sq Mean Sq F value Pr(>F) time temp:time Residuals 분산분석표에서 온도의 효과를 검정하는 F-통계량의 값은 이고 p-값은 0.004이다. 따라서 5% 유의수준 으로 귀무가설을 기각하며 온도에 따라서 제품의 수명이 유의하게 다르다. 가열시간의 효과를 검정하는 F-통계량의 값은 이고 p-값은 이다. 따라서 5% 유의수준으로 귀무가 설을 기각할 수 없으며 가열시간에 따라서 제품의 수명이 다르지 않다. 온도와 가열시간의 상호작용 효과를 검정하는 F-통계량의 값은 이고 p-값은 이다. 따라서 5% 유의수준으로 귀무가설을 기각할 수 없으며 상호작용은 유의하지 않다. 6

Microsoft PowerPoint - IPYYUIHNPGFU

Microsoft PowerPoint - IPYYUIHNPGFU 분산분석 분산분석 (ANOVA: ANALYSIS OF VARIANCE) 두개이상의모집단의차이를검정 예 : 회사에서세종류의기계를설치하여동일한제품을생산하는경우, 각기계의생산량을조사하여평균생산량을비교 독립변수 : 다른변수에의해영향을주는변수 종속변수 : 다른변수에의해영향을받는변수 요인 (Factor): 독립변수 예에서의요인 : 기계의종류 (I, II, III) 요인수준