= Classical Tests =
== Exercises ==
 * load the data (file: session4data.rdata)
 * make a new summary data frame (per subject and time) containing:
    * the number of trials
    * the number correct trials (absolute and relative)
    * the mean TTime and the standard deviation of TTime
    * the respective standard error of the mean
 * keep the information about Sex and Age\_PRETEST
 * make a plot with time on the x-axis and TTime on the y-axis showing the means and the 95\% confidence intervals (geom\_pointrange())
 * add the number of trials and the percentage of correct ones using geom\_text()
== Exercises - Solutions ==
   * load the data (file: session4data.rdata)
   * make a new summary data frame (per subject and time) containing:
{{{#!highlight r
> sumdf <- data %>%
+     group_by(Subject,Sex,Age_PRETEST,testid) %>%
+     summarise(count=n(),
+               n.corr = sum(Stim.Type=="hit"),
+               perc.corr = n.corr/count,
+               mean.ttime = mean(TTime),
+               sd.ttime = sd(TTime),
+               se.ttime = sd.ttime/sqrt(count))
> head(sumdf)
Subject Sex Age_PRETEST testid count n.corr perc.corr mean.ttime 
1       1   f        3.11  test1    95     63 0.6631579   8621.674 
2       1   f        3.11      1    60     32 0.5333333   9256.367 
3       1   f        3.11      2    59     32 0.5423729   9704.712 
4       1   f        3.11      3    60     38 0.6333333  14189.550 
5       1   f        3.11      4    59     31 0.5254237  13049.831 
6       1   f        3.11      5    59     33 0.5593220  14673.525 
Variables not shown: sd.ttime se.ttime (dbl)  
}}}
== four possible situations ==
|| ||Situation || ||
|| ||H_0 is true ||H_0 is false ||
||Conclusion ||H_0 is not rejected || Correct decision ||Type II error ||
||||H_0 is rejected ||Type I error || Correct decision ||

