---
title: "Vitamin E and IQ - confirmatory study"
author: Jana Hozzová
output:
  pdf_document: default
  html_notebook: default
---

# Task

This document will guide you through the process of pre-registration and subsequent analysis. It's a part of exercise in Experiments seminar dedicated to quantitative evaluation of experiments. We will work with a made-up example of influence of vitamin E on IQ.

First, you need to answer the questions in pre-registration. Some options are pre-filled for you, others you can decide absolutely arbitrarily. Modify the scripts according to chosen options.

After you pre-registration is done, export it to pdf (click Knit). Then take data files with experiments' results and do your analysis. There are three folders, control, treatment and treatment-high. There are five control data files and 15x5 treatment files. This time, take treatment file with the first number corresponding with the day you were born modulo 15 (e.g. 28th \% 15 = 13). The second number in the filename goes from 1 to 5, you can use those if you are trying to replicate your study. After than is done, you can again export the whole notebook to pdf. 

If you have time and want to reproduce your experiments with different dataset, or explore different hypotheses, or try out different analyses, fill free to do so just start fresh with empty template (or copy your notebook to a NEW file), ask for a new dataset and do it all over.

*Text in italics are instructions and examples*, normal text is what is pre-filled for you.

This document is written in R markdown, for details see [here] (https://blog.rstudio.com/2014/08/01/the-r-markdown-cheat-sheet/).

# Pre-registration

## What's the main question being asked or hypothesis being tested in this study?

*Decide for yourself. Some examples of hypotheses:*

- *Vitamin E influences IQ.*
- *Vitamin E does not influence IQ.*
- *Vitamin E increases IQ.*
- *Vitamin E does not increase IQ.*
- *Vitamin E increases IQ by at least two IQ points.*
- *Vitamin E decreases IQ by at least two IQ points.*
- *Vitamin E increases IQ, but less than two points.*
- *Vitamin E increases IQ by five points exactly.*

*Explicitly write both your H0 and H1.*

## Describe the key dependent variable(s) specifying how they will be measured.

Dependent variable is IQ, measured by a Wechsler Adult Intelligence Scale test in controlled environment.

## How many and which conditions will participants be assigned to?

*Choose one of two following options or think up your own:*

*Two conditions: control group and treatment group (100% of recommended daily intake take in form of a daily pill from company GoodVitamins for 30 days).*

*Three conditions: control group, treatment group with low doses of vitamin E (100% of recommended daily intake taken in form of a daily pill from company GoodVitamins for 30 days) and treatment group with high doses of vitamin E (300% of recommended daily intake, taken in form of three daily pills from GoodVitamins for 30 days).*

## Specify exactly which analyses you will conduct to examine the main question/hypothesis.

*You can use Welch's t-test for basic testing, but feel free to use something else if you find it more suitable to test your hypothesis. In any case specify the statistical test you want to use and its parameters, e.g. direction of effect, $\alpha$ levels, etc.*

## Any secondary hypotheses and analyses?

*Go for it, come up with another hypothesis that you would find interesting to examine. For example:*

- *Vitamin E influences IQ more with higher dose.*
- *Vitamin E does not influence IQ with lower dose, but it does with higher dose.*

## How many observations will be collected or what will determine sample size? No need to justify decision, but be precise about exactly how the number will be determined.

*We were talking about this. You can determine the sample by:*

- *available resources (e.g. you have funding only to test 200 people together)*
- *estimated effect size and subsequent power analysis (you can pretend you estimate the effect size on other research or try and determine it in exploratory analysis beforehand - then you need to fill out exploratory template and do the analysis there before continuing this)*
- *the smallest effect size of interest and subsequent power analysis (you can argue that you are not interested in effects smaller than two IQ points)*

```{r}
if(!require(pwr)){install.packages('pwr')} #install necessary package if not available
library(pwr) #load the necessary package
wantedStatisticalPower <- 0.8 #assign the number between <0,1>
points <- 2 #estimated effect size in IQ points,
            #if you expect lower IQ in treatment, if should be a negative number
SD <- 15 # standard deviaiton of IQ in population
estimatedEffectSize <- points/SD # estimated effect size in Cohen's d
alpha <- 0.05 # Type 1 error rate
direction <- "greater" # or "two.sided" or "less" depending on your hypothesis
sampleSize <- pwr.t.test(d=estimatedEffectSize,sig.level=alpha,type="two.sample",
                         alternative=direction, power=wantedStatisticalPower)$n

sampleSize <- ceiling(sampleSize) #round it up
```

```{r}
#if you have determined your sample size in any other way, reassign the value here
sampleSize <- sampleSize
```


*In any case, report here the estimated statistical power of your study in a form of power curve.*

```{r, echo=FALSE}
N<-sampleSize
d<-estimatedEffectSize
p_upper <- alpha
direction<-direction
#ncp<-(input$d*sqrt(N/2)) #Calculate non-centrality parameter d
plot_power_d <- (function(d, N, p_upper, direction)
{
   power<-pwr.t.test(d=d, n=N,sig.level=p_upper,type="two.sample",
                     alternative=direction)$power
}
)
par(bg = "aliceblue")
plot(-10,xlab=substitute(paste("Cohen's ", delta)), ylab="Power", axes=FALSE,
main=(paste("Power for t-test with N = ",ceiling(N)," per group")), xlim=c(0,2),  ylim=c(0, 1))
abline(v = seq(0,2, 0.2), h = seq(0,1,0.1), col = "lightgray", lty = 1)
axis(side=1, at=seq(0,2, 0.2), labels=seq(0,2,0.2))
axis(side=2, at=seq(0,1, 0.2), labels=seq(0,1,0.2))
curve(plot_power_d(d=x, N=N, p_upper=p_upper, direction=direction), 0, 2, type="l", lty=1, lwd=2, ylim=c(0,1), xlim=c(0,2), add=TRUE)
points(x=d, y=plot_power_d(d, N, p_upper, direction=direction), cex=2, pch=19, col=rgb(1, 0, 0,0.5))
```


## How will you treat data exclusions and outliers?

If anybody does not finish the IQ test, we omit their data from  the study and they will not be included in studied dataset. If anybody in treatment group stops taking vitamin E before one-month period, we omit their data from the study and they will not be included in studied dataset. By outliers, we mean results outside 2.5 times standard deviation. These will be excluded from the datasets.

## Will you use sequential analysis (gradual acquisition of data) or test the data in multiple ways? How will you control your false positive rate?

*When you feel confident enough, try at least multiple testing in some form. You will encounter these techniques rather often in your own research. Some examples:*

- *We will collect data with sequential analysis. First, we will gather data from 400 people per group at first and test. If we are not able to reject H0 at that point, we will gather data from addidtional 450 people per group and test. Our total Type 1 error rate is set to 0.05, and we will control it with Bonferroni correction/Pocock boundary/other method.*

- *We will test the data two times. First, we will use Welch's t-test to compare control and treatment group. Then we will compare people below 40 in control and treatment groups. Our total Type 1 error rate is set to 0.07, and we control it with Bonferroni correction/Holm correction.*

- *We will test the data multiple times. We will divide both groups in two parts according to the following: age (below and above 40), socio-economic status (income below 150000 Kč per year per person in household and income above that), achieved education (university degree of some sort or not). For each division, we will compare corresponding control and treatment groups (i.e. control below 40 with treatment below 40, control above 40 and treatment above 40). That means 12 p-values. We will control false discovery rate with Benjamini–Hochberg procedure and our controlled $\alpha = 0.17$.*

## Anything else you would like to pre-register?

*Is there anything else to mention and prevent reviewers to think you misused the statistics?*

## Have any data been collected for this study already?

No, no data have been collected for this study yet.

# Analysis

```{r, include=FALSE}
# COMMENT THIS CHUNK OUT WHEN YOU ARE DOING YOUR ANALYSIS
# OTHERWISE THE CODES BELOW WILL NOT EXECUTE
knitr::opts_chunk$set(eval = FALSE)
```

*Set the sample size and load data.*
```{r}
sampleSize
controlFile <- "" #change to the actual file name of control group data
treatmentFile <- "" #change to the actual file name of treatment group data

# read the data, omitting the first line, reading samplesize number of data points
control <- read.csv(controlFile, skip=1, nrows = sampleSize+1)
control <- control[,1]
treatment <- read.csv(treatmentFile, skip=1, nrows = sampleSize+1)
treatment <- treatment[,1]
```

*Exclude the outliers.*
```{r}
controlMean <- mean(control)
controlSD <- 15
treatmentMean <- mean(treatment)
treatmentSD <- sd(treatment)
control <- control[control > controlMean-2.5*controlSD & 
                     control < controlMean+2.5*controlSD]
treatment <- treatment[treatment > treatmentMean - 2.5*treatmentSD &
                         treatment < treatmentMean + 2.5*treatmentSD]
```

*Do the statistical test. Control your Type 1 error rate.*
```{r}
# for t-test, you don't need any special library, just type ?t.test into console for help
# for testing "bounded" hypotheses, i.e. vitamin E increases IQ by more than two points, 
#     you can use package TOSTER
# for Bonferroni correction, just divide your alpha 
#     a_i <- alpha/number_of_tests and use a_i in tests
# for FDR, there is library(astsa) and its function FDR, 
#     look into real-world example notebook for details
# for other tests and corrections, use google, 
#     R has packages for almost everything

#beware that the order of parameters control and treatment matters
#if your direction is "greater", that means you expect control > treatment 
#   and your estimated effect size is a positive number
#if your direction is "less", that means you expect control < treatment 
#   and your estimated effect size is a negative number
# if your direction is "two.sided", you made no claims 
#   about the direction of the effect, 
#   you just think the treatment and control are not the same
mean(control)
mean(treatment)
test <- t.test(control, y = treatment, type = "two.sample", alternative = direction)
test$p.value
test$p.value < alpha
```


*Do the other analyses.*
```{r}
# just do what you promised in pre-registration
```

  
# Interpretation

## Interpret the results.

- *Are you able to reject H0?*
- *If so, what that means?*
- *If not, what that means?*
- *What about the effect size? What's its practical significance?*


## What would be your next steps?

- *hypotheses to examine*
- *tests to conduct *
- *exploratory analyses to do*
- *limitations of your study (e.g. small sample size)*
- *modifications to the future study setup (e.g. including possible confounding variables)*

# Optional Exploratory Analysis of Available Data
*Almost always, something will surprise you in your evaluation and you haven't even thought about that so it's not in your pre-registration. It's okay, we can do a little exploratory analysis with it, as long as we don't consider our hypotheses generated from that to be true. But you can, for example, make your pre-registration more detailed when you are reproducing your study or doing something similar next time.*

*So, if you have anything you would like to examine further in data you have available, you can do so here.*