pdf("g11.pdf");options(width=64)
#setwd("C:\\Users\\kolassa\\Class555")
setwd("~/Taught1/960-555/Data")
<pre>
#*************************************************************/
# Data from http://lib.stat.cmu.edu/datasets/Arsenic         */
# reformatted into ASCII on the course home page.  Data re-  */
# flect arsenic levels in toenail clippings; covariates in-  */
# clude age, sex (1=M), categorical measures of quantities   */
# used for drinking and cooking, arsenic in the water, and   */
# arsenic in the nails.  To make arsenic.dat from Arsenic, do*/
#antiword Arsenic|awk '((NR>39)&&(NR<61)){print}'>arsenic.dat*/
#*************************************************************/
arsenic<-read.table("arsenic.dat")
names(arsenic)<-c("age","sex","drink","cook","water",
   "nails")
arsenic$status<-(arsenic$water>=.0001)+0
arsenic$water[arsenic$water<.0001]<-.0001
arsenic$lnail<-log(arsenic$nails)
arsenic$lwater<-log(arsenic$water)
#<a name=block1> Block 1</a>
naildensity<-density(arsenic$nails)
plot(density(arsenic$nails,bw=naildensity$bw/10),xlab="Arsenic",
   type="l", main="Density estimate with inopportune bandwidths")
lines(density(arsenic$nails,bw=naildensity$bw*5),lty=2)
legend(1,2,legend=paste("Band width default",c("/10","*5")),
   lty=c(1,2))
#<a name=block2> Block 2</a>
library(NonparametricHeuristic)
t1<-fun.testreg(dists=c("rnorm","rcauchy","rlaplace","runif","rexp"))
t2<-fun.testreg(dists=c("rnorm","rcauchy","rlaplace","runif","rexp"),npergp=50)
#<a name=block3> Block 3</a>
library(quantreg)#For rq
rq(nails~water,data=arsenic)
#<a name=block4> Block 4</a>
rqo<-summary(rq(nails~water,data=arsenic))$coef
regp<-rqo[2,2]+diff(rqo[2,2:3])*(0:100)/100
m<-apply(abs(outer(arsenic$nails,rep(1,length(regp)))
   - outer(arsenic$water,regp)-rqo[1,1]),2,sum)
plot(regp,m,type="l",xlab="Regression Parameter", 
   ylab="Sum of absolute values of residuals",
   main="L1 fit for Arsenic Large Scale", 
   sub="Intercept fixed at optimum")
#<a name=block5> Block 5</a>
library(deming)#for theilsen
theilsen(arsenic$nails~arsenic$water,data=arsenic,conf=.9)$coef
#<a name=block6> Block 6</a>
library(MultNonParam);theil(arsenic$water,arsenic$nails)
#**************************************************************/
#* Data on blood presure change after captopril from          */
#* http://www.stat.ucla.edu/~rgould/datasets/bloodpressure.dat*/
#* Variables are                                              */
#* before systolic bp, after systolic bp, systolic bp change  */
#* before diastolic bp, after diastolic bp, diast. bp change  */
#* Example E from Cox and Snell, Applied Statistics: Princi-  */
#* ples and examples.  First line has column headings.  Fields*/
#* are separated by tabs.  SAS will give an error with this.  */
#**************************************************************/
bp<-as.data.frame(scan('bloodpressure.dat',skip=1,
  what=list(spb=0,spa=0,spd=0,dpb=0,dpa=0,dpd=0)))
#<a name=block7> Block 7</a>
#Expect warnings about nonunique estimators.  This is OK.
rqout<-rq(spa~spb,data=bp); summary(rqout)
#<a name=block8> Block 8</a>
rqoutt<-rq(spa~spb,tau=.2,data=bp); summary(rqoutt)
#<a name=block9> Block 9</a>
plot(bp$spb,bp$spa,main="Systolic Blood Pressure",
   xlab="Before Treatment",ylab="After Treatment")
abline(rqout); abline(rqoutt,lty=2)
legend(150,200,legend=c("Median Regression",
   ".2 Quantile Regresson"),lty=1:2)
#<a name=block10> Block 10</a>
library(MASS)#For rlm
rlmout<-rlm(spa~spb,data=bp)
library(quantreg)#For rq
rqout<-rq(spa~spb,data=bp)
plot(bp$spb,bp$spa, main="Blood Pressure After and Before Treatment",
   xlab="Systolic Blood Pressure Before Treatment (mm HG)",
   ylab="Systolic Blood Pressure After Treatment (mm HG)",
   sub="Linear Fits to Original Data")
