## Methods in Biostatistics 2 (140.652)

### Non-parametric test

```{r}
xA <- c(27.0,54.6,33.5,27.6,46.0,22.2,44.2,17.3,15.9,32.8)
xB <- c(17.4,20.5,13.9,14.8,27.9,10.6,33.7,15.4,25.0,24.1)

stripchart(list(xA,xB),vertical=T,pch=21,xlim=c(0.5,2.5))
segments(0.9,mean(xA),1.1,mean(xA),lwd=2,col="red")
segments(1.9,mean(xB),2.1,mean(xB),lwd=2,col="blue")
abline(h=mean(c(xA,xB)),lty=2)

t.test(xA,xB)
wilcox.test(xA,xB)
```

Let's see how one outlier influences the test statistics!

```{r}
sort(xA)
sort(xA,decreasing=T)
order(xA)
order(xA,decreasing=T)

xAnew <- xA
xAnew[2] <- 99.9
xA
xAnew

stripchart(list(xAnew,xB),vertical=T,pch=21,xlim=c(0.5,2.5))
segments(0.9,mean(xAnew),1.1,mean(xAnew),lwd=2,col="red")
segments(1.9,mean(xB),2.1,mean(xB),lwd=2,col="blue")
abline(h=mean(c(xAnew,xB)),lty=2)

t.test(xA,xB)$p.value
t.test(xAnew,xB)$p.value # why is the p-value now > 0.05?
wilcox.test(xA,xB)$p.value
wilcox.test(xAnew,xB)$p.value # why are the p-values the same?
```

### Permutation test

Let's write a function first, to make things easier. Note that the input are two vectors, generically called xa and xb.

```{r}
my.perm <- function(xa,xb){
  # combine the two vectots
  x <- c(xa,xb)
  # calculate the total length of observations
  n <- length(x)
  # calculate the length of the first vector
  na <- length(xa)
  # shuffle the order of the observations across grups
  z <- x[sample(1:n)]
  # return the t test statistic from the shuffled data
  return(t.test(z[1:na],z[-(1:na)],var.equal=T)$stat)
}
```

Example:

```{r}
# Let's try one with our xA and xB
my.perm(xA,xB)
# And onther one
my.perm(xA,xB)
# And onther one
my.perm(xA,xB)

# Let's do 1000 shuffles
nit <- 1000
# Create a vector of length nit (1000)
myres <- rep(NA,nit)
# Fill in the slots of that vector with the permutation test t statistics
for(j in 1:nit) myres[j] <- my.perm(xA,xB)

# Plot the data
hist(myres,breaks=21,col="lightgrey")
# Add the observed test statistic to the plot
tst <- t.test(xA,xB,var.equal=T)$stat 
abline(v=tst,col="red")
# Calculate the permutation p-value
mean(abs(myres)>abs(tst))
```

### End of code