R version 4.2.2 Patched (2022-11-10 r83330) -- "Innocent and Trusting"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> pdf("g11.pdf");options(width=64)
> #setwd("C:\\Users\\kolassa\\Class555")
> setwd("~/Taught1/960-555/Data")
> #*************************************************************/
> # 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)
> 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))
> library(NonparametricHeuristic)
> t1<-fun.testreg(dists=c("rnorm","rcauchy","rlaplace","runif","rexp"))
> t2<-fun.testreg(dists=c("rnorm","rcauchy","rlaplace","runif","rexp"),npergp=50)
> library(quantreg)#For rq
> rq(nails~water,data=arsenic)
Call:
rq(formula = nails ~ water, data = arsenic)

Coefficients:
(Intercept)       water 
  0.1141978  15.6043956 

Degrees of freedom: 21 total; 19 residual
> 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")
> library(deming)#for theilsen
> theilsen(arsenic$nails~arsenic$water,data=arsenic,conf=.9)$coef
  (Intercept) arsenic$water 
    0.1167158    14.2156504 
> library(MultNonParam);theil(arsenic$water,arsenic$nails)
$est
[1] NA

$ci
[1]  8.170732 16.535948

> #**************************************************************/
> #* 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)))
> #Expect warnings about nonunique estimators.  This is OK.
> rqout<-rq(spa~spb,data=bp); summary(rqout)

Call: rq(formula = spa ~ spb, data = bp)

tau: [1] 0.5

Coefficients:
            coefficients lower bd  upper bd 
(Intercept) -11.41026    -34.76720  29.00312
spb           0.94872      0.69282   1.11167
> rqoutt<-rq(spa~spb,tau=.2,data=bp); summary(rqoutt)

Call: rq(formula = spa ~ spb, tau = 0.2, data = bp)

tau: [1] 0.2

Coefficients:
            coefficients lower bd  upper bd 
(Intercept)  17.25000    -68.48262  23.34294
spb           0.75000      0.40129   1.20640
> 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)
> 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)
> 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)
> plot(arsenic$water,arsenic$nails,
+   main="Nail and Water Arsenic Levels", xlab="Water",ylab="Nails",
+   sub="Smoother fits.  Bandwidth chosen by inspection.")
> 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"))
> library(KernSmooth)
> lines(locpoly(arsenic$water,arsenic$nails,bandwidth=.05,degree=1),lty=3)
> plot(arsenic$water,arsenic$nails,main="Nail and Water Arsenic Levels",
+   xlab="Water",ylab="Nails", sub="Loess fits")
> 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)
> 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)                   
> plot(isoreg(arsenic$water,arsenic$nails), xlab="Water As",
+    ylab="Nail As", main="Nail Arsenic and Water Arsenic")
> 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","","")))
>