##### Lesson 99: Two-Sample Hypothesis Tests in R #####

# Setting the working directory 
setwd("path to your folder")
setwd(getwd())


# Read the data file 
nyc_trash_data = read.csv("DSNY_Monthly_Tonnage_Data.csv",header=T)


# Extract February refuse data for Manhattan's community district 9
feb_index = which((nyc_trash_data$BOROUGH=="Manhattan") & (nyc_trash_data$COMMUNITYDISTRICT==9) & (nyc_trash_data$MONTH==2))
feb_refuse_data = nyc_trash_data$REFUSETONSCOLLECTED[feb_index]


# Extract February refuse data for Manhattan's community district 9
aug_index = which((nyc_trash_data$BOROUGH=="Manhattan") & (nyc_trash_data$COMMUNITYDISTRICT==9) & (nyc_trash_data$MONTH==8))
aug_refuse_data = nyc_trash_data$REFUSETONSCOLLECTED[aug_index]


# Visualizing the distributions of the two samples
boxplot(cbind(feb_refuse_data,aug_refuse_data),horizontal=T, main="Refuse from Manhattan's Community District 9")
text(1500,1,"February Tonnage",font=2)
text(1500,2,"August Tonnage",font=2)

p_threshold = 2500 # tons of refuse
abline(v=p_threshold,lty=2,col="maroon")


# Preliminaries 
  # sample sizes
  n1 = length(feb_refuse_data)
  n2 = length(aug_refuse_data)
  
  # sample means 
  x1bar = mean(feb_refuse_data)
  x2bar = mean(aug_refuse_data)
  
  # sample variances 
  x1var = var(feb_refuse_data)
  x2var = var(aug_refuse_data)
  
  # sample proportions
  p1 = length(which(feb_refuse_data<p_threshold))/n1
  p2 = length(which(aug_refuse_data<p_threshold))/n2
  
  
# Hypothesis Test on the Difference in Means 

  #1. t-Test
  pooled_var = ((n1-1)/(n1+n2-2))*x1var + ((n2-1)/(n1+n2-2))*x2var
  
  t0 = (x1bar-x2bar)/sqrt(pooled_var*((1/n1)+(1/n2)))
  
  df = n1+n2-2
  
  pval = pt(t0,df=df)
  
  print(pooled_var)
  print(df)
  print(t0)
  print(pval)
  
  t.test(feb_refuse_data,aug_refuse_data,alternative="two.sided",var.equal = TRUE)

  #2. Welch's t-Test
  f = (((x1var/n1)+(x2var/n2))^2)/(((x1var/n1)^2/(n1-1))+((x2var/n2)^2/(n2-1)))
  
  t0 = (x1bar-x2bar)/sqrt((x1var/n1)+(x2var/n2))
  
  pval = pt(t0,df=f)
  
  print(f)
  print(t0)
  print(pval)
  
  t.test(feb_refuse_data,aug_refuse_data,alternative="two.sided",var.equal = FALSE)

  #3. Wilcoxon's Rank-sum Test
  wilcox.test(feb_refuse_data,aug_refuse_data,alternative = "two.sided")
  
  #4. Bootstrap
  N = 10000
  null_mean = matrix(0,nrow=N,ncol=1)
  null_mean_ratio = matrix(0,nrow=N,ncol=1)

  for(i in 1:N)
  {
    xboot = sample(feb_refuse_data,replace=T)
    yboot = sample(aug_refuse_data,replace=T)
    
    null_mean_ratio[i] = mean(xboot)/mean(yboot)

    if(mean(xboot)>mean(yboot)){null_mean[i]=1} 
  }
  
  pvalue_mean = sum(null_mean)/N

  hist(null_mean_ratio,font=2,main="Null Distribution Assuming H0 is True",xlab="Xbar/Ybar",font.lab=2)
  abline(v=1,lwd=2,lty=2)
  text(0.95,1000,paste("p-value=",pvalue_mean),col="red")
  
  
# Hypothesis Test on the Equality of Variances
  
  #1. F-Test
  f0 = x1var/x2var
  
  df_numerator = n1-1
  
  df_denominator = n2-1
  
  pval = 1-pf(f0,df1=df_numerator,df2=df_denominator)
  
  print(f0)
  print(df_numerator)
  print(df_denominator)
  print(pval)
  
  var.test(feb_refuse_data,aug_refuse_data,alternative = "two.sided")

  #2. Bootstrap
  N = 10000
  null_var = matrix(0,nrow=N,ncol=1)
  null_var_ratio = matrix(0,nrow=N,ncol=1)
  
  for(i in 1:N)
  {
    xboot = sample(feb_refuse_data,replace=T)
    yboot = sample(aug_refuse_data,replace=T)
    
    null_var_ratio[i] = var(xboot)/var(yboot)
    
    if(var(xboot)>var(yboot)){null_var[i]=1} 
  }
  
  pvalue_var = sum(null_var)/N
  
  hist(null_var_ratio,font=2,main="Null Distribution Assuming H0 is True",xlab="XVar/YVar",font.lab=2)
  abline(v=1,lwd=2,lty=2)
  text(2,500,paste("p-value=",pvalue_var),col="red")
  
  
# Hypothesis Test on the Difference in Proportions
  
  #1. Fisher's Exact Test  
  x1 = length(which(feb_refuse_data<p_threshold))
  p1 = x1/n1
  
  x2 = length(which(aug_refuse_data<p_threshold))
  p2 = x2/n2
  
  N = n1+n2
  t = x1+x2
  k = seq(from=0,to=n1,by=1)
  
  p = k
  for(i in 1:length(k))
  {
    p[i] = (choose(t,k[i])*choose((N-t),(n1-k[i])))/choose(N,n1)
  }
  
  plot(k,p,type="h",xlab="Number of successes in X1",ylab="P(X=k)",font=2,font.lab=2)
  points(k,p,type="o",lty=2,col="grey50")
  points(k[13:length(k)],p[13:length(k)],type="o",col="red",lwd=2)
  points(k[13:length(k)],p[13:length(k)],type="h",col="red",lwd=2)
  #abline(h=0,col='red',lwd=2)
  
  pvalue = sum(p[13:length(k)])
  print(pvalue)
  
  fisher_data = cbind(c(x1,x2),c((n1-x1),(n2-x2)))
  fisher.test(fisher_data,alternative="greater")  
    
  #2. Z-approximation
  p = (x1+x2)/(n1+n2)
  z = (p1-p2)/sqrt(p*(1-p)*((1/n1)+(1/n2)))
  pval = 1-pnorm(z)
  
  #3. Bootstrap
  N = 10000
  null_prop = matrix(0,nrow=N,ncol=1)
  null_prop_ratio = matrix(0,nrow=N,ncol=1)
  
  for(i in 1:N)
  {
    xboot = sample(feb_refuse_data,replace=T)
    yboot = sample(aug_refuse_data,replace=T)
    
    p1boot = length(which(xboot<p_threshold))/n1
    p2boot = length(which(yboot<p_threshold))/n2
    
    null_prop_ratio[i] = p1boot/p2boot
    
    if(p1boot>p2boot){null_prop[i]=1} 
  }
  
  pvalue_prop = sum(null_prop)/N
  
  hist(null_prop_ratio,font=2,main="Null Distribution Assuming H0 is True",xlab="P1/P2",font.lab=2)
  abline(v=1,lwd=2,lty=2)
  text(2,250,paste("p-value=",pvalue_prop),col="red")