Session03

stats03.R

# CHAPTER 6

# TWO SAMPLES

# COMPARING TWO VARIANCES 

# To test to see whether the sample variances are significantly different we use Fisher's F test
# You divide the larger variance by the smaller variance
# We need to look for the critical value of the variance ratio 'qf' (quantiles of the F distribution)

qf(0.975,9,9)
# this gives us 4.025994 which means that a calculated variance ratio needs to be greater than this value

f.test.data<-read.csv("P:\\MSc\\Stats\\Data\\f.test.data.csv")
attach(f.test.data)
names(f.test.data)

var(gardenB)
var(gardenC)
# this computes the two variances
# the larger variance is in GardenC

F.ratio<-var(gardenC)/var(gardenB)
# this computes the F ratio
# it shows that th variance in gardenC is more than 10 times as big as the variance in gardenB

# the critical value for F is 4.026
# since the test statistic is larger than the cirtical value, we reject the null hypothesis
# we now accept the alternative hypothesis that the two variances are significantly different

# it is better to represent the P values associated to the F ratio: 'pf'

2*(1-pf(F.ratio,9,9))
# 0.0016
# so the probability of obtaining an F ratio as large as this or larger is less than 0.002
# this not the probability that the null hypothesis is true
# the null hypothesis is assumed to be true in carrying out the test

var.test(gardenB,gardenC)
# this speeds up the procedure
# although this command will but the variable that comes first in the alphabet on top, the p value remains the same
 detach(f.test.data)
 
# COMPARING TWO MEANS

# we calculate a test statistic, we judge the probability by comparing this test statistic with a critical value
# the critical value is calculated on the assumption that the null hypothesis is true
# there are two simple tests for comparing two sample means:

# 1. STUDENT'S T TEST

# the t distribution is a family of curves in which the number of degrees of freedom specifies a particluar curve
# we first formulate a nul hypothesis, stating that there is not effective difference between the observed sample mean and the hypothesiised mean
# a t-test can be two-sided (means are not equivalent) or one-sided (whether the the observed mean is larger or smaller than the hypothesised mean
# so if we don't know which element is going to have the higher mean it is a two-tailed test
# the critical value is:

qt(9.75,18)
# 2.101
# this means our test statistic needs to be higher than 2.1

t.test.data<-read.csv("P:\\MSc\\Stats\\Data\\t.test.data.csv")
attach(t.test.data)
names(t.test.data)
#"gardenA" "gardenB"

ozone<-c(gardenA,gardenB)
label<-factor(c(rep("A",10),rep("B",10)))
boxplot(ozone~label,notch=T,xlab="Garden",
  ylab="Ozone pphm",col="lightblue")
# the notches of two plots don't overlap, so the medians are significantly different at the 5% level

s2A<-var(gardenA)
s2B<-var(gardenB)
s2A/s2B
# to calculate whether the two variances are significantly different
# they are identical (s2A/s2B=1)

# the value of the test statistic for Student's t test is the difference divided by the standard error of the difference:

(mean(gardenA)-mean(gardenB))/sqrt(s2A/10+s2B/10)
# we divide the sample size and not the degree of freedom
# -3.872983
# with t tests you can ignore the minus sign
# the two means are significantly different

# we now use a two-tailed test

2*pt(-3.872983,18)
# p<0.0015

# there is also a built-in function to do all the work for you

t.test(gardenA,gardenB)
# the result is exactly the same as we obtained long-hand

# WILCOXON RANK-SUM TEST

# a non-parametric alternative to Student's t test
# both sample are put into a single array with their sample names clearly attached
# the aggregate is sorted, a rank is assigned to each value 
# the ranks are then added up for each of the two samples and significance is assessed on size of the smaller sum of ranks

ozone<-c(gardenA,gardenB)
ozone

label<-c(rep("A",10),rep("B",10))
label
# this makes a list of the sample names A and B

combined.ranks<-rank(ozone)
combined.ranks
# 'rank' makes a vector containing the ranks smallest to largest within the combined vector

# the ties have been dealt with by averaging the appropriate ranks

tapply(combined.ranks,label,sum)
# this calculates the sum of the ranks for each garden
# we compare the smallest of the two values 
# reject the null hypothesis after observing that the 5% value in tables is 78 and our value is smaller than this

wilcox.test(gardenA,gardenB)
# the p value is much less than 0.05, so we reject the null hypothesis
# conclude that the mean concentrations of ozone in gardens A and B are significantly different

# the non-parametric test is much more approrpiate than the t-test when the errors are not normal
# if a difference is significant under a Wilcoxon test it would be even more significant under a t-test

# TEST ON PAIRED SAMPLES

streams<-read.csv("c:\\MSc\\Statistics\\Data\\streams.csv")
attach(streams)
names(streams)

# the two samples were taken from the same river

t.test(down,up)
# by ignoring the fact that they are paired we get no change on p

t.test(down,up,paired=T)
# now the difference between the means is significanttly high

d<-up-down
t.test(d)
# this is the same paired test done as a one-sample t-test

# THE BINOMIAL TEST

binom.test(1,9)
# in this case the p value is the exact probability of obtaining the  observed result

