### R code from vignette source 'helpSASR3.SASnw'

###################################################
### code chunk number 1: helpSASR3.SASnw:11-12
###################################################
options(continue="  ")


###################################################
### code chunk number 2: helpSASR3.SASnw:51-54
###################################################
options(digits=3)
options(width=72)  # narrows output to stay in the grey box
ds = read.csv("http://www.amherst.edu/~nhorton/sasr2/datasets/help.csv")


###################################################
### code chunk number 3: helpSASR3.SASnw:90-91
###################################################
with(ds, mean(cesd))


###################################################
### code chunk number 4: helpSASR3.SASnw:95-102
###################################################
library(mosaic)
mean(~ cesd, data=ds)
with(ds, median(cesd))
with(ds, range(cesd))
sd(~ cesd, data=ds)
var(~ cesd, data=ds)
favstats(~ cesd, data=ds)


###################################################
### code chunk number 5: helpSASR3.SASnw:107-110
###################################################
library(moments)
with(ds, skewness(cesd))
with(ds, kurtosis(cesd))


###################################################
### code chunk number 6: helpSASR3.SASnw:138-139
###################################################
with(ds, quantile(cesd, seq(from=0, to=1, length=11)))


###################################################
### code chunk number 7: bplotr (eval = FALSE)
###################################################
## with(ds, hist(cesd, main="", freq=FALSE))
## with(ds, lines(density(cesd), main="CESD", lty=2, lwd=2))
## xvals = with(ds, seq(from=min(cesd), to=max(cesd), length=100))
## with(ds, lines(xvals, dnorm(xvals, mean(cesd), sd(cesd)), lwd=2))


###################################################
### code chunk number 8: bplot
###################################################
with(ds, hist(cesd, main="", freq=FALSE))
with(ds, lines(density(cesd), main="CESD", lty=2, lwd=2))
xvals = with(ds, seq(from=min(cesd), to=max(cesd), length=100))
with(ds, lines(xvals, dnorm(xvals, mean(cesd), sd(cesd)), lwd=2))


###################################################
### code chunk number 9: cormat1
###################################################
cormat = cor(with(ds, cbind(cesd, mcs, pcs)))
cormat


###################################################
### code chunk number 10: cormat2
###################################################
cormat[c(2, 3), 1]


###################################################
### code chunk number 11: rugplotr (eval = FALSE)
###################################################
## with(ds, plot(cesd[female==1], mcs[female==1], xlab="CESD", ylab="MCS",
##    type="n", bty="n"))
## with(ds, text(cesd[female==1&substance=="alcohol"],
##    mcs[female==1&substance=="alcohol"],"A"))
## with(ds, text(cesd[female==1&substance=="cocaine"],
##    mcs[female==1&substance=="cocaine"],"C"))
## with(ds, text(cesd[female==1&substance=="heroin"],
##    mcs[female==1&substance=="heroin"],"H"))
## with(ds, rug(jitter(mcs[female==1]), side=2))
## with(ds, rug(jitter(cesd[female==1]), side=3))


###################################################
### code chunk number 12: splot
###################################################
with(ds, plot(cesd[female==1], mcs[female==1], xlab="CESD", ylab="MCS",
   type="n", bty="n"))
with(ds, text(cesd[female==1&substance=="alcohol"],
   mcs[female==1&substance=="alcohol"],"A"))
with(ds, text(cesd[female==1&substance=="cocaine"],
   mcs[female==1&substance=="cocaine"],"C"))
with(ds, text(cesd[female==1&substance=="heroin"],
   mcs[female==1&substance=="heroin"],"H"))
with(ds, rug(jitter(mcs[female==1]), side=2))
with(ds, rug(jitter(cesd[female==1]), side=3))


###################################################
### code chunk number 13: helpSASR3.SASnw:405-408
###################################################
tally(~ homeless + female, data=ds)
tally(~ homeless + female, format="percent", data=ds)
tally(~ homeless | female, data=ds)


###################################################
### code chunk number 14: helpSASR3.SASnw:424-429
###################################################
or = with(ds, (sum(homeless==0 & female==0)*
       sum(homeless==1 & female==1))/
      (sum(homeless==0 & female==1)*
       sum(homeless==1 & female==0)))
or


###################################################
### code chunk number 15: helpSASR3.SASnw:436-440
###################################################
library(epitools)
oddsobject = with(ds, oddsratio.wald(homeless, female))
oddsobject$measure
oddsobject$p.value


###################################################
### code chunk number 16: helpSASR3.SASnw:451-453
###################################################
chisqval = with(ds, chisq.test(homeless, female, correct=FALSE))
chisqval


###################################################
### code chunk number 17: helpSASR3.SASnw:461-462
###################################################
with(ds, fisher.test(homeless, female))


###################################################
### code chunk number 18: gridtable (eval = FALSE)
###################################################
## library(gridExtra)
## mytab = tally(~ racegrp + substance, data=ds)
## plot.new()
## grid.table(mytab)


###################################################
### code chunk number 19: helpSASR3.SASnw:555-557
###################################################
ttres = t.test(age ~ female, data=ds)
print(ttres)


###################################################
### code chunk number 20: helpSASR3.SASnw:589-592
###################################################
library(coin)
oneway_test(age ~ as.factor(female),
   distribution=approximate(B=9999), data=ds)


###################################################
### code chunk number 21: helpSASR3.SASnw:654-655
###################################################
with(ds, wilcox.test(age ~ as.factor(female), correct=FALSE))


###################################################
### code chunk number 22: kstest
###################################################
ksres = with(ds, ks.test(age[female==1], age[female==0]))
print(ksres)


###################################################
### code chunk number 23: helpSASR3.SASnw:787-796
###################################################
plotdens = function(x,y, mytitle, mylab) {
   densx = density(x)
   densy = density(y)
   plot(densx, main=mytitle, lwd=3, xlab=mylab, bty="l")
   lines(densy, lty=2, col=2, lwd=3)
   xvals = c(densx$x, rev(densy$x))
   yvals = c(densx$y, rev(densy$y))
   polygon(xvals, yvals, col="gray")
}


###################################################
### code chunk number 24: polyplotr (eval = FALSE)
###################################################
## mytitle = paste("Test of ages: D=", round(ksres$statistic, 3),
##    " p=", round(ksres$p.value, 2), sep="")
## with(ds, plotdens(age[female==1], age[female==0], mytitle=mytitle, 
##    mylab="age (in years)"))
## legend(50, .05, legend=c("Women", "Men"), col=1:2, lty=1:2, lwd=2)


###################################################
### code chunk number 25: kdeplot
###################################################
mytitle = paste("Test of ages: D=", round(ksres$statistic, 3),
   " p=", round(ksres$p.value, 2), sep="")
with(ds, plotdens(age[female==1], age[female==0], mytitle=mytitle, 
   mylab="age (in years)"))
legend(50, .05, legend=c("Women", "Men"), col=1:2, lty=1:2, lwd=2)


###################################################
### code chunk number 26: helpSASR3.SASnw:873-878
###################################################
library(survival)
survobj = survdiff(Surv(dayslink, linkstatus) ~ treat, 
   data=ds)
print(survobj)
names(survobj)