causal-inference-week5.Rmd

---
title: "casual inference week 5"
author: "Zihao Wang"
date: "2017�6�15�"
output: html_document
---

#task 1:
for K fixed, draw the plot between the statistics power and the sample size N.
```{r}
simulation<-function(N,K,n,Ybar,sigma,M){
  
Ybar.sim<-matrix(nrow = K,ncol = M)
variance.sim<-matrix(nrow = K,ncol = M)
for (i in 1:K){
for (j in 1:M){
  Ybar.sim[i,j]=rnorm(1,Ybar[i],sigma[i]*sqrt(K/N))
  variance.sim[i,j]=rchisq(1,n[i]-1)*(sigma[i]^2)/(n[i]-1)
}
}

Vneyman.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M) {
    Vneyman.sim[i,j]<-variance.sim[i+1,j]/n[i]+variance.sim[1,j]/n[1]
  }
}

theta.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M) {
    theta.sim[i,j]<-(Ybar.sim[i+1,j]-Ybar.sim[1,j])/sqrt(Vneyman.sim[i,j])
  }
}

pvalue.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M){
    pvalue.sim[i,j]=2*pnorm(-abs(theta.sim[i,j]))
  }
}

pcombine.sim<-as.numeric()
for(j in 1:M) {pcombine.sim[j]=min(pvalue.sim[,j])}

pcombine.sim

}
```
```{r}
simulation.mean<-function(N,K,n,Ybar,sigma,M){
  
Ybar.sim<-matrix(nrow = K,ncol = M)
variance.sim<-matrix(nrow = K,ncol = M)
for (i in 1:K){
for (j in 1:M){
  Ybar.sim[i,j]=rnorm(1,Ybar[i],sigma[i]*sqrt(K/N))
  variance.sim[i,j]=rchisq(1,n[i]-1)*(sigma[i]^2)/(n[i]-1)
}
}

Vneyman.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M) {
    Vneyman.sim[i,j]<-variance.sim[i+1,j]/n[i]+variance.sim[1,j]/n[1]
  }
}

theta.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M) {
    theta.sim[i,j]<-(Ybar.sim[i+1,j]-Ybar.sim[1,j])/sqrt(Vneyman.sim[i,j])
  }
}

pvalue.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M){
    pvalue.sim[i,j]=2*pnorm(-abs(theta.sim[i,j]))
  }
}

pcombine.sim<-as.numeric()
for(j in 1:M) {pcombine.sim[j]=mean(pvalue.sim[,j])}

pcombine.sim

}
```

```{r}
simulation.median<-function(N,K,n,Ybar,sigma,M){
  
Ybar.sim<-matrix(nrow = K,ncol = M)
variance.sim<-matrix(nrow = K,ncol = M)
for (i in 1:K){
for (j in 1:M){
  Ybar.sim[i,j]=rnorm(1,Ybar[i],sigma[i]*sqrt(K/N))
  variance.sim[i,j]=rchisq(1,n[i]-1)*(sigma[i]^2)/(n[i]-1)
}
}

Vneyman.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M) {
    Vneyman.sim[i,j]<-variance.sim[i+1,j]/n[i]+variance.sim[1,j]/n[1]
  }
}

theta.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M) {
    theta.sim[i,j]<-(Ybar.sim[i+1,j]-Ybar.sim[1,j])/sqrt(Vneyman.sim[i,j])
  }
}

pvalue.sim<-matrix(nrow = (K-1),ncol = M)
for(i in 1:(K-1)){
  for(j in 1:M){
    pvalue.sim[i,j]=2*pnorm(-abs(theta.sim[i,j]))
  }
}

pcombine.sim<-as.numeric()
for(j in 1:M) {pcombine.sim[j]=median(pvalue.sim[,j])}

pcombine.sim

}
```
```{r}
#Here,I choose the numbe of the arms,K, equals to 10. And under the alternative hypothesis, the last Ybar equals to 1.
K=10
i=1 
prob.min.5<-as.numeric()
prob.mean.5<-as.numeric()
prob.median.5<-as.numeric()
for(i in 1:100){
  pcombine.min<-simulation(100*i,K,rep(100*i/10,10),rep(0,K),rep(1,K),10000)
  cutoff.min<-sort(pcombine.min)[500]
  pcombine.min<-simulation(100*i,K,rep(100*i/10,10),c(rep(0,K-1),1),rep(1,K),10000)
  prob.min.5[i]<-length(which(pcombine.min<cutoff.min))/length(pcombine.min)
  
  pcombine.mean<-simulation.mean(100*i,K,rep(100*i/10,10),rep(0,K),rep(1,K),10000)
  cutoff.mean<-sort(pcombine.mean)[500]
  pcombine.mean<-simulation.mean(100*i,K,rep(100*i/10,10),c(rep(0,K-1),1),rep(1,K),10000)
  prob.mean.5[i]<-length(which(pcombine.mean<cutoff.mean))/length(pcombine.mean)
  
  pcombine.median<-simulation.median(100*i,K,rep(100*i/10,10),rep(0,K),rep(1,K),10000)
  cutoff.median<-sort(pcombine.median)[500]
  pcombine.median<-simulation.median(100*i,K,rep(100*i/10,10),c(rep(0,K-1),1),rep(1,K),10000)
  prob.median.5[i]<-length(which(pcombine.median<cutoff.min))/length(pcombine.median)
  
}
```
```{r}
plot(1:100,prob.min.5,xlab = "Sample size",ylab = "power",main = "Power versus Sample size")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,epsilon=1,method=minimum")
plot(1:100,prob.mean.5,xlab = "Sample size",ylab = "power",main = "Power versus Sample size")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,epsilon=1,method=mean")
plot(1:100,prob.median.5,xlab = "Sample size",ylab = "power",main = "Power versus Sample size")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,epsilon=1,method=median")

```