== Common symbols ==
||n ||number of observations (sample size)||
||K ||number of samples (each having n elements)||
||s ||standard deviation (sample)||
||r ||sample correlation coefficient||
||Z ||standard normal deviate ||
== Alternatives ==
<img alt='sesssion2/twosided.pdf}' src='-1' />
<img alt='sesssion2/greater.pdf}' src='-1' />
<img alt='sesssion2/less.pdf}' src='-1' />
== Alternatives ==
The p-value is the probability of the sample estimate (of the respective estimator) under the null.
== Z-test for a population mean ==
The z-test is a something like a t-test (it is like you would know almost everything about the perfect conditions. It uses the normal distribution as test statistic and is therefore a good example.
To investigate the significance of the difference between an assumed population mean {{{#!latex \mu_0}}} and a sample mean {{{#!latex \bar{x}}}}.
* It is necessary that the population variance {{{#!latex \sigma^2}}} is known. 
* The test is accurate if the population is normally distributed. If the population is not normal, the test will still give an approximate guide.
== Z-test for a population mean ==
* Write a function which takes a vector, the population standard deviation and the population mean as arguments and which gives the Z score as result. 
   * name the function ztest or my.z.test - not z.test because z.test is already used
   * set a default value for the population mean
* add a line to your function that allows you to process numeric vectors containing missing values!
* the function pnorm(Z) gives the probability of {{{#!latex x \leq Z {{{#!latex . Change your function so that it has the p-value (for a two sided test) as result. 
* now let the result be a named vector containing the estimated difference, Z, p and the n.
You can always test your function using simulated values: rnorm(100,mean=0) gives you a vector containing 100 normal distributed values with mean 0.
== Z-test for a population mean ==
Write a function which takes a vector, the population standard deviation and the population mean as arguments and which gives the Z score as result.
{{{#!highlight r
> ztest <- function(x,x.sd,mu=0){
+     sqrt(length(x)) * (mean(x)-mu)/x.sd
+ }
> set.seed(1)
> ztest(rnorm(100),x.sd = 1)
[1] 1.088874
}}}
== Z-test for a population mean ==
Add a line to your function that allows you to also process numeric vectors containing missing values!
{{{#!highlight r
> ztest <- function(x,x.sd,mu=0){
+     x <- x[!is.na(x)]
+     if(length(x) < 3) stop("too few values in x")
+     sqrt(length(x)) * (mean(x)-mu)/x.sd
+ }
}}}
== Z-test for a population mean ==
The function pnorm(Z) gives the probability of {{{#!latex x \leq Z {{{#!latex . Change your function so that it has the p-value (for a two sided test) as result. 
{{{#!highlight r
> ztest <- function(x,x.sd,mu=0){
+     x <- x[!is.na(x)]
+     if(length(x) < 3) stop("too few values in x")
+     z <- sqrt(length(x)) * (mean(x)-mu)/x.sd
+     2*pnorm(-abs(z))
+ }
> set.seed(1)
> ztest(rnorm(100),x.sd = 1)
[1] 0.2762096
}}}
== Z-test for a population mean ==
Now let the result be a named vector containing the estimated difference, Z, p and the n.
{{{#!highlight r
> ztest <- function(x,x.sd,mu=0){
+     x <- x[!is.na(x)]
+     if(length(x) < 3) stop("too few values in x")
+     est.diff <- mean(x)-mu
+     z <- sqrt(length(x)) * (est.diff)/x.sd
+     round(c(diff=est.diff,Z=z,pval=2*pnorm(-abs(z)),n=length(x)),4)
+ }
> set.seed(1)
> ztest(rnorm(100),x.sd = 1)
diff        Z     pval        n 
}}}
== Z-test for a population mean ==
* Z-test for two population means (variances known and equal)
* Z-test for two population means (variances known and unequal)
To investigate the statistical significance of the difference between an assumed population mean {{{#!latex \mu_0}}} and a sample mean {{{#!latex \bar{x}}}}. There is a function z.test() in the BSDA package.
* It is necessary that the population variance {{{#!latex \sigma^2}}} is known. 
* The test is accurate if the population is normally distributed. If the population is not normal, the test will still give an approximate guide.
== Simulation Exercises ==
* Now sample 100 values from a Normal distribution with mean 10 and standard deviation 2 and use a z-test to compare it against the population mean 10. What is the p-value?
* Now do the sampling and the testing 1000 times, what would be the number of statistically significant results? Use replicate() (which is a wrapper of tapply()) or a for() loop! Record at least the p-values and the estimated differences! Use table() to count the p-vals below 0.05. What type of error do you associate with it? What is the smallest absolute difference with a p-value below 0.05?
* Repeat the simulation above, change the sample size to 1000 in each of the 1000 samples! How many p-values below 0.05? What is now the smallest absolute difference with a p-value below 0.05?
== Simulation Exercises -- Solutions ==
   * Now sample 100 values from a Normal distribution with mean 10 and standard deviation 2 and use a z-test to compare it against the population mean 10. What is the p-value? What the estimated difference?
{{{#!highlight r
> ztest(rnorm(100,mean=10,sd=2),x.sd=2,mu=10)["pval"]
pval 
> ztest(rnorm(100,mean=10,sd=2),x.sd=2,mu=10)["diff"]
diff 
> ztest(rnorm(100,mean=10,sd=2),x.sd=2,mu=10)[c("pval","diff")]
pval   diff 
}}}
== Simulation Exercises -- Solutions ==
   * Now do the sampling and the testing 1000 times, what would be the number of statistically significant results? Use replicate() (which is a wrapper of tapply()) or a for() loop. Record at least the p-values and the estimated differences! Transform the result into a data frame.
using replicate()\footnotesize
{{{#!highlight r
> res <- replicate(1000, ztest(rnorm(100,mean=10,sd=2),x.sd=2,mu=10))
> res <- as.data.frame(t(res))
> head(res)
diff       Z   pval   n
1 -0.2834 -1.4170 0.1565 100
2  0.2540  1.2698 0.2042 100
3 -0.1915 -0.9576 0.3383 100
4  0.1462  0.7312 0.4646 100
5  0.1122  0.5612 0.5747 100
6 -0.0141 -0.0706 0.9437 100
}}}
== Simulation Exercises -- Solutions ==
   * Now do the sampling and the testing 1000 times, what would be the number of statistically significant results? Use replicate() (which is a wrapper of tapply()) or a for() loop. Record at least the p-values and the estimated differences! Transform the result into a data frame.
using replicate() II\footnotesize
{{{#!highlight r
> res <- replicate(1000, ztest(rnorm(100,mean=10,sd=2),x.sd=2,mu=10),
+                        simplify = F)
> res <- as.data.frame(Reduce(rbind,res))
> head(res)
diff       Z   pval   n
init -0.0175 -0.0874 0.9304 100
}}}
== Simulation Exercises -- Solutions ==
   * Now do the sampling and the testing 1000 times, what would be the number of statistically significant results? Use replicate() (which is a wrapper of tapply()) or a for() loop. Record at least the p-values and the estimated differences! Transform the result into a data frame.
using for() \scriptsize
{{{#!highlight r
> res <- matrix(numeric(2000),ncol=2)
> for(i in seq.int(1000)){
+     res[i,] <- ztest(rnorm(100,mean=10,sd=2),x.sd=2,mu=10)[c("pval","diff")] }
> res <- as.data.frame(res)
> names(res) <- c("pval","diff")
> head(res)
pval    diff
1 0.0591 -0.3775
2 0.2466  0.2317
3 0.6368  0.0944
4 0.5538 -0.1184
5 0.9897 -0.0026
6 0.7748  0.0572
}}}
== Simulation Exercises -- Solutions ==
   * Use table() to count the p-vals below 0.05. What type of error do you associate with it? What is the smallest absolute difference with a p-value below 0.05?
{{{#!highlight r
> table(res$pval < 0.05)
FALSE  TRUE 
960    40 
> tapply(abs(res$diff),res$pval < 0.05,summary)
> min(abs(res$diff[res$pval<0.05]))
[1] 0.3928
}}}
== Simulation Exercises -- Solutions ==
   * Repeat the simulation above, change the sample size to 1000 in each of the 1000 samples! How many p-values below 0.05? What is now the smallest absolute difference with a p-value below 0.05?
{{{#!highlight r
> res2 <- replicate(1000, ztest(rnorm(1000,mean=10,sd=2),
+                              x.sd=2,mu=10))
> res2 <- as.data.frame(t(res2))
> head(res2)
diff       Z   pval    n
1 -0.0731 -1.1559 0.2477 1000
2  0.0018  0.0292 0.9767 1000
3  0.0072  0.1144 0.9089 1000
4 -0.1145 -1.8100 0.0703 1000
5 -0.1719 -2.7183 0.0066 1000
6  0.0880  1.3916 0.1640 1000
}}}
== Simulation Exercises -- Solutions ==
   * Repeat the simulation above, change the sample size to 1000 in each of the 1000 samples! How many p-values below 0.05? What is now the smallest absolute difference with a p-value below 0.05?
{{{#!highlight r
> table(res2$pval < 0.05)
FALSE  TRUE 
946    54 
> tapply(abs(res2$diff),res$pval < 0.05,summary)
}}}
== Simulation Exercises Part II ==
* Concatenate the both resulting data frames from above using rbind()
* Plot the distributions of the pvals and the difference per sample size. Use ggplot2 with an appropriate geom (density/histogram)
* What is the message?
== Simulation Exercises -- Solutions ==
   * Concatenate the both resulting data frames from above using rbind()
   * Plot the distributions of the pvals and the difference per sample size. Use ggplot2 with an appropriate geom (density/histogram)
{{{#!highlight r
> res <- rbind(res,res2)  
> require(ggplot2)
> ggplot(res,aes(x=pval)) +
+     geom_histogram(bin=0.1,fill="forestgreen") +
+     facet_grid(~ n)
> ggsave("hist.png")
}}}
== Simulation Exercises -- Solutions ==
<img alt='sesssion2/hist.png' src='-1' />
== Simulation Exercises -- Solutions ==
   * Plot the distributions of the pvals and the difference per sample size. Use ggplot2 with an appropriate geom (density/histogram)
{{{#!highlight r
> ggplot(res,aes(x=diff,colour=factor(n))) +
+     geom_density(size=3)
> ggsave("dens.png")
}}}
== Simulation Exercises -- Solutions ==
<img alt='sesssion2/dens.png' src='-1' />
== Simulation Exercises -- Solutions ==
<img alt='sesssion2/point.png' src='-1' />
== Simulation Exercises -- Solutions ==
<img alt='sesssion2/dens2d.png' src='-1' />
{{{#!highlight r
> set.seed(1)
> x <- rnorm(12)
> t.test(x,mu=0) ## population mean 0
t = 1.1478, df = 11, p-value = 0.2754
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
sample estimates:
mean of x 
}}}
{{{#!highlight r
> t.test(x,mu=1) ## population mean 1
t = -3.125, df = 11, p-value = 0.009664
alternative hypothesis: true mean is not equal to 1
95 percent confidence interval:
sample estimates:
mean of x 
}}}
== t-tests ==
A t-test is any statistical hypothesis test in which the test statistic follows a \emph{Student's t distribution} if the null hypothesis is supported.
   * \emph{one sample t-test}: test a sample mean against a population mean
== One Sample t-test ==
== Two Sample t-tests ==
There are two ways to perform a two sample t-test in R:
   * given two vectors x and y containing the measurement values from the respective groups t.test(x,y)
   * given one vector x containing all the measurement values and one vector g containing the group membership {{{#!latex t.test(x \sim g)}}} (read: x dependend on g)
== Two Sample t-tests: two vector syntax ==
{{{#!highlight r
> set.seed(1)
> x <- rnorm(12)
> y <- rnorm(12)
> g <- sample(c("A","B"),12,replace = T)
> t.test(x,y)
t = 0.5939, df = 20.012, p-value = 0.5592
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
sample estimates:
mean of x  mean of y 
}}}
== Two Sample t-tests: formula syntax ==
{{{#!highlight r
> t.test(x ~ g)
t = -0.6644, df = 6.352, p-value = 0.5298
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
sample estimates:
mean in group A mean in group B 
}}}
== Welch/Satterthwaite vs. Student ==
   * if not stated otherwise t.test() will not assume that the variances in the both groups are equal
   * if one knows that both populations have the same variance set the var.equal argument to TRUE to perform a student's t-test
== Student's t-test ==
{{{#!highlight r
> t.test(x, y, var.equal = T)
t = 0.5939, df = 22, p-value = 0.5586
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
sample estimates:
mean of x  mean of y 
}}}
== t-test ==
   * the t-test, especially the Welch test is appropriate whenever the values are normally distributed
   * it is also recommended for group sizes {{{#!latex \geq 30}}} (robust against deviation from normality)
== Exercises ==
* use a t-test to compare TTime according to Stim.Type, visualize it. What is the problem?
* now do the same for Subject 1 on pre and post test (use filter() or indexing to get the resp. subsets)
* use the following code to do the test on every subset Subject and testid, try to figure what is happening in each step:\tiny
{{{#!highlight r
tob <- t.test(x$TTime ~ x$Stim.Type)
tmp <- data.frame(
Subject = unique(x$Subject),
testid = unique(x$testid),
pval = tob$p.value,
alternative = tob$alternative,
res <- Reduce(rbind,tmp.l)
}}}
* make plots to visualize the results. 
* how many tests have a statistically significant result? How many did you expect? Is there a tendency? What could be the next step?
== Exercises - Solutions ==
   * use a t-test to compare TTime according to Stim.Type, visualize it. What is the problem?
{{{#!highlight r
> t.test(data$TTime ~ data$Stim.Type)
t = -6.3567, df = 9541.891, p-value = 2.156e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
sample estimates:
mean in group hit mean in group incorrect 
> ggplot(data,aes(x=Stim.Type,y=TTime)) +
+     geom_boxplot()
}}}
== Exercises - Solutions ==
   * now do the same for Subject 1 on pre and post test (use filter() or indexing to get the resp. subsets)
{{{#!highlight r
> t.test(data$TTime[data$Subject==1 & data$testid=="test1"] ~
+        data$Stim.Type[data$Subject==1 & data$testid=="test1"])
t = -0.5846, df = 44.183, p-value = 0.5618
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
sample estimates:
mean in group hit mean in group incorrect 
> t.test(data$TTime[data$Subject==1 & data$testid=="test2"] ~
+        data$Stim.Type[data$Subject==1 & data$testid=="test2"])
t = -1.7694, df = 47.022, p-value = 0.08332
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
sample estimates:
mean in group hit mean in group incorrect 
}}}
== Exercises - Solutions ==
   * make plots to visualize the results
{{{#!highlight r
> ggplot(data,aes(x=testid,y=TTime)) +
+    geom_boxplot(aes(fill=Stim.Type)) +
+    facet_wrap(~Subject)
> ggplot(data,aes(x=factor(Subject),y=TTime)) +
+    geom_boxplot(aes(fill=Stim.Type)) +
+    facet_wrap(~testid)
}}}
== Exercises - Solutions ==
   * how many tests have an statistically significant result? How many did you expect? 
{{{#!highlight r
> table(res$pval < 0.05)
FALSE  TRUE 
165    22 
> prop.table(table(res$pval < 0.05))
FALSE      TRUE 
> 
}}}
== Exercises - Solutions ==
   *  What could be the next step?
{{{#!highlight r
tob <- t.test(x$TTime ~ x$Stim.Type)
tmp <- data.frame(
Subject = unique(x$Subject),
testid = unique(x$testid),
pval = tob$p.value,
alternative = tob$alternative,
res <- Reduce(rbind,tmp.l)
}}}