RDX2 X   g13.1 source function(device=""){ $ oldpal <- palette(gray(c(0,.6))) 4 if(device!="") hardcopy(width=2.75, height=2.75, * device=device) 6 oldpar <- par(mar = par()$mar + c(-0.5, 0, -2, 0))  on.exit(par(oldpar)) # require(DAAG); data(socsupport) B not.na <- apply(socsupport[,9:19],1,function(x)!any(is.na(x)))  not.na[36] <- FALSE F ss.pcp1 <- princomp(socsupport[not.na, 9:19], cor=TRUE)$scores[,1]  attach(socsupport) ? plot(BDI[not.na] ~ ss.pcp1, col=as.numeric(gender[not.na]), ( pch=as.numeric(gender[not.na]), 7 xlab ="First principal component", ylab="BDI") topleft <- par()$usr[c(1,4)] 4 legend(topleft[1], topleft[2], col=1:2, pch=1:2, * legend=levels(gender), cex=0.7)  detach(socsupport)  if(device!="")dev.off()  palette(oldpal)  invisible()  } device  srcref  srcfile encoding native.enc timestampAҳJ8@ class POSIXt POSIXct filename ~/r-book/ed2/figures/figs13.R wd /Users/johnm/r/ch13  srcfile  srcref $$  srcref **  srcref 66  srcref   srcref   srcref ##  srcref   BB  srcref     srcref   FF  srcref     srcref  77  srcref     srcref **  srcref   srcref   srcref   srcref   srcref { <- oldpal palette gray c?333333 if !=  hardcopy width@ height@  oldpar par mar + $ -?@ on.exit require DAAG data socsupport  not.na apply [! :@"@3? function x ! any is.na' function(x)!any(is.na(x)) $"@B   ss.pcp1$ princomp$!"%@"@3 cor  scores? attach! plot ~$ BDI"+ col as.numeric$ gender" pch4$5" xlab First principal component ylab BDI  topleft$ usr?@ legend$9?$9@3%?@6%?@; levels5 cex?ffffff detach!  dev.off invisible g13.2# .function(df=nsw74psid1, device="", lty=c(1,1), , colr=c("black","gray"), lwd=1:2){ 1 if(device!="")hardcopy(width=5.5, height=1.6, 6 device=device, pointsize=9) M oldpar <- par(mfrow=c(1,5), mar=c(4.1,2.6,2.1,0.1), oma=c(0,2.1,2.1,1.1), , mgp=c(2.25,.5,0), pty="s")  on.exit(par(oldpar))  attach(df) k <- 0 3 for(i in c("educ","age","re74","re75","re78")){  k <- k+1  y <- df[,i]  d0 <- density(y[trt==0])  d1 <- density(y[trt==1]) " ylim <- range(c(d0$y, d1$y)) ## if(k %in% c(1,4)) P## plot(d0$x, d0$y, ylim=ylim, main="", xlab=i, type="l", ylab="Density", '## lwd=lwd[1], col=colr[1]) ## else G plot(d0$x, d0$y, ylim=ylim, main="", xlab=i, type="l", ylab="", 2 lwd=lwd[1], col=colr[1], cex.lab=1.5) < lines(d1$x, d1$y, lty=lty[2], lwd=lwd[2], col=colr[2])  } ;# plot(0:1, 0:1, type="n", axes=FALSE, xlab="", ylab="") 2 mtext(side=2, line=0.1, "Density", outer=TRUE) = par(fig=c(0,1,0,1), mar=rep(0,4), oma=rep(0,4), new=TRUE) xpos <- mean(par()$usr[1:2])  ypos <- par()$usr[4]  par(xpd=TRUE) ? legend(x=xpos, y=ypos, lty=lty, col=colr, lwd=lwd, bty="n", = legend=c("control", "treated"), horiz=T, inset=-3,  xjust=0.75, yjust=0)  if(device!="")dev.off()  detach(df)  } df nsw74psid1  lty?? colr black gray lwd%?@ ,,,,  srcref 66  srcref ,,  srcref     srcref !!  srcref ""    srcref #0  srcref 2222  srcref 33==  srcref 44    srcref 55  srcref 66  srcref 79  srcref ::  srcref ;;  srcref  @? pointsize@"  mfrow?@@ffffff@@? oma@@?񙙙 mgp@? pty s/B  k for i educ age re74 re75 re78 #3#333  srcref $$  srcref %%  srcref &&  srcref ''  srcref ((""  srcref - .2 2  srcref //<<  srcref  LL?  y$BN  d0 density$O == trt  d1Q$ORS?  ylim rangePOTO0P'POUU main 7N type l8 F$F?3$E? cex.lab? linesT'TOD$D@F$F@3$E@ mtext side@ line? Density outer  fig?? rep@I`@ new   xpos mean$:%?@  ypos$:@ xpd ;'bOdDD3EFF bty n; control treated horiz T inset@ xjust? yjust ?>B g13.3 4function(device="", df=nsw74psid1, maxf=30, bw=TRUE) {  library(splines) @ if(device!="")hardcopy(width=2.25, height=2.25, device=device, 8 trellis=TRUE, pointsize=c(8,5)) , oldpar <- par(mar = c(4.1, 4.1, 1.1, 1.1), ' mgp = c(2.25, 0.5, 0) )  on.exit(par(oldpar))  common <- multilap(maxf=maxf) * offset <- mean(sapply(nsw74psid1[,8:10], 7 function(x)unique(sort(x))[2])) F if(bw)trellis.par.set(superpose.symbol=list(col="gray", pch=c(1,3)), C superpose.line=list(col=c("black", "gray"), 2 lwd=c(1,2)), lty=c(1,1)) & A.lm <- lm(log(re78+offset) ~ trt + 3 (ns(age,4) + educ + log(re74+offset) + 8 log(re75+offset)) + I(re74>0) + I(re75>0)+ : (black + hisp + marr + nodeg), subset=common,  data = df)  print(summary(A.lm)$coef) 6 trt <- factor(df$trt[common], labels=c("ctl","trt")) 6 print(xyplot(resid(A.lm) ~ fitted(A.lm), groups=trt, 4 auto.key=list(columns=2, lines=TRUE), % type=c("p","smooth")))  if(device!="")dev.off() invisible() } BC maxf@> bw  @@  srcref AA  srcref BC88  srcref DE''  srcref FF  srcref GG  srcref HI77  srcref JL22  srcref MQ  srcref RR  srcref SS66  srcref TV%%  srcref WW  srcref XX    srcref  library splines @@ trellis G@ @ @ffffff@ffffff?񙙙?񙙙J@?  common multilapmm  offsetc sapply$C%@ @$&'$ unique sort'@ function(x)unique(sort(x))[2]n trellis.par.set superpose.symbol list3 gray6?@ superpose.linez3 black grayF?@D??  A.lm lm1 log re78tS ( ns age@ educ~ re74t~ re75t I > black hisp marr nodeg subsetr B print summary| coef S factor$BSr labels ctl trt xyplot1 resid| fitted| groupsS auto.keyz columns@Z X p smooth ?@ g13.4: Dfunction(df=nsw74psid1, device="", maxf=30, maxf.score=30, bw=TRUE){ B if(device!="")hardcopy(width=4.25, height=2.35, trellis=TRUE, F device=device, pointsize=c(8,5), color=!bw) 2 oldpar <- par(mar=c(2.1,4.1,2.6,2.1), cex=0.7, , mgp=c(2.25,.5,0), pty="s")  on.exit(par(oldpar)) ! common <- multilap(maxf=maxf) , offset <- mean(sapply(nsw74psid1[,8:10], 9 function(x)unique(sort(x))[2])) @ disc.glm <- glm(trt ~ (ns(age,4) + educ + log(re74+offset) + E log(re75+offset)) + I(re74>0) + I(re75>0)+ A (black + hisp + marr + nodeg), subset=common, 3 data = df, family=binomial) Pscores <- predict(disc.glm) & nsw74psidB <- nsw74psid1[common, ]  trt <- nsw74psidB$trt  par(fig=c(0, 0.5, 0, 1)) ( if(bw)colr <- c("black","gray") else 2 colr <- trellis.par.get()$superpose.symbol$col  par(xpd=TRUE) I xchop <- overlapDensity(Pscores[trt==0], Pscores[trt==1], lty=c(1,1), < compare.numbers=FALSE, col=colr, > ratio=c(1/maxf.score, maxf.score)) 4 mtext(side=3, line=2.85, adj=-0.2, "A", cex=0.8)  par(xpd=FALSE)  par(fig=c(0.5, 1, 0, 1))  par(new=TRUE) 8 plot(0:1,0:1, xlab="", ylab="",axes=FALSE, type="n") 5 mtext(side=3, line=2.85, adj=-0.35, "B", cex=0.8) 6 overlap <- Pscores > xchop[1] & Pscores < xchop[2] @ print(with(subset(nsw74psidB, overlap), table(re78>0, trt))) > AS.glm <- glm(I(re78>0) ~ trt + Pscores, data=nsw74psidB, 9 subset = overlap, family=binomial) ? AS.lm <- lm(log(re78+47) ~ trt + Pscores, data=nsw74psidB, 4 subset = overlap&nsw74psidB$re78>0)  print(summary(AS.glm)$coef)  print(summary(AS.lm)$coef)  print(sum(overlap)) . trt <- factor(trt, labels=c("ctl", "trt")) 3 linecols <- c(rgb(0,0.5,0.5), rgb(0.65,0,0.65)) K if(bw) trellis.par.set(superpose.symbol=list(col=c("gray40", "gray40"), ) pch=c(1,3)), F superpose.line=list(col=c("black", "gray"), 5 lty=c(1,1), lwd=c(1,2)), = layout.heights=list(axis.top=0.6))  else 8 trellis.par.set(superpose.symbol=list(pch=c(1,3)), C superpose.line=list(lty=c(1,1), col=linecols, ( lwd=c(1.5,1.5))) / print(trellis.par.get()$superpose.line$lty) / gph <- xyplot(resid(AS.lm) ~ fitted(AS.lm), 9 groups=trt[overlap&nsw74psidB$re78>0], ' type=c("p","smooth"), A auto.key=list(columns=2, lines=TRUE), aspect=1) 7 print(gph, position=c(0.5, 0, 1, 1), newpage=FALSE)  if(device!="")dev.off()  }BC m@> maxf.score@>n " \D\DDD  srcref ]^FF  srcref _`,,  srcref aa  srcref bb!!  srcref cd99  srcref eh//  srcref ii    srcref jj&&  srcref kk  srcref ll  srcref mn22  srcref oo  srcref pr>>  srcref ss44  srcref tt  srcref uu  srcref vv  srcref ww88  srcref xx55  srcref yy66  srcref zz@@  srcref {|55  srcref }~44  srcref   srcref   srcref   srcref ..  srcref 33  srcref ((  srcref //  srcref AA  srcref 77  srcref   srcref  @@q G@ @ color(n @@ffffff@@=?ffffffJ@?K s rsmm tcu$C%@ @$&'$vw'@ function(x)unique(sort(x))[2]  disc.glm glm1S@~t~tr B family binomial  Pscores predict  nsw74psidB$Cr SS_??n E black gray E trellis.par.gety3e   xchop overlapDensity$RS$RS?D?? compare.numbers 3E ratio /?[\@]@ adj?ə A=?陙e _???a 0%?%?7 8  axes X n[\@]@?ffffff B=?陙  overlap &$? <$@ with tableS  AS.glm1S   AS.lm}1~@GS  sum SS ctl trt  linecols rgb????nxyz3 gray40 gray406?@{z3 black grayD??F?@ layout.heightsz axis.top?333333xyz6?@{zD??3F??{D  gph1$SX p smoothz@Z  aspect? position??? newpage  ? gdump 7function(fnam=NULL, prefix="~/r-book/ed2/figures/figs", K xtras=c("hardcopy","renum.obj","renum.files"), splitchar="/ch"){  if(is.null(fnam)){  path <- getwd() 7 pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] 5 fnam <- paste(prefix, pathtag[length(pathtag)], ! ".R", sep="")  } - else fnam <- paste(prefix, fnam, sep="/") C objnames <- c(objects(pattern="^g", envir=sys.frame(0)), xtras) & cat("\nDump to file:", fnam, "\n")  print(objnames)  dump(objnames, fnam)  } fnam prefix ~/r-book/ed2/figures/figs xtras hardcopy renum.obj renum.files splitchar /ch KKKK  srcref --  srcref CC  srcref &&  srcref   srcref   srcref  is.null   srcref   srcref 77  srcref !!  srcref   path getwd  pathtag [[ strsplit /ch fixed ?  paste$ length .R sep   /  objnames objects pattern ^g envir sys.frame cat  Dump to file:   dump gsave 7function(fnam=NULL, prefix="~/r-book/ed2/figures/figs", K splitchar="/ch", xtras=c("hardcopy","renum.obj","renum.files")){  if(is.null(fnam)){  path <- getwd() 7 pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] 5 fnam <- paste(prefix, pathtag[length(pathtag)], % ".RData", sep="")  } - else fnam <- paste(prefix, fnam, sep="/") C objnames <- c(objects(pattern="^g", envir=sys.frame(0)), xtras) & cat("\nDump to file:", fnam, "\n")  print(objnames) " save(list=objnames, file=fnam)  } ~/r-book/ed2/figures/figs /ch hardcopy renum.obj renum.files KKKK  srcref --  srcref CC  srcref &&  srcref   srcref ""  srcref    srcref   srcref 77  srcref %%  srcref    /ch ? $ .RData   /  ^g  Dump to file:   savez file4 =function(width=3.75, height=3.75, color=FALSE, trellis=FALSE, H device=c("","pdf","ps"), path="~/r-book/ed2/Art/", file=NULL, D format=c("nn-nn", "name"), split="\\.", pointsize=c(8,4),  fonts=NULL,  horiz=FALSE, ...){ ) if(!trellis)pointsize <- pointsize[1]  funtxt <- sys.call(1) B nam <- strsplit(as.character(funtxt), "(", fixed=TRUE)[[1]][1] 3 suffix <- switch(device, ps=".eps", pdf=".pdf") N if(is.character(path) & nchar(path)>1 & substring(path, nchar(path))!="/") & path <- paste(path, "/", sep="") - if(is.null(file)) if(format[1]=="nn-nn"){ C if(!is.null(split))dotsplit <- strsplit(nam, split)[[1]] else  dotsplit <- nam 8 if(length(dotsplit)==1)dotsplit <- c("", dotsplit) E nn2 <- paste(if(nchar(dotsplit[2])==1)"0" else "", dotsplit[2],  sep="")  if(nchar(dotsplit[1])>0){ O numstart <- which(unlist(strsplit(dotsplit[1], "")) %in% paste(0:9))[1] / nn1 <- substring(dotsplit[1], numstart) G nn1 <- paste(if(nchar(nn1) == 1) "0" else "", nn1, "-", sep="")  } else nn1 <- "" % file <- paste(nn1, nn2, sep="")  } else file <- nam L if(nchar(file)>4 & substring(file, nchar(file)-nchar(suffix)+1)==suffix)  suffix <- "" - file <- paste(path, file, suffix, sep="") : print(paste("Output will be directed to file:", file))  dev.out <- device[1] 6 dev.fun <- switch(dev.out, pdf=pdf, ps=postscript)  if(trellis){  library(lattice)  if(device=="ps") 1 trellis.device(file=file, device=dev.fun, ? color = color, horiz=horiz, fonts=fonts, < width=width, height=height, ...) else < trellis.device(file=file, device=dev.fun, fonts=fonts, D color = color, width=width, height=height, ...) R trellis.par.set(list(fontsize=list(text=pointsize[1], points=pointsize[2]))) } else  if (dev.out!=""){  print(c(width, height))  if(device=="ps") F dev.fun(file=file, paper="special", horiz=horiz, fonts=fonts, M width=width, height=height, pointsize=pointsize[1], ...) else 6 dev.fun(file=file, paper="special", fonts=fonts, L width=width, height=height, pointsize=pointsize[1], ...)  } D if(trellis)trellis.par.set(list(fontsize=list(text=pointsize[1], < points=pointsize[2])))  }@@ q   pdf ps ~/r-book/ed2/Art/ format nn-nn name split \.G@ @ fontsg  ...   srcref ))  srcref   srcref BB  srcref 33  srcref &&  srcref   srcref   srcref --  srcref ::  srcref   srcref 66  srcref   srcref <<  srcref (q G$G?  funtxt sys.call?  nam$ as.character ( ??  suffix switch ps .eps pdf .pdf is.character nchar? substring /  / R$? nn-nn ----  srcref   srcref 88  srcref   srcref   srcref %%  srcref (  dotsplit? R?    nn2R$@? 0 $@ $?   srcref  O O  srcref  / /  srcref  G G  srcref   numstart$ which %in% unlist$? %@"?  nn1$? R? 0  -      @R?     Output will be directed to file:  dev.out$?  dev.fun postscriptq   srcref   srcref DD  srcref RR  srcref o latticeR ps trellis.deviceggxz fontsizez text$G? points$G@    srcref   srcref FF  srcref R ps paper specialggG$G? specialG$G?qxzz$G?$G@ renum.obj+ 9function(from.prefix="g", to.prefix="g",from=2:5, to=1:4, L prenum=paste(rev(strsplit(path, "/ch", fixed=TRUE)[[1]])[1], ".", ? sep=""), details=FALSE, doit=FALSE, archive=TRUE){ " if(length(from) != length(to)) J stop(paste("to vector has length ", to, "; from vector has length ",  from,sep=""))  path <- getwd() 6# pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] '# endbit <- pathtag[length(pathtag)] 5 from.prefix <- paste(from.prefix, prenum, sep="") 1 to.prefix <- paste(to.prefix, prenum, sep="") 0 namsfrom <- paste(from.prefix, from, sep="") * namsto <- paste(to.prefix, to, sep="") ) tmpfrom <- paste("tmp", from, sep="") 6 from2to <- paste(namsto, " <- ", namsfrom, sep="") & cat(paste(from2to, collapse="\n")) cat("\n") if(doit)  for(i in 1:length(to)){ ! if(!exists(namsfrom[i])) A stop(paste("From file", namsfrom[i], "was not found"))  if(!exists(namsto[i])) = stop(paste("To file", namsto[i], "was not found")) C from2tmp <- paste(tmpfrom[i], " <- ", namsfrom[i], sep="") # if(details)print(from2tmp) 7 eval(parse(text=from2tmp), envir=sys.frame(0))  }  if(doit)  for(i in 1:length(to)){ @ tmp2to <- paste(namsto[i], " <- ", tmpfrom[i], sep="") 8 rmtxt <- paste("rm(", tmpfrom[i], ")", sep="")  if(details){  print(tmp2to)  print(rmtxt) } 6 eval(parse(text=tmp2to), envir=sys.frame(0)) 5 eval(parse(text=rmtxt), envir=sys.frame(0)) }  if(archive & doit){ / save(list=namsfrom, file="objfrom.RData") N print("From objects have been archived to 'objfrom.RData'", quote=FALSE)  }  } from.prefix g to.prefix g from%@@ to%?@ prenum$ rev /ch ?? .  details  doit  archive  ????  srcref   srcref   srcref 55  srcref 11  srcref 00  srcref **  srcref ))  srcref 66  srcref &&  srcref     srcref   srcref      srcref   srcref  stop to vector has length  ; from vector has length         namsfrom   namsto   tmpfrom tmp   from2to  <-   collapse    MN%?   srcref  A A  srcref  = =  srcref  C C  srcref  # #  srcref  7 7  srcref ( exists$N From file$N was not found($N To file$N was not found  from2tmp$N  <- $N  eval parseMN%?     srcref   @ @  srcref   8 8  srcref     srcref  6 6  srcref  5 5  srcref   tmp2to$N  <- $N   rmtxt rm($N )      srcref     srcref     srcref    srcref //  srcref NN  srcref z objfrom.RData 2From objects have been archived to 'objfrom.RData' quote  renum.files Hfunction(from.prefix="~/r-book/ed2/Art/", to.prefix="~/r-book/ed2/Art/", ' from=20:7, to=21:8, doit=F){  path <- getwd() 5 pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] & endbit <- pathtag[length(pathtag)] 2 if(nchar(endbit)==2)chap <- paste(endbit) else $ chap <- paste("0",endbit,sep="") 8 from.prefix <- paste(from.prefix, chap, "-", sep="") 4 to.prefix <- paste(to.prefix, chap, "-", sep="")  for(i in 1:length(from)){ 8 if (from[i]<=9) ltext <- paste("0",from[i],sep="") " else ltext <- paste(from[i]) 4 if (to[i]<=9) rtext <- paste("0",to[i],sep="") , else rtext <- paste(to[i]) 8 txt<-paste("mv ", from.prefix, ltext, ".eps", " ", 2 to.prefix, rtext, ".eps", sep="") ; backup<-paste("cp ", from.prefix, ltext, ".eps", " ", & "archive", sep="")  if(doit)system(backup)  if(doit)system(txt)  print(backup)  print(txt)  }  } ~/r-book/ed2/Art/ ~/r-book/ed2/Art/%@4@%@5@  F ''''  srcref   srcref 55  srcref &&  srcref  !$$  srcref ""88  srcref ##44  srcref $1  srcref    /ch ?  endbit$R @  chap    0     -    - MN%?  $$  srcref %&""  srcref '(    srcref )*22  srcref +,&&  srcref --  srcref ..  srcref //  srcref 00  srcref  <=$N@"  ltext 0$N  $N $N@"  rtext 0$N  $N  txt mv  .eps   .eps   backup cp  .eps   archive  system