```{r}
prob.min.2<-as.numeric()
for(i in 100:1000){
  pcombine.min<-simulation(i,K,rep(i,10),rep(0,K),rep(1,K),10000)
  cutoff.min<-sort(pcombine.min)[500]
  pcombine.min<-simulation(N,K,rep(i,10),c(rep(0,K-1),1),rep(1,K),10000)
  prob.min.2[i]<-length(which(pcombine.min<cutoff.min))/length(pcombine.min)
  
}
```

```{r}
plot(seq(100,1000,1),prob.min.2[100:1000],xlab = "Sample size",ylab = "power",main = "Power versus Sample size")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,epsilon=1,method=minimum")
```

```{r}
K=4
i=1 
prob.min<-as.numeric()
prob.mean<-as.numeric()
prob.median<-as.numeric()
for(i in 1:100){
  pcombine.min<-simulation(100*i,K,rep(100*i/K,K),rep(0,K),rep(1,K),10000)
  cutoff.min<-sort(pcombine.min)[500]
  pcombine.min<-simulation(100*i,K,rep(100*i/K,K),c(rep(0,K-1),1),rep(1,K),10000)
  prob.min[i]<-length(which(pcombine.min<cutoff.min))/length(pcombine.min)
  
  pcombine.mean<-simulation.mean(100*i,K,rep(100*i/K,K),rep(0,K),rep(1,K),10000)
  cutoff.mean<-sort(pcombine.mean)[500]
  pcombine.mean<-simulation.mean(100*i,K,rep(100*i/K,K),c(rep(0,K-1),1),rep(1,K),10000)
  prob.mean[i]<-length(which(pcombine.mean<cutoff.mean))/length(pcombine.mean)
  
  pcombine.median<-simulation.median(100*i,K,rep(100*i/K,K),rep(0,K),rep(1,K),10000)
  cutoff.median<-sort(pcombine.median)[500]
  pcombine.median<-simulation.median(100*i,K,rep(100*i/K,K),c(rep(0,K-1),1),rep(1,K),10000)
  prob.median[i]<-length(which(pcombine.median<cutoff.min))/length(pcombine.median)
  
}
```
```{r}
#par(mfrow=c(3,1))
plot(1:100,prob.min,xlab = "Sample size",ylab = "power",main = "Power versus Sample size")
legend("bottomright","(x,y)", pch = 1,,title = "K=4,epsilon=1,method=minimum")
plot(1:100,prob.mean,xlab = "Sample size",ylab = "power",main = "Power versus Sample size")
legend("bottomright","(x,y)", pch = 1,,title = "K=4,epsilon=1,method=minimum")
plot(1:100,prob.median,xlab = "Sample size",ylab = "power",main = "Power versus Sample size")
legend("bottomright","(x,y)", pch = 1,,title = "K=4,epsilon=1,method=minimum")

```
From the figure above,I find that no matter the number of the arms is large or small(here I choose 4 and 10), the plot between the probability and the K has an increasing trend when using minimum to calculate pcombine. However, for mean and median to produce pcombine,the plot displays a random fluctuation. Transendentally, we expect a increasing trend for this plot,which means minimum behaves better than median and mean.


