set.seed(19104)
y0 <- rnorm(50,23, 13)
y1 <-  rnorm(50, 28, 6)
y <- c(y0,y1)

x <- c(rep(0,50), rep(1,50))

dat <- cbind.data.frame(x,y)
head(dat)

table(dat$x)

mean(dat$y[x==1])
sd(dat$y[x==1])

mean(dat$y[x==0])
sd(dat$y[x==0])

#Test statistic
diff.means <- mean(dat$y[dat$x==1]) - mean(dat$y[dat$x==0])
diff.means


diff.means <- rep(NA, 10000)

for(i in 1:10000){
y0 <- rnorm(50,23, 13)
y1 <-  rnorm(50, 28, 6)
diff.means[i] <-  mean(y1) - mean(y0)
}

mean(diff.means)
sd(diff.means)

plot(density(diff.means),  main="Sampling Distribution of the Difference in Means")


sd(diff.means)

plot(density(diff.means), main="Sampling Distribution of the Difference in Means")

eval <- seq(0,20,.001)
points(eval, dnorm(eval, mean=5, sd=sd(diff.means)), type="l", col="firebrick")

se <- sqrt( (13^2/50)  + (6^2/50))
se

#Confirmed to be 2ish

(5.48-0)/2.02

pnorm(2.71, lower.tail=F)*2

se <- sqrt((var(dat$y[dat$x==0])/50) + (var(dat$y[dat$x==1])/50) )

num <- ((var(dat$y[dat$x==0])/50) + (var(dat$y[dat$x==1])/50))^2
denom1 <- ((var(dat$y[dat$x==0])/50)^2)/49
denom2 <- ((var(dat$y[dat$x==1])/50)^2)/49
df <- num/(denom1 +denom2)
df

5.48/se

pt(3.1, df=df, lower.tail=F)*2

t.test(dat$y, mu=25)

head(dat)
x <- t.test(dat$y ~ dat$x)
x

names(x)
x$stderr

eval <- seq(-5,15,.001)
plot(eval, dnorm(eval, mean=5, sd=2.02), type="l", col="firebrick", lwd=2, main="Power Analysis")
legend("topleft", c("True Sampling Distribution","Null Sampling Distribution"), lty=c(1,1), col=c("firebrick","darkblue"))

eval <- seq(-5,15,.001)
plot(eval, dnorm(eval, mean=5, sd=2.02), type="l", col="firebrick", lwd=2, main="Power Analysis")
points(eval, dnorm(eval, mean=0, sd=2.02), type="l", lwd=2, col="darkblue")
abline(v=1.96*2.02, lty=2, col="darkblue")
abline(v=-1.96*2.02, lty=2, col="darkblue")
legend("topleft", c("True Sampling Distribution","Null Sampling Distribution"), lty=c(1,1), col=c("firebrick","darkblue"))

pnorm(1.96*2.02, mean=5, sd=2.02) - pnorm(-1.96*2.02, mean=5, sd=2.02)


false.negative <- rep(NA, 10000)

for(i in 1:10000){
y0 <- rnorm(50,23, 13)
y1 <-  rnorm(50, 28, 6)
y <- c(y0,y1)
x <- c(rep(0,50), rep(1,50))
dat.sim <- cbind.data.frame(x,y)
test <- t.test(dat.sim$y ~ dat.sim$x)
false.negative[i] <- test$p.value>.05
}

mean(false.negative)


n <- 50
se <- sqrt( ((10^2)/n) +  ((10^2)/n) )

pnorm(1.96*se, mean=3, sd=se) - pnorm(-1*1.96*se, mean=3, sd=se)

true.diff <- seq(0,15,.01)
false.negative <- pnorm(1.96*se, mean=true.diff, sd=se) - pnorm(-1*1.96*se, mean=true.diff, sd=se)

plot(true.diff, false.negative, xlab="True difference between groups", ylab="False negative rate", type="l")


ctrl <- t.test(dat$y[dat$x==0])
treat <- t.test(dat$y[dat$x==1])

plot(ctrl$conf.int, c(0,0), ylim=c(-0.5,1.5), xlim=c(20,35),
     type="n", axes=F, xlab="Y",ylab="")
segments(ctrl$conf.int[1], 0,ctrl$conf.int[2],0)
segments(treat$conf.int[1], 1, treat$conf.int[2],1)
axis(side=2, at=c(0,1), labels=c("Control","Treatment"))
axis(side=1)

set.seed(19146)
y0 <- rnorm(50,23, 13)
y1 <-  rnorm(50, 28, 6)
y <- c(y0,y1)
x <- c(rep(0,50), rep(1,50))
dat.sim <- cbind.data.frame(x,y)


ctrl <- t.test(dat.sim$y[dat.sim$x==0])
treat <- t.test(dat.sim$y[dat.sim$x==1])

plot(ctrl$conf.int, c(0,0), ylim=c(-0.5,1.5), xlim=c(20,35),
     type="n", axes=F, xlab="Y",ylab="")
segments(ctrl$conf.int[1], 0,ctrl$conf.int[2],0)
segments(treat$conf.int[1], 1, treat$conf.int[2],1)
axis(side=2, at=c(0,1), labels=c("Control","Treatment"))
axis(side=1)

t.test(dat.sim$y ~ dat.sim$x)


reject.null.cis <- rep(NA)
reject.null.p <- rep(NA)

