code 8.1

id<-c(1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8)
policy<-c("no","no","no","no","no","no","no","no","yes","yes","yes","yes","yes","yes","yes","yes")
outcome<-c(76,80,77,80,77,75,77,80,73,82,80,85,85,78,87,82)

p<-data.frame(id, policy, outcome)

tapply(outcome, policy, summary)

## $no
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   75.00   76.75   77.00   77.75   80.00   80.00 
## 
## $yes
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    73.0    79.5    82.0    81.5    85.0    87.0

var.test(p$outcome[p$policy=="no"], p$outcome[p$policy=="yes"])

## 
##  F test to compare two variances
## 
## data:  p$outcome[p$policy == "no"] and p$outcome[p$policy == "yes"]
## F = 0.19366, num df = 7, denom df = 7, p-value = 0.04583
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.03877187 0.96732397
## sample estimates:
## ratio of variances 
##           0.193662

t.test(p$outcome[p$policy=="no"], p$outcome[p$policy=="yes"], 
    conf.level=0.95, alternative="less", var.equal=FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  p$outcome[p$policy == "no"] and p$outcome[p$policy == "yes"]
## t = -2.1555, df = 9.6133, p-value = 0.02881
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##        -Inf -0.5838202
## sample estimates:
## mean of x mean of y 
##     77.75     81.50

code 8.2

# code 8.1의 데이터 활용

t.test(p$outcome[p$policy=="no"], p$outcome[p$policy=="yes"], paired=TRUE, 
        conf.level=0.95, alternative="less", var.equal=FALSE)

## 
##  Paired t-test
## 
## data:  p$outcome[p$policy == "no"] and p$outcome[p$policy == "yes"]
## t = -2.6576, df = 7, p-value = 0.01629
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf -1.076649
## sample estimates:
## mean of the differences 
##                   -3.75

code 8.3

df1<-7
df2<-7
s1<-1.9821
s2<-4.5040
alpha<-0.05

f.test=s1^2/s2^2

cu=qf(0.975, df1, df2)
cl=qf(0.025, df1, df2)

cu

## [1] 4.994909

cl

## [1] 0.2002038

f<-seq(0, 7, by=0.01)

f.dist<-data.frame(f)
f.dist$pdf<-with(f.dist, df(f, df1, df2))
f.dist$p.value.l<-with(f.dist, ifelse(f<cl, pdf, NA))
f.dist$p.value.u<-with(f.dist, ifelse(f>cu, pdf, NA))

p.area.l<-data.frame(f.dist$f, f.dist$p.value.l)
names(p.area.l)<-c("f", "pdf")

p.area.u<-data.frame(f.dist$f, f.dist$p.value.u)
names(p.area.u)<-c("f", "pdf")


library(ggplot2)

ggplot(f.dist, aes(x=f, y=pdf)) + geom_line() +
    geom_area(data=p.area.l, fill="blue", alpha=0.2) + 
    geom_area(data=p.area.u, fill="blue", alpha=0.2) +
    annotate("text", x=cl, y=0, label=round(cl,2), vjust=1, size=4) +
    annotate("text", x=cu, y=0, label=round(cu,2), vjust=1, size=4) 

## Warning: Removed 680 rows containing missing values (position_stack).

## Warning: Removed 500 rows containing missing values (position_stack).