#task3:
different alternative types
##under a different alternative hypothesis, draw the usual statistics power plot(Ybar versus probability).
```{r}
epsilon<-seq(0.01,1,0.01)
Prob.min<-as.numeric()
for(i in 1:100){
  pcombine.min<-simulation(10000,10,rep(10,1000),rep(0,10),rep(1,10),10000)
  cutoff.min<-sort(pcombine.min)[500]
  y<-mvtnorm::rmvnorm(1,rep(1,9),diag(9))
  ynorm<-y*epsilon[i]/sqrt(sum(y^2))
  Ybar<-c(0,ynorm)
  pcombine.min<-simulation(10000,10,rep(10,1000),Ybar,rep(1,10),10000)
  Prob.min[i]<-length(which(pcombine.min<cutoff.min))/length(pcombine.min)
  
}

```
```{r}
plot(epsilon,Prob.min,xlab = "Module",ylab = "Power",main = "Power versus module under fixing-vector-module hypothesis ")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,method=minimum")
```

```{r}
epsilon<-seq(0.01,1,0.01)
Prob.mean<-as.numeric()
for(i in 1:100){
  pcombine.mean<-simulation.mean(10000,10,rep(10,1000),rep(0,10),rep(1,10),10000)
  cutoff.mean<-sort(pcombine.mean)[500]
  y<-mvtnorm::rmvnorm(1,rep(1,9),diag(9))
  ynorm<-y*epsilon[i]/sqrt(sum(y^2))
  Ybar<-c(0,ynorm)
  pcombine.mean<-simulation.mean(10000,10,rep(10,1000),Ybar,rep(1,10),10000)
  Prob.mean[i]<-length(which(pcombine.mean<cutoff.mean))/length(pcombine.mean)
  
}

epsilon<-seq(0.01,1,0.01)
Prob.median<-as.numeric()
for(i in 1:100){
  pcombine.median<-simulation.median(10000,10,rep(10,1000),rep(0,10),rep(1,10),10000)
  cutoff.median<-sort(pcombine.median)[500]
  y<-mvtnorm::rmvnorm(1,rep(1,9),diag(9))
  ynorm<-y*epsilon[i]/sqrt(sum(y^2))
  Ybar<-c(0,ynorm)
  pcombine.median<-simulation.median(10000,10,rep(10,1000),Ybar,rep(1,10),10000)
  Prob.median[i]<-length(which(pcombine.median<cutoff.median))/length(pcombine.median)
  
}

```
```{r}
plot(epsilon,Prob.mean,xlab = "Module",ylab = "Power",main = "Power versus module under fixing-vector-module hypothesis ")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,method=mean")
plot(epsilon,Prob.median,xlab = "Module",ylab = "Power",main = "Power versus module under fixing-vector-module hypothesis ")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,method=median")

```
Utilizing a different alternative hypothesis(the module of the vector Ybar equals to $epsilon$), I get similar "S" curves under three different ways to construct pcombine. However, all of these displays slightly fluctuation in the figures, which, I assume, comes from the randomness of generating Ybar vector.
Therefore, instead of plot probability directly versus $epsilon$, I use the mean of the probability coming from calculating the corresponding probability 10 times.


