과제2
김윤아
2024-06-03
1. 다음 데이터는 4가지 두통약인 아스피린(Aspirin), 타이레놀(Tylenol), 탁센(Naproxen), 게보린(Acetaminophen) 의 진통효과를 확인하기 위한 비교실험을 진행한 결과이다.
데이터는 진통효과의 크기를 나타내며 값이 클수록 진 통효과가 크다는 것을 의미한다.
실험은 4개의 병원(A, B, C, D)에서 나누어 진행되었으며 다음과 같은 데이터를 얻었다.
관심요인인 두통약 종류는 고정효과이며, 병원은 블록효과로 고려하였다.
(1) 어떠한 실험설계 방법을 사용하여 데이터를 수집하였는지 적으시오.
#또한, 처리의 수, 블록의 수, 블록 내의 처리 쌍의 수들을 계산하고, 그에 맞는 분할표 형태로 데이터를 표현하시오. #이 데이터는 블록설계를 사용하여 수집된것을 보임 #각 블록 내 처리수는 4, 블록의 수는 4, 블록내 처리 쌍의 수는 6이다
# 데이터 생성
hospital<-as.factor(c("A", "B", "C", "D"))
painkiller<-as.factor(c("Aspirin","Tylenol","Naproxen","Acetaminophen"))
y=c(17,18,19,NA,21,22,NA, 22, 15, NA, 17, 15, NA, 13, 13, 14)
# 데이터 프레임 생성
data <- data.frame(
painkiller = rep(painkiller, each = length(hospital)),
hospital = rep(hospital, times = length(painkiller)),
response = y
)
# 데이터 프레임 출력
print(data)## painkiller hospital response
## 1 Aspirin A 17
## 2 Aspirin B 18
## 3 Aspirin C 19
## 4 Aspirin D NA
## 5 Tylenol A 21
## 6 Tylenol B 22
## 7 Tylenol C NA
## 8 Tylenol D 22
## 9 Naproxen A 15
## 10 Naproxen B NA
## 11 Naproxen C 17
## 12 Naproxen D 15
## 13 Acetaminophen A NA
## 14 Acetaminophen B 13
## 15 Acetaminophen C 13
## 16 Acetaminophen D 14(2) 실험설계에 맞는 통계적 모형을 적으시오.
(3) 두통약에 따른 상자그림을 그리고, 두통약에 따른 평균 진통효과를 계산하여라.
# 데이터 시각화
plot(response ~ painkiller, data = data, main = "Response by Painkiller", xlab = "Painkiller", ylab = "Response")
plot(response ~ hospital, data = data, main = "Response by Hospital", xlab = "Hospital", ylab = "Response")
(4) 분산분석표를 작성하고, 두통약에 따라 진통효과에 차이가 있는지 확인하기 위한 가설을 설정하고, 유의수준 0.05 에서 가설검정하여라
data <- na.omit(data)
# 분산분석 수행
anovafit <- aov(response ~ painkiller, data = data)
# 분산분석 결과 출력
summary(anovafit)## Df Sum Sq Mean Sq F value Pr(>F)
## painkiller 3 113.7 37.89 50.52 1.52e-05 ***
## Residuals 8 6.0 0.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- 한 원예가는 A라는 기존 땅콩 종자와 2개의 다른 땅콩 종자(B, C)의 수확량 비교를 위해 실험을했다. 땅의 기울기와 토양의 질소 함유량을 고려한 라틴방격법을 사용한 실험계획을 하여 다음의 측정치를 얻었다. [40점]
통계적 모형을 적으시오. Y ijk=μ+αi+βj+γk+ϵ ijk
땅콩 종자별 상자그림을 그리시오.
# 데이터 입력
data <- data.frame(
Nitro = factor(rep(1:3, each=3)),
Slope = factor(c(1, 2, 3, 3, 1, 2, 2, 3, 1)),
Seed = c("C", "A", "B", "A", "B", "C", "B", "C", "A"),
Yield = c(26, 19, 28, 23, 20, 24, 28, 17, 28)
)
# 데이터 구조 확인
str(data)## 'data.frame': 9 obs. of 4 variables:
## $ Nitro: Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 3 3 3
## $ Slope: Factor w/ 3 levels "1","2","3": 1 2 3 3 1 2 2 3 1
## $ Seed : chr "C" "A" "B" "A" ...
## $ Yield: num 26 19 28 23 20 24 28 17 28# 라이브러리 로딩
library(ggplot2)
# 상자그림 그리기
ggplot(data, aes(x=Seed, y=Yield, fill=Seed)) +
geom_boxplot() +
labs(title="Peanut", x="Seed type", y="Yield") +
theme_minimal()
- 분산분석을 위한 R코드를 작성하고 분산분석표를 만드시오.
# 분산분석 실행
model <- lm(Yield ~ Seed + Slope + Nitro, data=data)
anova_result <- anova(model)
print(anova_result)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Seed 2 14 7 0.1228 0.8906
## Slope 2 6 3 0.0526 0.9500
## Nitro 2 8 4 0.0702 0.9344
## Residuals 2 114 57- 수확량이 땅콩 종자별, 기울기, 질소함유량에 따라 다른지 유의수준 0.05에서 검정하시오. 가설 설정 종자 유형(Seed)에 따른 수확량 차이
귀무 가설(H0): 모든 종자 유형의 평균 수확량이 동일하다. 대립 가설(H1): 적어도 한 종자 유형의 평균 수확량이 다른 종자 유형과 다르다. 종자의 유형의 F-값은 0.1228이고 p-값은 0.8906이므로 유의수준 0.05보다 크기 때문에 귀무가설을기각할수 없다 따라서 종자의 유형에 따른 수확량의 차이는 없다고 할 수 있습니다.
경사도(Slope)에 따른 수확량 차이
귀무 가설(H0): 모든 경사도의 평균 수확량이 동일하다. 대립 가설(H1): 적어도 한 경사도의 평균 수확량이 다른 경사도와 다르다.
경사도의 F 값은 0.0526이고 p-값은 0.9500이므로 귀무가설을 기각할수 없습니다. 따라서 경사도에 따른 평균 수확량은 차이가 없다고 할 수 있습니다.
질소 함유량(Nitro)에 따른 수확량 차이
귀무 가설(H0): 모든 질소 함유량 수준의 평균 수확량이 동일하다. 대립 가설(H1): 적어도 한 질소 함유량 수준의 평균 수확량이 다른 수준과 다르다.
질소 함유량의 F 값은 0.0702이고 p-값이 0.9344이므로 질소 함유량이 수확량에 미치는 영향도 통계적으로 유의하지 않습니다. 따러서 질소 함유량에 따른 수확량의 차이는 없다고 할수 있습니다.
(5) 농장이라는 효과(α, β, γ)를 추가로 고려하기 위해 그라코-라틴정방설계(Graco-Latin square design)을 고려하고자 한다. 3 × 3 그라코-라틴정방의 형태를 표현하고 통계적 모형을 작성하시오.
data$Farm = factor(rep(c("F1", "F2", "F3"), each=3))
# 그라코-라틴정방설계 모형
model_graco <- lm(Yield ~ Seed + Slope + Nitro + Farm, data=data)
anova_graco <- anova(model_graco)
print(anova_graco)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Seed 2 14 7 0.1228 0.8906
## Slope 2 6 3 0.0526 0.9500
## Nitro 2 8 4 0.0702 0.9344
## Residuals 2 114 57- 바이러스성장에 관한 연구에서, 바이러스가 성장하는데 바이러스의 종류와 배양시간, 배양액종류가 영향을 미친 다고 한다. 완전확률화설계 하에 23 요인 설계로 실험을 하여 바이러스 개체수를 반응값으로 해서 다음의 데이터를 얻었다.
- 위 실험에 대한 통계적 모형을 설계하시오. Y ijk=μ+αi+βj+γk+(αβ) ij+(αγ) ik+(βγ)jk+(αβγ) ijk+ϵ ijkl
𝑌𝑖𝑗𝑘:주어진 실험 조건에서의 바이러스 개체수입니다. 𝜇: 전체 평균입니다. 𝛼𝑖: i번째 바이러스 종류의 효과입니다. 𝛽𝑗: j번째 배양시간의 효과입니다. 𝛾𝑘: k번째 배양액 종류의 효과입니다. (𝛼𝛽)𝑖𝑗, (𝛼𝛾)𝑖𝑘, (𝛽𝛾)𝑗𝑘: 두 요인 간의 상호작용 효과입니다. (αβγ) ijk: 세 요인의 상호작용 효과입니다. ϵ ijkl: 실험 오차(잡음)를 나타냅니다.
- 완전모형에 대한 R코드를 이용해 분산분석표를 작성하시오.
# 데이터 프레임 생성
data <- data.frame(
virus_type = factor(rep(c("C0", "C1"), each=12)),
culture_time = factor(rep(rep(c("A0", "A1"), each=6), times=2)),
culture_medium = factor(rep(c("B0", "B1"), each=6, times=2)),
virus_count = c(3, 8, 9, 10, 14, 16, 4, 7, 9, 15, 19, 20, 10, 12, 16, 20, 21, 22, 10, 13, 16, 26, 27, 25)
)
# 데이터 구조 확인
str(data)## 'data.frame': 24 obs. of 4 variables:
## $ virus_type : Factor w/ 2 levels "C0","C1": 1 1 1 1 1 1 1 1 1 1 ...
## $ culture_time : Factor w/ 2 levels "A0","A1": 1 1 1 1 1 1 2 2 2 2 ...
## $ culture_medium: Factor w/ 2 levels "B0","B1": 1 1 1 1 1 1 2 2 2 2 ...
## $ virus_count : num 3 8 9 10 14 16 4 7 9 15 ...# 완전모형 분석을 위한 선형모델 적합
model <- lm(virus_count ~ virus_type * culture_time * culture_medium, data=data)
# 분산분석표 출력
anova_result <- anova(model)
print(anova_result)## Analysis of Variance Table
##
## Response: virus_count
## Df Sum Sq Mean Sq F value Pr(>F)
## virus_type 1 294.00 294.000 8.1253 0.009886 **
## culture_time 1 37.50 37.500 1.0364 0.320821
## virus_type:culture_time 1 0.17 0.167 0.0046 0.946564
## Residuals 20 723.67 36.183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- 배양액 종류에 대한 주효과의 크기 구하시오.
# 배양액 종류의 주효과 크기 구하기
summary(model)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.0000000 2.455719 4.07212723 0.0005943299
## virus_typeC1 6.8333333 3.472911 1.96760966 0.0631367529
## culture_timeA1 2.3333333 3.472911 0.67186671 0.5093559204
## virus_typeC1:culture_timeA1 0.3333333 4.911438 0.06786879 0.9465640192(4) 배양액과 바이러스 종류에 따른 상호작용효과에 대한 상호작용그림을 그리고 상호작용효과가 있는지 유의수준 0.05에서 가설검정하시오.
# 필요한 라이브러리 로딩
library(ggplot2)
# 상호작용 그림 그리기
ggplot(data, aes(x=virus_type, y=virus_count, color=culture_medium)) +
geom_point(position=position_dodge(width=0.1)) +
geom_line(aes(group=culture_medium), position=position_dodge(width=0.1)) +
labs(title="바이러스 종류와 배양액 종류의 상호작용 효과",
x="바이러스 종류",
y="바이러스 개체수") +
theme_minimal()
# 분산분석 수행
anova_result <- anova(model)
print(anova_result)## Analysis of Variance Table
##
## Response: virus_count
## Df Sum Sq Mean Sq F value Pr(>F)
## virus_type 1 294.00 294.000 8.1253 0.009886 **
## culture_time 1 37.50 37.500 1.0364 0.320821
## virus_type:culture_time 1 0.17 0.167 0.0046 0.946564
## Residuals 20 723.67 36.183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1귀무 가설 (H0): 모든 바이러스 종류에 대해 평균 바이러스 개체수는 같다. 대립 가설 (H1): 적어도 한 바이러스 종류의 평균 바이러스 개체수는 다르다. 바이러스 종류에 따른 F값이 8.1253이고 p-값이 0.009886이므로 귀무가설을 기각할 수 있습니다. 따라서 유의수준 0.05에서 바이러스 종류에 따른 바이러스 개체의 수에 차이가 있다고 할수 있습니다.
귀무 가설 (H0): 배양 시간에 따른 평균 바이러스 개체수의 차이는 없다. 대립 가설 (H1): 배양 시간에 따라 평균 바이러스 개체수가 다르다.
배양시간에 따른 F값은 1.0364이고 p-값은 0.320821이므로 귀무가설을 기각할 수 없습니다. 따라서 유의수준 0.05에서 배양 시간에 따른 바이러스 개체의 수에 차이가 없다고 할수 있습니다.
귀무 가설 (H0): 바이러스 종류와 배양 시간의 상호작용은 바이러스 개체수에 영향을 미치지 않는다. 대립 가설 (H1): 바이러스 종류와 배양 시간의 상호작용은 바이러스 개체수에 영향을 미친다.
바이러스 종류와 배양시간의 상호작용의 F-값은 0.0046이고 p-값은 0.946564이므로 귀무가설을 기각할 수 없습니다. 따라서 유의수준 0.05에서 두요인의 상호작용이 바이러스 개체수에 영향을 미치지 않는다고 할수 있습니다.
- A,B두요인에대해다음과같이교락법을실행하였다.각반복별로교락된효과를찾으시오.
data <- data.frame(
rep = rep(1:3, times = c(3, 2, 1)),
block = c(1, 2, 3, 4, 5, 6),
result = c("ab", "b", "ab", "ab", "a", "ab")
)
# 라이브러리 로딩
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# 각 반복에서 결과 분석
analysis <- data %>%
group_by(rep) %>%
summarize(
unique_results = length(unique(result)),
results = toString(unique(result))
)
# 출력
print(analysis)## # A tibble: 3 × 3
## rep unique_results results
## <int> <int> <chr>
## 1 1 2 ab, b
## 2 2 2 ab, a
## 3 3 1 ab# 각 반복별로 가능한 교락 효과 추정
# 예를 들어, 결과의 유니크한 패턴이나 빈도를 기반으로 분석
confounded_effects <- function(rep_data) {
if (length(unique(rep_data$result)) > 1) {
return("Possible interaction or main effect confounding")
} else {
return("No clear confounding")
}
}
# 반복별 교락 효과 추정
confounding_analysis <- data %>%
group_by(rep) %>%
summarize(
confounding = confounded_effects(cur_data())
)## Warning: There was 1 warning in `summarize()`.
## ℹ In argument: `confounding = confounded_effects(cur_data())`.
## ℹ In group 1: `rep = 1`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.# 출력
print(confounding_analysis)## # A tibble: 3 × 2
## rep confounding
## <int> <chr>
## 1 1 Possible interaction or main effect confounding
## 2 2 Possible interaction or main effect confounding
## 3 3 No clear confounding첫 번째 반복에서는 여러 결과값이 관찰되었으며, 이는 주요 요인의 효과나 상호작용 효과가 교락될 수 있음을 의미합니다. 첫번째반복에서는 결과값이 “ab”, “b”, “ab”로 나타난 것으로 보아 여러 요인의 효과가 혼재되어 있을 가능성이 있습니다. 두 번째 반복 역시 다양한 결과값(“ab”, “a”)이 관찰되어, 여러 요인 간의 상호작용이나 주요 효과의 교락이 있을 가능성을 나타냅니다. 이를 미루어 봤을 때 특정 요인의 효과를 명확히 구분하기 어려울 수 있습니다. 세 번째 반복에서는 결과값이 “ab” 하나로 일관되므로, 교락 효과가 명확하지 않습니다. 세번째 반복에서는 특정한 교락 효과가 없거나, 한 가지 결과만 나타났기 때문에 교락 효과를 추론하기 어렵다는 것을 의미합니다.