#Title: 0928 Homework1
#Name: Szu-Yu Chen
#Date: 4 Oct,2020


library(PairedData)

## Loading required package: MASS

## Loading required package: gld

## Loading required package: mvtnorm

## Loading required package: lattice

## Loading required package: ggplot2

## 
## Attaching package: 'PairedData'

## The following object is masked from 'package:base':
## 
##     summary

data(Anorexia, package="PairedData")
dta01 <- Anorexia
str(dta01)

## 'data.frame':    17 obs. of  2 variables:
##  $ Prior: num  83.8 83.3 86 82.5 86.7 79.6 76.9 94.2 73.4 80.5 ...
##  $ Post : num  95.2 94.3 91.5 91.9 100.3 ...

head(dta01)

##   Prior  Post
## 1  83.8  95.2
## 2  83.3  94.3
## 3  86.0  91.5
## 4  82.5  91.9
## 5  86.7 100.3
## 6  79.6  76.7

#Display data mean, SD, and Correlation
colMeans(dta01)

##    Prior     Post 
## 83.22941 90.49412

print(apply(dta01, 2, sd), 3)

## Prior  Post 
##  5.02  8.48

cov(dta01)

##          Prior     Post
## Prior 25.16721 22.88268
## Post  22.88268 71.82684

print(cor(dta01),3)

##       Prior  Post
## Prior 1.000 0.538
## Post  0.538 1.000

#Paired t-test output
t.test(dta01$Post, dta01$Prior, pair=T)

## 
##  Paired t-test
## 
## data:  dta01$Post and dta01$Prior
## t = 4.1849, df = 16, p-value = 0.0007003
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   3.58470 10.94471
## sample estimates:
## mean of the differences 
##                7.264706

# Reshape data_wide to long
dta011 <- tidyr::gather(dta01, key = "Treatment", value = "Response", Prior, Post)
head(dta011)

##   Treatment Response
## 1     Prior     83.8
## 2     Prior     83.3
## 3     Prior     86.0
## 4     Prior     82.5
## 5     Prior     86.7
## 6     Prior     79.6

#Plot the data

ggplot(dta011, aes(Treatment, Response)) + 
  geom_point(shape = 1,
             colour = "black") + 
  stat_summary(fun.data = mean_cl_boot, 
               geom = "pointrange",
               colour = "black") + 
               coord_flip() +
    labs(x = "Treatment", y = "Response") +
  theme_bw()

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

# output
summary(lm(Response ~ Treatment -1, data=dta011))

## 
## Call:
## lm(formula = Response ~ Treatment - 1, data = dta011)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.294  -2.454   1.106   4.004  11.106 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## TreatmentPost    90.494      1.689   53.58   <2e-16 ***
## TreatmentPrior   83.229      1.689   49.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.964 on 32 degrees of freedom
## Multiple R-squared:  0.994,  Adjusted R-squared:  0.9936 
## F-statistic:  2649 on 2 and 32 DF,  p-value: < 2.2e-16

# generate subject
dta011 <- dplyr::mutate(dta011, Subject=rep(paste0("S", 1:17), 2), Sbject=paste0("S", 1:34))
# within subject design output
summary(aov(Response ~ Treatment + Error(Subject/Treatment), data = dta011))

## 
## Error: Subject
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 16   1142   71.38               
## 
## Error: Subject:Treatment
##           Df Sum Sq Mean Sq F value Pr(>F)    
## Treatment  1  448.6   448.6   17.51  7e-04 ***
## Residuals 16  409.8    25.6                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#the end

```