# teset one mean 
t.test(var1,mu=mu)

# test 2 means
t.test(var1,var2) # both numeric variables

# graphs
# histogram
hist(var1)
# boxplot
boxplot(var1)
# barplot (categorical data only)
# you may have to arrange your data like mine here with table()
# people making a complete stop at a 4-way stop sign on D and Hayes St. in Moscow, Idaho.
# 45 did and 100 did not
stop=c(rep(1,each=45),rep(0,each=100))
stop

# make a table to be able to graph it

tstop=table(stop)
tstop
barplot(tstop,main='Proportion of Complete \n Stops at Hayes and D St.',col=grey.colors(2))
# make values something other than 0, 1
stop2=c(rep("Yes",each=45),rep("No",each=100))
stop2
tstop2=table(stop2)
tstop2

# graph with barplot(). You will need to do a title
barplot(tstop2,main='Proportion of Complete \n Stops at Hayes and D St.',col=cm.colors(2))



# 1 sample proportion test
# alt: can be 'g' for Ha>, 'l' for Ha<, 'ne' for Ha not equal
test.1prop=function(x,n,p0,alt,alpha) { 
  phat=x/n; p0=p0; q0=1-p0
  sep0=sqrt(p0*q0/n); zcalc=(phat-p0)/sep0
  alt=alt; alpha=alpha; sig=sprintf("%.0f%%",alpha*100)
  if(alt=='g'){
    pvalue=pnorm(zcalc,lower.tail=F)
  }    else if(alt=='l'){  
    pvalue=pnorm(zcalc)
  }    else if(alt=='ne'){
    pvalue=2*pnorm(-abs(zcalc))
  }
  altname=c(if(alt=='g'){
    altname=c("greater than")
  }    else if(alt=='l'){  
    altname=c("less than")
  }    else if(alt=='ne'){
    altname=c("not equal to")
  })
  cat(" ","\n","Hypothesis Test for 1-sample proportion","\n",
      "Number of successes =",x,"\n","Alternative hypothesis =",altname,"\n","Sample size =",n,
      "\n","Hypothesized proportion =",p0,"\n", "SE =",sep0,"\n","Significance Level =",
      sig,"\n","Zcalc =",zcalc,"\n","pvalue =",pvalue,"\n","","\n")
  results=list(phat=phat,n=n,x=x,p0=p0,sep0=sep0,alt=alt,zcalc=zcalc,pvalue=pvalue)
}  
# H0: p=.11 Ha: p not=.11
x=15; n=182; p0=.11
test.1prop(x=x,n=n,p0=p0,alt='ne',alpha=0.05)

# 1 sample proportion CI
# you must run the whole function from ci.1prop through the last } for this to work
ci.1prop=function(x,n,alpha){ 
  phat=x/n; qhat=1-phat; alpha=alpha; zstar=qnorm(1-alpha/2)
  cl=(1-alpha); cl.perc=sprintf("%.0f%%",cl*100)
  se=sqrt(phat*qhat/n); moe=zstar*se
  lower=phat-moe; upper=phat+moe
  lowerpct=sprintf("%.2f%%",lower*100); upperpct=sprintf("%.2f%%",upper*100)
  cat(" ","\n","CI for 1-sample proportion","\n","Confidence Level =",cl.perc,
      "\n","Critical Value of Z* =",zstar,"\n","SE =",se,"\n","Point Estimate of p =",phat,"\n",
      "Margin of Error =",moe,"\n","Lower Bound =",lower,"Upper Bound =",upper,"\n","Lower Bound (%) =",
      lowerpct,"Upper Bound (%) =",upperpct,"\n","","\n")
  results=list(zstar=zstar,se=se,phat=phat,moe=moe,lower=lower,upper=upper,
               lowerpct=lowerpct,upperpct=upperpct)
}   
# CI with x=15, n=182
ci.1prop(x=15,n=182,0.02)

# 2 sample proportion test 
# you must run the whole function from test.2prop through the last } for this to work
test.2prop=function(x1,x2,n1,n2,alt,alpha) { 
  phat1=x1/n1; qhat1=1-phat1; phat2=x2/n2; qhat2=1-phat2; pt.est=phat1-phat2
  pooled.phat=(x1+x2)/(n1+n2); pooled.qhat=1-pooled.phat
  se=sqrt((pooled.phat*pooled.qhat)*(1/n1+1/n2)); zcalc=pt.est/se
  alpha=alpha; alt=alt; sig=sprintf("%.0f%%",alpha*100)
  if(alt=='g'){
    pvalue=pnorm(zcalc,lower.tail=F)
  }    else if(alt=='l'){  
    pvalue=pnorm(zcalc)
  }    else if(alt=='ne'){
    pvalue=2*pnorm(-abs(zcalc))
  }
  cat(" ","\n","Hypothesis Test for 2-sample proportion","\n","p1 =",phat1,"p2 =",phat2,"\n","Point Estimate of p1-p2 =",pt.est,"\n",
      "Pooled p =",pooled.phat,"SE =",se,"\n","Alternative hypothesis =",alt,"Significance Level =",sig,"\n",
      "\n","Zcalc =",zcalc,"\n","pvalue =",pvalue,"\n")
  results=list(phat1=phat1,n1=n1,phat2=phat2,n2=n2,pt.est=pt.est,pooled.phat=pooled.phat,
               se=se,alt=alt,zcalc=zcalc,pvalue=pvalue)
}  
# a
# H0: p1=p2 Ha: p1>p2
test.2prop(x1=15,x2=92,n1=182,n2=2300,alt='g',0.05)

# 2 sample proportion CI
# you must run the whole function from ci.2prop through the last } for this to work
ci.2prop=function(x1,x2,n1,n2,alpha){ 
  phat1=x1/n1; qhat1=1-phat1; phat2=x2/n2; qhat2=1-phat2 
  pt.est=phat1-phat2; alpha=alpha; zstar=qnorm(1-alpha/2)
  se=sqrt((phat1*qhat1/n1)+(phat2*qhat2/n2)); moe=zstar*se
  lower=pt.est-moe; upper=pt.est+moe
  lowerpct=sprintf("%.2f%%",lower*100); upperpct=sprintf("%.2f%%",upper*100)
  print("Results of the Confidence Interval of the Difference of 2 Proportions")
  cat(" Critical Value of Z* =",zstar,"\n","SE =",se,"\n","Point Estimate of p1-p2 =",
      pt.est,"\n","Margin of Error =",moe,"\n","Lower Bound =",lower,
      "Upper Bound =",upper,"\n","Lower Bound (%) =",lowerpct,
      "Upper Bound (%) =",upperpct,"\n")
  results=list(zstar=zstar,se=se,pt.est=pt.est,moe=moe,lower=lower,upper=upper,
               lowerpct=lowerpct,upperpct=upperpct)
}   
# CI with x1=15,x2=92,n1=182,n2=2300
ci.2prop(x1=15,x2=92,n1=182,n2=2300,0.02)