```{r}
epsilon<-seq(0.01,1,0.01)
Prob.min<-as.numeric()
for(i in 1:100){
  sum=0
  pcombine.min<-simulation(10000,10,rep(10,1000),rep(0,10),rep(1,10),10000)
  cutoff.min<-sort(pcombine.min)[500]
  for(j in 1:10){
  y<-mvtnorm::rmvnorm(1,rep(1,9),diag(9))
  ynorm<-y*epsilon[i]/sqrt(sum(y^2))
  Ybar<-c(0,ynorm)
  
  pcombine.min<-simulation(10000,10,rep(10,1000),Ybar,rep(1,10),10000)
  Prob<-length(which(pcombine.min<cutoff.min))/length(pcombine.min)
  sum=sum+Prob
  }
  Prob.min[i]=sum/10
}
  
epsilon<-seq(0.01,1,0.01)
Prob.mean<-as.numeric()
for(i in 1:100){
  sum=0
  pcombine.mean<-simulation.mean(10000,10,rep(10,1000),rep(0,10),rep(1,10),10000)
  cutoff.mean<-sort(pcombine.mean)[500]
  for(j in 1:10){
  y<-mvtnorm::rmvnorm(1,rep(1,9),diag(9))
  ynorm<-y*epsilon[i]/sqrt(sum(y^2))
  Ybar<-c(0,ynorm)
  pcombine.mean<-simulation.mean(10000,10,rep(10,1000),Ybar,rep(1,10),10000)
  Prob<-length(which(pcombine.mean<cutoff.mean))/length(pcombine.mean)
  sum=sum+Prob
  }
  Prob.mean[i]=sum/10
}

epsilon<-seq(0.01,1,0.01)
Prob.median<-as.numeric()
for(i in 1:100){
  sum=0
  pcombine.median<-simulation.median(10000,10,rep(10,1000),rep(0,10),rep(1,10),10000)
  cutoff.median<-sort(pcombine.median)[500]
  for(j in 1:10){
  y<-mvtnorm::rmvnorm(1,rep(1,9),diag(9))
  ynorm<-y*epsilon[i]/sqrt(sum(y^2))
  Ybar<-c(0,ynorm)
  
  pcombine.median<-simulation.median(10000,10,rep(10,1000),Ybar,rep(1,10),10000)
  Prob<-length(which(pcombine.median<cutoff.median))/length(pcombine.median)
  sum=sum+Prob
  }
  Prob.median[i]=sum/10
}

```
```{r}
plot(epsilon,Prob.min,xlab = "Module",ylab = "Power",main = "Power versus module under fixed-vector-module hypothesis")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,method=minimum,Using mean to remove randomness")
plot(epsilon,Prob.mean,xlab = "Module",ylab = "Power",main = "Power versus module under fixed-vector-module hypothesis") 
legend("bottomright","(x,y)", pch = 1,,title = "K=10,method=mean,Using mean to remove randomness")
plot(epsilon,Prob.median,xlab = "Module",ylab = "Power",main = "Power versus module under fixed-vector-module hypothesis")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,method=median,Using mean to remove randomness")
```
Here, I find that after removing part of randomness via using the mean of power *** is the best in this case since it has the most smooth plot.

##Under a different alternative hypothesis(fixing vector module), draw thw plot between the sample size N and the corresponding probability. K is a constant,10 in this case.
```{r}
K=10
i=1
y<-mvtnorm::rmvnorm(1,rep(1,9),diag(9))
ynorm<-y*1/sqrt(sum(y^2))
Ybar<-c(0,ynorm)
probability.min<-as.numeric()
probability.mean<-as.numeric()
probability.median<-as.numeric()
for(i in 1:100){
  pcombine.min<-simulation(100*i,K,rep(100*i/10,10),rep(0,K),rep(1,K),10000)
  cutoff.min<-sort(pcombine.min)[500]
  pcombine.min<-simulation(100*i,K,rep(100*i/10,10),c(0,Ybar),rep(1,K),10000)
  probability.min[i]<-length(which(pcombine.min<cutoff.min))/length(pcombine.min)
  
  pcombine.mean<-simulation.mean(100*i,K,rep(100*i/10,10),rep(0,K),rep(1,K),10000)
  cutoff.mean<-sort(pcombine.mean)[500]
  pcombine.mean<-simulation.mean(100*i,K,rep(100*i/10,10),c(0,Ybar),rep(1,K),10000)
  probability.mean[i]<-length(which(pcombine.mean<cutoff.mean))/length(pcombine.mean)
  
  pcombine.median<-simulation.median(100*i,K,rep(100*i/10,10),rep(0,K),rep(1,K),10000)
  cutoff.median<-sort(pcombine.median)[500]
  pcombine.median<-simulation.median(100*i,K,rep(100*i/10,10),c(0,Ybar),rep(1,K),10000)
  probability.median[i]<-length(which(pcombine.median<cutoff.median))/length(pcombine.median)
  
}
```

```{r}
#par(mfrow=c(3,1))
plot(1:100,probability.min,xlab = "Sample Size",ylab = "Power",main = "Power versus Sample Size under constant-vector-module hypothesis") 
legend("bottomright","(x,y)", pch = 1,,title = "K=10,epsilon=1,method=minimum")
plot(1:100,probability.median,xlab = "Sample",ylab = "Power",main = "Power versus Sample Size under constant-vector-module hypothesis")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,epsilon=1,method=mean")
plot(1:100,probability.mean,xlab = "Sample",ylab = "Power",main = "Power versus Sample Size under constant-vector-module hypothesis")
legend("bottomright","(x,y)", pch = 1,,title = "K=10,epsilon=1,method=median")
```