abline(rlmout,lty=2)
abline(rqout,lty=3)
abline(lm(spa~spb,data=bp),lty=4)
legend(min(bp$spb),max(bp$spa),
 legend=c("Huber","Median","OLS"),lty=2:4)
#<a name=block11> Block 11</a>
bp$spa[1]<-120
rqout<-rq(spa~spb,data=bp); rlmout<-rlm(spa~spb,data=bp)
plot(bp$spb,bp$spa,
   xlab="Systolic Blood Pressure Before Treatment (mm HG)",
   ylab="Systolic Blood Pressure After Treatment (mm HG)",
   sub="Linear Fits to Perturbed Data")
abline(rlmout,lty=2);abline(rqout,lty=3);
abline(lm(spa~spb,data=bp),lty=4)
legend(min(bp$spb),max(bp$spa),
 legend=c("Huber","Median","OLS"),lty=2:4)
#*************************************************************/
# Data from http://lib.stat.cmu.edu/datasets/Arsenic         */
# reformatted into ASCII on the course home page.  Data re-  */
# flect arsenic levels in toenail clippings; covariates in-  */
# clude age, sex (1=M), categorical measures of quantities   */
# used for drinking and cooking, arsenic in the water, and   */
# arsenic in the nails.  To make arsenic.dat from Arsenic, do*/
#antiword Arsenic|awk '((NR>39)&&(NR<61)){print}'>arsenic.dat*/
#*************************************************************/
arsenic<-read.table("arsenic.dat")
names(arsenic)<-c("age","sex","drink","cook","water",
   "nails")
arsenic$status<-(arsenic$water>=.0001)+0
arsenic$water[arsenic$water<.0001]<-.0001
arsenic$lnail<-log(arsenic$nails)
arsenic$lwater<-log(arsenic$water)
#<a name=block12> Block 12</a>
plot(arsenic$water,arsenic$nails,
  main="Nail and Water Arsenic Levels", xlab="Water",ylab="Nails",
  sub="Smoother fits.  Bandwidth chosen by inspection.")
#<a name=block13> Block 13</a>
lines(ksmooth(arsenic$water,arsenic$nails,"normal",bandwidth=0.10),lty=1)
lines(ksmooth(arsenic$water,arsenic$nails,"box",bandwidth=0.10),lty=2)
legend(0.0,2.0,lty=1:3, legend=c("Smoother, Normal Kernel",
   "Smoother, Box Kernel", "Local Polynomial"))
#<a name=block14> Block 14</a>
library(KernSmooth)
lines(locpoly(arsenic$water,arsenic$nails,bandwidth=.05,degree=1),lty=3)
#<a name=block15> Block 15</a>
plot(arsenic$water,arsenic$nails,main="Nail and Water Arsenic Levels",
  xlab="Water",ylab="Nails", sub="Loess fits")
#<a name=block16> Block 16</a>
x<-min(arsenic$water)+(0:50)*diff(range(arsenic$water))/50
lines(x,y=predict(loess(nails~water,data=arsenic),x))
lines(x,y=predict(loess(nails~water,data=arsenic,span=1),x),lty=2)
#<a name=block17> Block 17</a>
lines(x,y=predict(loess(nails~water,data=arsenic,span=10),x),lty=3)
legend(0.0,2.0,legend=c("Default span .75",
  "Span 1","Span 10"),lty=1:3)                   
#<a name=block18> Block 18</a>
plot(isoreg(arsenic$water,arsenic$nails), xlab="Water As",
   ylab="Nail As", main="Nail Arsenic and Water Arsenic")
#<a name=block19> Block 19</a>
plot(arsenic$water,arsenic$nails, main="Nail Arsenic and Water Arsenic")
hgrid<-min(arsenic$water)+(0:50)*diff(range(arsenic$water))/50
lines(predict(spl<-smooth.spline(arsenic$water,arsenic$nails),hgrid))
lines(predict(smooth.spline(arsenic$water,arsenic$nails,spar=.1),hgrid),
  lty=3)
lines(predict(smooth.spline(arsenic$water,arsenic$nails,spar=.5),hgrid),
  lty=2)
legend(1250,1110,lty=1:3,col=1:3,
   legend=paste("Smoothing",c(round(spl$spar,3),.5,.1),
   c("Default","","")))
</pre>