# BINOMIAL TESTS TO COMPARE TWO PROPORTIONS

prop.test(c(4,196),c(40,3270))
# p=0.4696 so there is no evidence in favour
# this would occour more than 45% of the time by chance

# CHI-SQUARED CONTINGENCY TABLES

# a great deal of statistical information comes in the form of counts
# in statistics the contingencies are all the events that could possibly happen
# in order to make sense of the counts or frequencies we nee a model
# the simplest is that two categorical variables are independent
# if and only the two variables are independent then the probability of having one trait form one and one trait from the other is the product of two probabilities
# we then calculate the expected frequency: the probability multiplied by the total sample

#           Blue eyes Brown eyes
# Fair hair   38          11
# Dark hair   14          51 

# we work out the four expected frequencies in this case
# we need to figure out whether the expected frequencies are significantly different from the observed frequencies
# we assess the significance of the differences between using a chi-squared test
# we add the four components of chi-squared to get the test statistic
# we compare this value to the test statistic with the relevant critical value
# for this we need the number of degrees of freedom and the degree of certainty with which to work 

qchisq(0.95,1)
# this is the critical value
# if the calculated value of the test statistic is greater than the critical value we reject the null hypothesis

count<-matrix(c(38,14,11,51),nrow=2)
count
#qchisq.test(count) this uses Yate's correction as default
qchisq.test(count,correct=F)
# this does the whole procedure

# FISHER'S EXACT TEST

# used for the analysis of contingency tables based on small samples in which one or more of the expected frequencies are less than 5

#              Tree A Tree B Row totals
# With ants      6      2      8
# Without ants   4      8     12
# Column totals 10     10     20

factorial(8)*factorial(12)*factorial(10)*factorial(10)/ 
 (factorial(6)*factorial(2)*factorial(4)*factorial(8)*factorial(20))

# to compute the probability of outcomes that are more extreme than this:
# e.g. only one ant colony was found on tree B

factorial(8)*factorial(12)*factorial(10)*factorial(10)/
 (factorial(7)*factorial(3)*factorial(1)*factorial(9)*factorial(20))
 
# if not ant colonies were found on tree B:

factorial(8)*factorial(12)*factorial(10)*factorial(10)/
 (factorial(8)*factorial(2)*factorial(0)*factorial(10)*factorial(20))

# we need to add these probabilities together:

0.07501786 + 0.009526078 + 0.000352279

# we need to allow for extreme counts in the opposite direction by doubling this probability

2*(0.07501786+0.009526078+0.000352279)

x <- as.matrix(c(6,4,2,8))
dim(x) <- c(2,2)
x
fisher.test(x)
# this does all the computation above

# alternatively the function may be provided with two vectors containing factor levels instead of a matrix of counts:

table <- read.csv("c:\\MSc\\Statistics\\Data\\fisher.csv")
attach(table)
head(table)
fisher.test(tree,nests)

# CORRELATION AND COVARIANCE

# with two continuous variables x and y: are their values correlated with each other?
# we calculate the value of covariance of x and y
# this is the expectation of the vector product xXy: the expectation of the product minus the product of the two expectations

data<-read.csv("c\\MSc\\Statistics\\Data\\twosample.csv")
attach(data)
plot(x,y,pch=21,col="blue",bg="orange")

# we now need the variance of x and y:

var(x)
var(y)

# and the covariance

var(x,y)

# and the correlation coefficient

var(x,y)/sqrt(var(x)*var(y))

cor(x,y)
# simplified

# CORRELATION AND THE VARIANCE OF DIFFERENCES BETWEEN VARIABLES

paired <- read.csv("c:\\MSc\\Statistics\\Data\\water.table.csv ")
attach(paired)
names(paired)

cor(Summer, Winter)
# this shows a strong positive correlation (0.88201)

cor.test(Summer, Winter)
# this determines the significance of correlation
# it is highly significant (p=0.00165)

varS <- var(Summer)
varW <- var(Winter)
varD <- var(Summer-Winter)
# this determines the relationship between the correlation coefficient and the three variables including the variance of the differences

(varS+varW-varD)/(2*sqrt(varS)*sqrt(varW))
# we get a correlation coefficient of 0.88201

varD
varS + varW
# to see whether the variance of the difference is equal to the sum of the component variances
# 0.07821: it is not
# they would be equal only if the two samples were independent
# we know that the two variables are positively correlated, so the variance of difference should be less than the sum of variances

varS + varW - 2 * 0.8820102 * sqrt(varS) * sqrt(varW)
# 0.01015 
# now it makes sense

# SCALE-DEPENDENT CORRELATIONS

data <- read.csv("c:\\MSc\\Statistics\\Data\\productivity.csv")
attach(data)
names(data)
plot(productivity,mammals,pch=16,col="blue")
# it is a clear positive correlation

cor.test(productivity,mammals,method="spearman")
# the correlation is highly significant

plot(productivity,mammals,pch=16,col=as.numeric(region))
# this allows to look at different regions separately using a different colour
# increasing productivity is associated with reduced mammal species richness within each region
# need to be careful when looking at correlations across different scales
# warning: things that are positively correlated over short time scales may turn out to be negatively correlated in the long term and vice versa