Testing.Rmd

---
title: "R Notebook"
output: html_notebook
---

# Validity analysis 

## Simplest scenario: no prediction, unpaired t-tests 


Lets's start from the simplest scenario and build up from here. 
Sample twice 100 data points from the same Normal distribution, and perform unpaired ttest on the 2 vectors of sampled values. 
Repeat 10000 times and check if the distribution of pvalues is uniform.

```{r}
library(stats)
library(dplyr)

results=list()
n1=100
m=10000
for (i in 1:m){
  set.seed(i)
vec1=rnorm(n1, mean = 0, sd=1)
vec2=rnorm(n1, mean=0, sd=1)

test=t.test(vec1,vec2, alternative="two.sided", paired=FALSE)
results[[i]]=list()
results[[i]]=test
}
```

```{r}
source("get_pvalues.R")
res_pvalues=as.numeric(get_pvalues(results))
```
```{r}
hist(res_pvalues, col="cadetblue")
```

Now, do same thing but sampling from multivariate normal distribution with  2x2 Identity matrix as covariance matrix. 
Sample twice 100 data points for 2 variables (100x2 matrices) and 10 subjects.
For each matrix obtained compute sample covariance, select the same element from each cov matrix (here, [2,1]) and perform unpaired ttest on the 2 vectors of extracted values.  Repeat 10000 times and check for uniformity in the distribution of pvalues.

```{r}
library(MASS)
mu=replicate(2,0)
sig=diag(x=1,2,2)
ttests=list()
s1=10
for (j in 1:m){
  set.seed(j)
mat1 <-lapply(1:s1, mvrnorm, n=100, mu=mu, Sigma=sig)
mat2 <-lapply(1:s1, mvrnorm, n=100, mu=mu, Sigma=sig)

sc1=list()
sc2=list()
for (i in 1:s1){
  sc1[[i]]=cov(as.matrix(mat1[[i]]))
  sc2[[i]]=cov(as.matrix(mat2[[i]]))
}

v1=list()
v2=list()
for(i in 1:s1){
  v1[i]=sc1[[i]][2,1]
  v2[i]=sc2[[i]][2,1]
}
test_1=t.test(as.numeric(v1), as.numeric(v2), alternative = "two.sided", paired = FALSE)
ttests[[j]]=test_1
}
pvals1=get_pvalues(ttests)
```

```{r}
hist(as.numeric(pvals1) , col = "deeppink")
```

Exactly as above, but with higher dimensionality. Here, I sample from multivariate normal distribution with 15x15 Identity covariance matrix. 

```{r}
mu=replicate(15,0)
sig=diag(x=1,15,15)
ttests=list()
s2=10
for (j in 1:m){
  set.seed(j)
mat1 <-lapply(1:s2, mvrnorm, n=100, mu=mu, Sigma=sig)
mat2 <-lapply(1:s2, mvrnorm, n=100, mu=mu, Sigma=sig)

sc1=list()
sc2=list()
for (i in 1:s2){
  sc1[[i]]=cov(as.matrix(mat1[[i]]))
  sc2[[i]]=cov(as.matrix(mat2[[i]]))
}

v1=list()
v2=list()
for(i in 1:s2){
  v1[i]=sc1[[i]][2,1]
  v2[i]=sc2[[i]][2,1]
}
test_1=t.test(as.numeric(v1), as.numeric(v2), alternative = "two.sided", paired = FALSE)
ttests[[j]]=test_1
}
pvals=get_pvalues(ttests)
```

```{r}
hist(as.numeric(pvals), col="darksalmon")
```

```{r}
library(ggplot2)
library(reshape2)
df=data.frame(pvals1=as.numeric(res_pvalues), pvals2=as.numeric(pvals1), pvals3=as.numeric(pvals))
plot8= ggplot(melt(df),aes(value,fill=variable)) + geom_histogram(position="dodge", color="black", binwidth = 0.05)+
  #scale_x_discrete(labels = c("Drug", "Gas", "Gun", "Hang", "Jump", "Other")) +
  #scale_fill_discrete(labels=c("M", "F")) +
  facet_wrap( ~ variable, labeller = "label_value", dir="v")
plot8
```



