code 11.1
high<-c(100,91,66,66,83,75,88,93,88,56,NA)
med<-c(84,54,89,41,40,65,77,69,73,88,75)
low<-c(47,73,67,35,68,53,78,57,52,NA,NA)
work<-data.frame(high, med, low)
temp<-c()
work2<-c()
#ANOVA를 할 수 있도록 데이터프레임 변환
for(i in 1:nrow(work))
{
for(j in 1:ncol(work))
{
compet<-work[i,j]
level<-names(work)[j]
temp<-data.frame(compet, level)
work2<-rbind(work2, temp)
}
}
#서열형 변수로 level 재구성
work2$level<-ordered(work2$level, levels=c("low", "med", "high"))
ordered(work2$level, levels=c("low", "med", "high"))## [1] high med low high med low high med low high med low high med
## [15] low high med low high med low high med low high med low high
## [29] med low high med low
## Levels: low < med < highhead(work2)## compet level
## 1 100 high
## 2 84 med
## 3 47 low
## 4 91 high
## 5 54 med
## 6 73 lowresult.11.1<-aov(compet ~ level, data=work2) #ANOVA결과 저장
summary(result.11.1) #ANOVA결과 요약 출력## Df Sum Sq Mean Sq F value Pr(>F)
## level 2 2252 1126.2 4.831 0.0161 *
## Residuals 27 6294 233.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 3 observations deleted due to missingnesscode 11.2
gender<-c("f","f","f","m","m","m")
age50<-c(74,76,82,77,83,92)
age40<-c(84,81,90,89,93,97)
age30<-c(85,57,87,76,84,89)
age20<-c(66,60,72,92,90,77)
welfare<-data.frame(gender, age50, age40, age30, age20)
temp<-c()
welfare2<-c()
#ANOVA가 가능하도록 데이터프레임 변환
for(i in 1:nrow(welfare))
{
for(j in 1:(ncol(welfare)-1))
{
wel<-welfare[i,j+1]
agelevel<-names(welfare)[j+1]
gen<-welfare[i,1]
temp<-data.frame(wel, agelevel, gen)
welfare2<-rbind(welfare2, temp)
}
}
head(welfare2)## wel agelevel gen
## 1 74 age50 f
## 2 84 age40 f
## 3 85 age30 f
## 4 66 age20 f
## 5 76 age50 f
## 6 81 age40 f#서열형 변수로 설정
welfare2$agelevel<-ordered(welfare2$agelevel,
levels=c("age20", "age30", "age40", "age50"))
#ANOVA 결과 저장
result.11.2<-aov(wel ~ agelevel + gen + agelevel*gen, data=welfare2)
#ANOVA 결과 출력
summary(result.11.2)## Df Sum Sq Mean Sq F value Pr(>F)
## agelevel 3 532.1 177.4 2.638 0.08511 .
## gen 1 651.0 651.0 9.681 0.00672 **
## agelevel:gen 3 198.5 66.2 0.984 0.42521
## Residuals 16 1076.0 67.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#ANOVA 결과 출력. 전체 F-test값 출력
result.11.2.1<-lm(wel ~ agelevel + gen + agelevel*gen, data=welfare2)
summary(result.11.2.1)##
## Call:
## lm(formula = wel ~ agelevel + gen + agelevel * gen, data = welfare2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.333 -4.000 0.000 5.167 10.667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.167 2.367 32.174 5.7e-16 ***
## agelevel.L 9.541 4.735 2.015 0.06102 .
## agelevel.Q -9.000 4.735 -1.901 0.07549 .
## agelevel.C -3.280 4.735 -0.693 0.49844
## genm 10.417 3.348 3.111 0.00672 **
## agelevel.L:genm -8.870 6.696 -1.325 0.20389
## agelevel.Q:genm 6.167 6.696 0.921 0.37074
## agelevel.C:genm -3.950 6.696 -0.590 0.56344
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.201 on 16 degrees of freedom
## Multiple R-squared: 0.5622, Adjusted R-squared: 0.3706
## F-statistic: 2.935 on 7 and 16 DF, p-value: 0.0353code 11.3
install.packages(“car”)
performance<-c("l","l","l","h","h","h")
commit.h<-c(1.4,2.4,2.2,2.4,NA,NA)
commit.m<-c(2.1,1.7,NA,2.5,1.8,2)
commit.l<-c(0.7,1.1,NA,0.5,0.9,1.3)
service<-data.frame(performance, commit.h, commit.m, commit.l)
service2<-c()
for(i in 1:nrow(service))
{
for(j in 1:(ncol(service)-1))
{
performance<-service[i,j+1]
commitment<-names(service)[j+1]
p.pay<-service[i,1]
temp<-data.frame(performance, commitment, p.pay)
service2<-rbind(service2, temp)
}
}
service2$commitment<-ordered(service2$commitment,
c("commit.l", "commit.m", "commit.h"))
service2$p.pay<-ordered(service2$p.pay, c("l", "h"))
head(service2)## performance commitment p.pay
## 1 1.4 commit.h l
## 2 2.1 commit.m l
## 3 0.7 commit.l l
## 4 2.4 commit.h l
## 5 1.7 commit.m l
## 6 1.1 commit.l lresult.11.3<-aov(performance ~ commitment + p.pay + commitment*p.pay, data=service2)
summary(result.11.3)## Df Sum Sq Mean Sq F value Pr(>F)
## commitment 2 4.306 2.1531 13.250 0.00289 **
## p.pay 1 0.093 0.0926 0.570 0.47202
## commitment:p.pay 2 0.075 0.0377 0.232 0.79803
## Residuals 8 1.300 0.1625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 4 observations deleted due to missingness## 1, 2, 3종 제곱합
library(car)
result.lm<-lm(performance ~ commitment + p.pay + commitment*p.pay, data=service2)
Anova(result.lm, type=2) #car 패키지의 Anova 함수를 사용, 2/3종 제곱합 산출## Anova Table (Type II tests)
##
## Response: performance
## Sum Sq Df F value Pr(>F)
## commitment 4.3960 2 13.5262 0.002713 **
## p.pay 0.0926 1 0.5697 0.472022
## commitment:p.pay 0.0754 2 0.2321 0.798034
## Residuals 1.3000 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(result.lm, type=3)## Anova Table (Type III tests)
##
## Response: performance
## Sum Sq Df F value Pr(>F)
## (Intercept) 34.680 1 213.4154 4.729e-07 ***
## commitment 4.190 2 12.8914 0.003145 **
## p.pay 0.120 1 0.7385 0.415160
## commitment:p.pay 0.075 2 0.2321 0.798034
## Residuals 1.300 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# interaction.plot사용하여 평균반응그림 그리기
service3<-subset(service2, !is.na(performance)) #NA가 있으면 계산이 안 됨
interaction.plot(service3$commitment, service3$p.pay, service3$performance)
# 무식하게 그려보자
library(ggplot2)
library(plyr)
head(service2)## performance commitment p.pay
## 1 1.4 commit.h l
## 2 2.1 commit.m l
## 3 0.7 commit.l l
## 4 2.4 commit.h l
## 5 1.7 commit.m l
## 6 1.1 commit.l lmean.plot<-ddply(service2, c("p.pay", "commitment"), summarize,
performance=mean(performance, na.rm=T))
mean.plot$commitment<-as.numeric(mean.plot$commitment)
ggplot(service2, aes(x=commitment, y=performance, color=p.pay)) + geom_point() +
geom_line(data=mean.plot)## Warning: Removed 4 rows containing missing values (geom_point).
code 11.4
# code 11.3의 자료를 사용
result.11.4<-aov(wel ~ agelevel, data=welfare2)
summary(result.11.4)## Df Sum Sq Mean Sq F value Pr(>F)
## agelevel 3 532.1 177.37 1.842 0.172
## Residuals 20 1925.5 96.27TukeyHSD(result.11.4, conf.level=0.90)## Tukey multiple comparisons of means
## 90% family-wise confidence level
##
## Fit: aov(formula = wel ~ agelevel, data = welfare2)
##
## $agelevel
## diff lwr upr p adj
## age30-age20 3.500000 -10.365991 17.365991 0.9251770
## age40-age20 12.833333 -1.032658 26.699324 0.1400860
## age50-age20 4.500000 -9.365991 18.365991 0.8561531
## age40-age30 9.333333 -4.532658 23.199324 0.3761354
## age50-age30 1.000000 -12.865991 14.865991 0.9979786
## age50-age40 -8.333333 -22.199324 5.532658 0.4724762plot(TukeyHSD(result.11.4, conf.level=0.90))
?TukeyHSD## starting httpd help server ... donecode 11.5
install.packages(“lsr”)
library(lsr)
etaSquared(result.11.2, type = 3, anova = TRUE )## eta.sq eta.sq.part SS df MS F
## agelevel 0.22311513 0.3375744 548.3333 3 182.77778 2.7178852
## gen 0.26490684 0.3769693 651.0417 1 651.04167 9.6809170
## agelevel:gen 0.08075208 0.1557198 198.4583 3 66.15278 0.9836844
## Residuals 0.43782107 NA 1076.0000 16 67.25000 NA
## p
## agelevel 0.079137044
## gen 0.006717264
## agelevel:gen 0.425208100
## Residuals NA