TreatmentSizeK.Rmd

---
title: "PowerKAnalysis"
author: "Xiaohui Li Yuxin Ma"
date: "6/4/2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(Iso)
source("PvalueFunction.R")
source("RatioPower.R")
```

```{r}
# Setting the unit number to 1000
N = 10000
```



# A. Test the power shape for effective size
## k = 10 
```{r}
K = 10
Y = rep(0, K)
Y_j = Y
power_10 = NULL
for (j in 1:100){
  Y_j[2] = 0.01 * j
  r1_j = Y_j[2] /1
  power_10[j] = ratio_power(Y_j, r1_j, 1, N, K, compute_pcombine, compute_power)

  }
plot(seq(0.01,1,0.01), power_10, type = "l", xlab = "r1_value", ylab ="power ", main = "Power with increasing r1 given K")

```

##k = 50
```{r}
K = 50
Y = rep(0, K)
Y_j = Y
power_50 = NULL
for (j in 1:100){
  Y_j[2] = 0.01 * j
  r1_j = Y_j[2] /1
  power_50[j] = ratio_power(Y_j, r1_j , 1, N, K, compute_pcombine, compute_power)
}
plot(seq(0.01,1,0.01), power_50, type = "l", xlab = "r1_value", main = "Power with effective Y1 given K = 50" )

```

## K = 100
```{r}
K = 100
Y = rep(0, K)
Y_j = Y
power_100 = NULL
for (j in 1:100){
  Y_j[2] = 0.01 * j
  r1_j = Y_j[2] /1
  power_100[j] = ratio_power(Y_j, r1_j , 1, N, K, compute_pcombine, compute_power)
}
plot(seq(0.01,1,0.01), power_100, type = "l", xlab = "r1_value", main = "Power with effective Y1 given K = 100" )




```

## Compare the power pattern given different size K
```{r}
plot(seq(0.01,1,0.01), power_10, type = "l", xlab = "r1_value", ylab ="power ", main = "Power with increasing r1 given different K", col = "turquoise" )
lines(seq(0.01,1,0.01), power_50, col = "blue")
lines(seq(0.01,1,0.01), power_100, col = "steelblue")
legend("bottomright", c("K = 10", "K = 50", "K = 100"), col = c("Turquoise","Blue","steelblue"), lty=1)
```


```{r}
K = 10
Y = rep(0, K)
Y_j = Y
power_10 = NULL
for (j in 1:100){
  Y_j[2] = 0.01 * j
  r1_j = Y_j[2] /1
  power_10[j] = ratio_power(Y_j, r1_j, 1, N, K, compute_pcombine, compute_power)

  }
plot(seq(0.01,1,0.01), power_10, type = "l", xlab = "r1_value", ylab ="power ", main = "Power with increasing r1 given K")
```






# B. Caculate the power with K increasing 
## r1 = 1
```{r}
power_K = NULL
for (k in 2:201){
  Y_k = rep(0,k)
  Y_k[2] = 1
  power_K[k] = ratio_power(Y_k , 1 , 1, N, k, compute_pcombine, compute_power)
}

plot(power_K, type = "l", xlab = "K size", main = "Power with k size given r1 = 1" )
```

## r1 = 2
```{r}
power_K = NULL
for (k in 2:200){
  Y_k = rep(0,k)
  Y_k[2] = 1
  power_K[k] = ratio_power(Y_k , 0.5 , 1, N, k, compute_pcombine, compute_power)
}

plot(seq(2,201,1), power_K, type = "l", xlab = "K size", main = "Power with k size given r1 = 0.5" )
```


# C. 80% Power critical ratio for treatment size K

```{r}
start.time = Sys.time()

critical_ratio = NULL
for (k in 2:100){
  Yk = rep(0,k)
  Yk[2] = 1 
  
  critical_ratio[k] = compute_binary_ctratio(Yk, N, k, 0.005)
}

# plot the critical ratio against the k treatments number
plot(critical_ratio, type = "l", xlab = "K size", main = "Critical ratio value with k size" )

end.time = Sys.time()
time.taken_binaryctratio = end.time - start.time
time.taken_binaryctratio

x <- seq(2,98,1)
critical_ratio = critical_ratio[1:97]
iso1 =  ufit(critical_ratio, x, type = "b")
plot(x, critical_ratio, type = "l", xlab = "K size", ylab = "Critical ratio", main = "Isotonic Regression Plot")
lines(iso1$y, col="red")



```


## Analyze the increasing rate of critical ratio with K - log K and squared K form
### Plot the critical ratio against log k and squared k
```{r}
logk = NULL
sqrtk = NULL
for (k in 2:100){
  logk[k] = log(k)
  sqrtk[k] = sqrt(k)
}
# Plot the critical ratio against log k 
  plot(logk, critical_ratio, type = "l", xlab = "log K size", main = "Critical ratio value with log k" )

# Plot the critical ratio against squared k
plot(sqrtk, critical_ratio, type = "l", xlab = "squared K size", main = "Critical ratio value with squared k" )
```

### Linear regression to test 
```{r}
logreg = lm(critical_ratio~logk)
summary(logreg)

sqrtreg = lm(critical_ratio ~ sqrtk)
summary(sqrtreg)
```