## One step further: Use average covariance from opposite condition to predict the covariance matrix for the missing potential outcome.


Let's get back to the simple scenario of sampling from same multivariate normal distribution with 2x2 identity matrix.
n=100 s=10 cov=2x2 identity matrix. 

After finding the sample covariance matrix for each subject in each condition, find the average covariance matrix from the opposite condition and use it as prediction for the covariance matrix of the missing potential outcome. 
Build the treatment and control vectors, extract same matrix element from each matrix of the 2 vectors and perform paired ttest.
Repeat 10000 times. 
To see if using a correction when computing the average gives an unbiased estimator, I perform both with and without correction on the same data points.
How are the pvalues distributed?

```{r}
source("base_analysis.R")
source("sample_covariance.R")
source("get_triangle.R")
source("base_analysis_unbiased.R")
mu=replicate(2,0)
sig=diag(x=1,2,2)
ttests_1=list()
ttests_2=list()
s3=10
for (j in 1:m){
  set.seed(j)
mat1 <-lapply(1:s3, mvrnorm, n=100, mu=mu, Sigma=sig)
mat2 <-lapply(1:s3, mvrnorm, n=100, mu=mu, Sigma=sig)

analysis_1=base_analysis(mat1,mat2,channels=2, n_row=2, n_col=1,paired=TRUE)
analysis_2=base_analysis_unbiased(mat1,mat2,channels=2, n_row=2, n_col=1, paired=TRUE)
ttests_1[[j]]=analysis_1$test
ttests_2[[j]]=analysis_2$test

}
pvals=get_pvalues(ttests_1)
pvals_bis=get_pvalues(ttests_2)
```


```{r}
hist(as.numeric(pvals), breaks = seq(0,1,by=0.05), main="Histogram of pvalues, without adjustment term", col="gold")
#abline(h=qbinom(0.95, length(pvals), 0.05), col="red")
#abline(h=qbinom(v, length(pvals), 0.05), col="blue")
```
```{r}
hist(as.numeric(pvals_bis), breaks = seq(0,1,by=0.05), main="Histogram of pvalues, with adjustment term", col="gold")
```
```{r}
df2=data.frame("no adjustment"=as.numeric(pvals), "adjusted"=as.numeric(pvals_bis))
plot9=ggplot(melt(df2), aes(value,fill=variable))+
  geom_histogram(position = "identity", color="black", binwidth = 0.05, alpha=0.8)+
  facet_wrap( ~ variable, labeller = "label_value", dir="v")
plot9
```


How many rejections are there?
```{r}
pv_indexed=tibble(pv=pvals, index=seq(1,m))
pv_ind_ordered=pv_indexed%>%
  arrange(pv)
count(filter(pv_ind_ordered, pv_ind_ordered$pv<0.05))
```
```{r}
pv_indexed=tibble(pv=pvals_bis, index=seq(1,m))
pv_ind_ordered=pv_indexed%>%
  arrange(pv)
count(filter(pv_ind_ordered, pv_ind_ordered$pv<0.05))
```
BH adjustment
```{r}
i=seq(along=pvals)
i2=seq(along=pvals_bis)
m=10000
alpha=0.05
k2 <- max( which( sort(as.numeric(pvals)) < i/m*alpha) )
k3 <- max( which( sort(as.numeric(pvals_bis)) < i2/m*alpha) )
cutoff2 <- sort(as.numeric(pvals))[k2]
cutoff3 <- sort(as.numeric(pvals_bis))[k3]
cat("k =",k2,"p-value cutoff=",cutoff2)
print("")
cat("k =",k3,"p-value cutoff=",cutoff3)
```
```{r}
fdr <- p.adjust(pvals, method="fdr")
#mypar(1,1)
plot(pvals,fdr,log="xy")
abline(h=alpha,v=cutoff)
```