for(i in 1:10000){
y0 <- rnorm(50,23, 13)
y1 <-  rnorm(50, 28, 6)
y <- c(y0,y1)
x <- c(rep(0,50), rep(1,50))
dat.sim <- cbind.data.frame(x,y)

ctrl <- t.test(dat.sim$y[dat.sim$x==0])
treat <- t.test(dat.sim$y[dat.sim$x==1])

#Reject Null with CIs
reject.null.cis[i] <- treat$conf.int[1] > ctrl$conf.int[2]

#Reject Null with real test
p <- t.test(dat.sim$y ~ dat.sim$x)$p.value
reject.null.p[i] <- p<.05

}

table(reject.null.cis, reject.null.p)



set.seed(19104)
y0 <- rnorm(50,23, 13)
y1 <-  rnorm(50, 28, 6)
y <- c(y0,y1)
x <- c(rep(0,50), rep(1,50))
dat <- cbind.data.frame(x,y)
head(dat)

plot(dat$x, dat$y)
points(c(0,1), c(mean(dat$y[dat$x==0]), mean(dat$y[dat$x==1])), col="firebrick", pch=16)

set.seed(19103)
x <- rnorm(10,mean=0, sd=5)
y <- x*6*rnorm(10)
a <- cbind.data.frame(x,y)
plot(a$x, a$y, pch=16)

plot(a$x[a$x<0], a$y[a$x<0], pch=16, xlim=c(-12,7), col="darkblue")
points(a$x[a$x>0], a$y[a$x>0], pch=16, col="darkgreen")
segments(-20,mean(a$x[a$x<0]),0, mean(a$x[a$x<0]), col="darkblue")
segments(0,mean(a$x[a$x>0]),20, mean(a$x[a$x>0]), col="darkgreen")
mean(a$x[a$x>0]) -  mean(a$x[a$x<0])

set.seed(19103)
x <- rnorm(10,mean=0, sd=5)
y <- x*3 + rnorm(10, mean=0, sd=1)
b <- cbind.data.frame(x,y)

plot(b$x, b$y, pch=16, xlim=c(-12,7), col="darkblue")
points(b$x[b$x<0], b$y[b$x<0], pch=16, col="darkblue")
points(b$x[b$x>0], b$y[b$x>0], pch=16, col="darkgreen")
segments(-20,mean(b$x[b$x<0]),0, mean(b$x[b$x<0]), col="darkblue")
segments(0,mean(b$x[b$x>0]),20, mean(b$x[b$x>0]), col="darkgreen")
mean(b$x[b$x>0]) -  mean(b$x[b$x<0])

a
#Determine deviations from the column averages
a$x.dev <- a$x - mean(a$x)
a$y.dev <- a$y - mean(a$y)
a$product <- a$x.dev*a$y.dev
a
sum(a$product)/9

plot(a$x, a$y, pch=16)
abline(v=mean(a$x), lty=2)
abline(h=mean(a$y), lty=2)


set.seed(19103)
x <- rnorm(10,mean=0, sd=5)
y <- x*-6*rnorm(10)
c <- cbind.data.frame(x,y)
plot(a$x, a$y, pch=16)

plot(c$x, c$y, pch=16)
abline(v=mean(c$x), lty=2)
abline(h=mean(c$y), lty=2)

c$x.dev <- c$x - mean(c$x)
c$y.dev <- c$y - mean(c$y)
c$product <- c$x.dev*c$y.dev
c
sum(c$product)/9


plot(b$x, b$y, pch=16)
abline(v=mean(b$x), lty=2)
abline(h=mean(b$y), lty=2)

b$x.dev <- b$x - mean(b$x)
b$y.dev <- b$y - mean(b$y)
b$product <- b$x.dev*b$y.dev
b
sum(b$product)/9

set.seed(19106)
x <- rnorm(10,mean=0, sd=5)
y <- x*rnorm(10, sd=2)
d <- cbind.data.frame(x,y)


plot(d$x, d$y, pch=16)
abline(v=mean(d$x), lty=2)
abline(h=mean(d$y), lty=2)

d$x.dev <- d$x - mean(d$x)
d$y.dev <- d$y - mean(d$y)
d$product <- d$x.dev*d$y.dev
d
sum(d$product)/9

cov(a$x, a$y)
#The covariance function does the same thing we did above for us...

cov(a$x*100, a$y*100)

plot(a$x*100, a$y*100, pch=16)

x <- seq(1,10,1)
y <- seq(1,10,1)
plot(x,y, pch=16)

cov(x,y)
var(x)
var(y)
var(x)*var(y)
sqrt(var(x)*var(y))
cov(x,y)/sqrt(var(x)*var(y))

x[3] <- x[3]+.3
plot(x,y)
cov(x,y)
var(x)
var(y)
var(x)*var(y)
sqrt(var(x)*var(y))
cov(x,y)/sqrt(var(x)*var(y))



x <- seq(1,10,1)
y <- seq(10,1,-1)
plot(x,y, pch=16)

cov(x,y)
var(x)
var(y)
var(x)*var(y)
sqrt(var(x)*var(y))
cov(x,y)/sqrt(var(x)*var(y))

library(MASS)
cors <- seq(-.9,.7,.2)
cors

#par(mfrow=c(3,3))
for(i in 1:length(cors)){
sigma<-rbind(c(1,cors[i]), c(cors[i],1))
mu <- c(0,0)
d <- mvrnorm(n=100, mu=mu, Sigma=sigma)
plot(d[,1], d[,2], main=paste("Cor =",cors[i] ))
}