## Use average covariance from opposite condition, based on k nearest neighbors for the covariate, to predict the covariance matrix for the missing potential outcome.

```{r}
X_conv=as.matrix(X_conv_subset[,2])
X_short=as.matrix(X_short_subset[,2])
```

```{r}
library(RANN)
source("knn_analysis_unbiased.R")
source("knn_analysis.R")
mu=replicate(2,0)
sig=diag(x=1,2,2)
ttests_3=list()
ttests_4=list()
na=20
nb=30
M=10000

unique_X_s<-unique(X_short)%>%sort(decreasing = FALSE)
max=max(unique_X_s)
min=min(unique_X_s)
X_s<-data.frame(X=sample(seq(min,max,by=0.5),size=20,replace=TRUE))

unique_X_c<-unique(X_conv)%>%sort(decreasing = FALSE)
max=max(unique_X_c)
min=min(unique_X_c)
X_c<-data.frame(X=sample(seq(min,max,by=0.5),size=30,replace=TRUE))


for (j in 1:M){
  set.seed(j)
mat_short <-lapply(1:na, mvrnorm, n=100, mu=mu, Sigma=sig)
mat_conv <-lapply(1:nb, mvrnorm, n=100, mu=mu, Sigma=sig)

test_3=knn_analysis(X_s$X, X_c$X, mat_short, mat_conv, channels=2, n_row=2, n_col=1, K=5)
test_4=knn_analysis_unbiased(X_s$X, X_c$X, mat_short, mat_conv,channels=2,n_row=2, n_col=1, K=5)
ttests_3[[j]]=test_3$test
ttests_4[[j]]=test_4$test
}
pvals_3=get_pvalues(ttests_3)
pvals_4=get_pvalues(ttests_4)
``` 
```{r}
df3=data.frame("no adjustment"=as.numeric(pvals_3), "adjusted"=as.numeric(pvals_4))
plot10=ggplot(melt(df3), aes(value,fill=variable))+
  geom_histogram(position = "identity", color="black", binwidth = 0.05, alpha=0.8)+
  facet_wrap( ~ variable, labeller = "label_value", dir="v")
plot10
```


```{r}
hist(as.numeric(pvals_3), breaks = seq(0,1,by=0.05), main="Histogram of pvalues, knn without adjustment term", col="gold")
```
```{r}
v=1-0.5*0.5
hist(as.numeric(pvals_4), breaks = seq(0,1,by=0.05), main="Histogram of pvalues, knn with adjustment term", col="gold")
abline(h=qbinom(0.95, length(pvals_4), 0.05), col="red")
abline(h=qbinom(v, length(pvals_4), 0.05), col="blue")
```

```{r}
pv_indexed=tibble(pv=pvals_3, index=seq(1,m))
pv_ind_ordered=pv_indexed%>%
  arrange(pv)
count(filter(pv_ind_ordered, pv_ind_ordered$pv<0.05))
```
```{r}
pv_indexed=tibble(pv=pvals_4, index=seq(1,m))
pv_ind_ordered=pv_indexed%>%
  arrange(pv)
count(filter(pv_ind_ordered, pv_ind_ordered$pv<0.05))
```
```{r}
i3=seq(along=pvals_3)
#i4=seq(along=pvals_4)
m=10000
alpha=0.05
k4 <- max( which( sort(as.numeric(pvals_3)) < i3/m*alpha) )
k5 <- max( which( sort(as.numeric(pvals_4)) < i3/m*alpha) )
cutoff4 <- sort(as.numeric(pvals_3))[k4]
cutoff5<- sort(as.numeric(pvals_4))[k5]
cat("k =",k4,"p-value cutoff=",cutoff4)
print("")
cat("k =",k5,"p-value cutoff=",cutoff5)
```

```{r}
fdr <- p.adjust(pvals_3, method="fdr")
#mypar(1,1)
plot(pvals_3,fdr,log="xy")
abline(h=alpha,v=cutoff)
```
 



Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Ctrl+Alt+I*.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Ctrl+Shift+K* to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.