Chapter 10: Multi-level Models, and Repeated Measures
Packages required: “DAAG”, “lattice”, “MEMSS”, “lme4”, “grid”

The script that follows is designed to be executed as it stands. It sets up
functions that create Ch 10 graphs, then runs those functions. If executed by
clicking on the RStudio “Compile Notebook …” command, it will be processed
through R Markdown, generating a document that includes code, output and
graphs generated by executing the functions. If some needed packages are
missing, warning messages will appear. See further:
http://rmarkdown.rstudio.com/r_notebook_format.html

g10.1 <-
function(device=""){
    if(device!="")hardcopy(width=3.25, height=2.5,
                           device=device, trellis=TRUE,
                           color=FALSE, pointsize=c(9,6))
    oldpar <- par(mar = par()$mar - c(1, 0, 1, 0))
    on.exit(par(oldpar))
    if(!require(lattice))return("Package 'lattice' must be installed.")        
    Site <- with(ant111b, reorder(site, harvwt, FUN=mean))
    print(stripplot(Site ~ harvwt, data=ant111b, scales=list(tck=0.5),
                    xlab="Harvest weight of corn"))
    if(device!="")dev.off()
  }

g10.2 <-
function(device=""){
    if(device!="")hardcopy(width=5.5, height=1.8, pointsize=c(8,4),
                           device=device, trellis=TRUE, color=FALSE)
    if(!require(DAAG))return("Package 'DAAG' must be installed.") 
    if(!require(lme4))return("Package 'lme4' must be installed.")  
    if(!require(lattice))return("Package 'lattice' must be installed.")     
    trellis.par.set(layout.heights=list(key.top=0.5, axis.top=1.15,
                  bottom.padding=0.15, main.key.padding=2.5))
    ant111b.lmer <- lmer(harvwt ~ 1 + (1|site), data=ant111b)
    prof.lmer <- profile(ant111b.lmer)
    if(!require(lattice))return("Package 'lattice' must be installed.")        
    print(xyplot(prof.lmer, conf=c(50, 80, 95, 99)/100, aspect=0.8,
          between=list(x=0.35)))
    if(device!="")dev.off()
  }

g10.3 <-
function(device="")
{
  if(device!="")hardcopy(width=3, height=1.85, device=device)
  oldpar <- par(mar = c(4.1, 6.1,1.1,1.1), mgp=c(2.5, 0.75,0), las=2)
  on.exit(par(oldpar))
  require(DAAG)
  classmeans <- with(science, aggregate(like,by=list(PrivPub,Class),mean))
  names(classmeans) <- c("PrivPub","Class","like")
  with(classmeans, {
    boxplot(split(like, PrivPub), xlab = "Class average of score",
            boxwex = 0.4, horizontal=TRUE)
    rug(like[PrivPub == "private"], side = 1)
    rug(like[PrivPub == "public"], side = 3)
})
  if(device!="")dev.off()
}

g10.4 <-
function(pnum=c(1,3), device="", term = "school:class"){
                                        #    leg <- c("Class effect vs #", "W/i class var vs #",
                                        #             "qnorm(); site effect", "qnorm(); w/i class resids")
    if(!require(DAAG))return("Package 'DAAG' must be installed.") 
    if(!require(lme4))return("Package 'lme4' must be installed.")   
    leg <- rep("",4)
    mtext3 <- function(item="A", txt=leg[1],
                       xleft=par()$usr[1]+1.25*par()$cxy[1]){
      mtext(side = 3, line = 0.4, item, cex=0.925, adj = 0)
      mtext(side = 3, line = 0.4, txt, at=xleft, adj = 0)
    }
    science$class <- factor(science$class)
    science$school <- factor(science$school)
    if(device!="")hardcopy(width=3.75, height=3.75,
                           device=device, pointsize=9)
    science <- na.omit(science)
    science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class),
                          data = science, na.action=na.exclude)
    oldpar <- par(mfrow=c(2,2), mar = c(4.6, 4.1, 2.6, 2.1), mgp=c(2.5,0.75,0))
    on.exit(par(oldpar))
    ranf <- ranef(obj = science1.lmer, drop=TRUE)[[term]]
    flist <- science1.lmer@flist[[term]]
    privpub <- science[match(names(ranf), flist), "PrivPub"]
    num <- unclass(table(flist))
    if(!all(names(num)==names(flist)))stop("!all(names(num)==names(flist))")
    plot(sqrt(num), ranf, xaxt="n", pch=pnum[as.numeric(privpub)],
         xlab="# in class (square root scale)",
         ylab="Estimate of class effect")
    lines(lowess(sqrt(num[privpub=="private"]),
                 ranf[privpub=="private"], f=1.1), lty=2)
    lines(lowess(sqrt(num[privpub=="public"]),
                 ranf[privpub=="public"], f=1.1), lty=3)
    numlabs <- pretty(num)
    axis(1, at=sqrt(numlabs), labels=paste(numlabs))
    mtext3(item="A", txt=leg[1], xleft=par()$usr[1]+1.25*par()$cxy[1])
    res <- residuals(science1.lmer)
    vars <- tapply(res, INDEX=list(flist),
                      FUN=var)
    vars <- vars*(num-1)/(num-2)
    plot(sqrt(num), vars, xaxt="n", pch=pnum[as.numeric(privpub)],
         xlab="# in class (square root scale)",
         ylab="Within class variance")
    lines(lowess(sqrt(num[privpub=="private"]),
                 as.vector(vars)[privpub=="private"], f=1.1), lty=2)
    lines(lowess(sqrt(num[privpub=="public"]),
                 as.vector(vars)[privpub=="public"], f=1.1), lty=3)
    axis(1, at=sqrt(numlabs), labels=paste(numlabs))
    mtext3(item="B", txt=leg[2], xleft=par()$usr[1]+1.25*par()$cxy[1])
    qqnorm(ranf, ylab="Ordered site effects", main="")
    mtext3(item="C", txt=leg[3], xleft=par()$usr[1]+1.25*par()$cxy[1])
    qqnorm(res, ylab="Ordered w/i class residuals", main="")
    mtext3(item="D", txt=leg[4], xleft=par()$usr[1]+1.25*par()$cxy[1])
    par(mfrow=c(1,1))
    par(fig=c(0,1,0.8,1), mar=c(0,0,0,0), mgp=c(0,0,0),
        oma=c(0,0,0.5,0), new=TRUE)
    plot.new()
    plot.window(xlim=c(0,1), ylim=c(0,1))
    legend(x="top", legend=c("Private     ", "Public"), pch=c(1,3),
           lwd=c(1,1), lty=2:3,
           xjust=0.5, yjust=0, horiz=TRUE, merge=FALSE, bty="n")
    if(device!="")dev.off()
  }

g10.5 <-
function(device=""){
    if(device!="")hardcopy(width=4.25, height=4.25, device=device)
    oldpar <- par(mar=c(0.5,1,1,1),mgp=c(0,.5,.5))
    on.exit(par(oldpar))
    plot(c(0,13),c(0,13),type="n",xlab="",ylab="",
         axes=F)
    eps <- 0.1
    vines<-function(x=1,y=1,subp=0){
      lines(c(x,x+1,x+1,x,x),c(y,y,y+1,y+1,y))
      points(c(x+.2,x+.8,x+.8,x+.2),c(y+.2,y+.2,y+.8,y+.8),pch=3)
      text(x+.5,y+.5,paste(subp))
    }
    k<-0
    for(i in c(1,3,5,7)){k<-k+1; vines(1,i,c(3,1,0,2)[k])}
    k<-0
    for(i in c(1,3,5,7)){k<-k+1; vines(4,i,c(2,1,0,3)[k])}
    k <- 0
    for(i in c(1,4,4,1)){k<-k+1
                         j<-c(9,9,11,11)[k]
                         vines(i,j,c(3,2,1,0)[k])
                       }
    lines(c(3,3,NA,3,3),c(0,2.85,NA,10.15,13),lty=2)
    lines(c(0,0,NA,0,0),c(0,2.85,NA,10.15,13),lty=2)
    lines(c(8,8,NA,8,8),c(0,4.5,NA,8.5,13),lty=2)
    lines(c(0,1.25,NA,6.75,8),rep(0,5),lty=2)
    lines(c(0,1.25,NA,6.75,8),rep(13,5),lty=2)
    lines(c(1,5,5,1,1)+c(-eps,eps,eps,-eps,-eps),
          c(9,9,12,12,9)+c(-eps,-eps,eps,eps,-eps),lwd=2)
    lines(c(1,2,2,1,1)+c(-eps,eps,eps,-eps,-eps),
          c(1,1,8,8,1)+c(-eps,-eps,eps,eps,-eps),lwd=2)
    lines(c(1,2,2,1,1)+3+c(-eps,eps,eps,-eps,-eps),
          c(1,1,8,8,1)+c(-eps,-eps,eps,eps,-eps),lwd=2)
    text(0,6.5,"6 meters height artifical shelter belt", srt=90)
    text(4,0,"9 meters height shelter belt")
    text(4,13,"19 meters height shelter belt")
    text(8,6.5,"Willow shelter belt",srt=90)
    text(8.5,10.5,"0 Unshaded   \n1 Shaded Aug-Dec   \n2 Dec-Feb   \n3 Feb-May",
         adj=0)
    text(3,6.5,"16 meters height willow shelter belt", srt=90)
    lines(c(8.25,8.75),c(13,13-sqrt(3)*.5))

    lines(rep(8.25,2),c(12.6,13))
    lines(c(8.25,8.25+sqrt(3)/5),c(13,12.8))
    text(8.85,12.65,"N")
    if(device!="")dev.off()
  }

g10.6 <-
function(device="")
{
  if(device!="")hardcopy(width=3.75, height=2.25, device=device)
  oldpar<-par(mar=c(4.5,6.1,1.1,1.1), mfrow=c(1,1))
  on.exit(par(oldpar))
  if(!require(DAAG))return("Package 'DAAG' must be installed.")    
  kiwimeans <- with(kiwishade, aggregate(yield,by=list(block,shade),mean))
  names(kiwimeans) <- c("block","shade","meanyield")
  plotmeans.lm <- lm(meanyield~block+shade, data=kiwimeans)
  effects <- predict(plotmeans.lm, type="terms")
  kiwishade.lm <- lm(yield ~ block*shade, data=kiwishade)
  y <- c(effects[,"block"]/sqrt(2) * sqrt(16),
         effects[,"shade"]/sqrt(3) * sqrt(12),
         residuals(plotmeans.lm)/sqrt(6) * sqrt(4),
         residuals(kiwishade.lm)/
         sqrt(12))
  n <- rep(4:1, c(12, 12, 12, 48))
  gps <- rep(c("block effect\n(ms=86.2)", "shade\n(464.8)",
               "plot\n(20.9)", "vine\n(12.2)"), c(12, 12, 12, 48))
  gps <- factor(gps, levels = rev(c("block effect\n(ms=86.2)",
                       "shade\n(464.8)", "plot\n(20.9)",
                       "vine\n(12.2)")))
  gps.sd <- sapply(split(y,gps), sd)
  gps.av <- sapply(split(y,gps), mean)
  plot(range(y), range(n)+c(-0.5, 0.5), xlab="", ylab="", yaxt="n", type="n")
  text(y, n+0.15, "|")
  un <- 1:4
  sapply(un, function(j){lines(gps.av[j]+c(-gps.sd[j], gps.sd[j]),
                              rep(j-0.15,2), col="gray")
                        lines(rep(gps.av[j],2)-gps.sd[j],
                              j-0.15+c(-.06, .06), col="gray")
                        lines(rep(gps.av[j],2)+gps.sd[j],
                              j-0.15+c(-.06, .06), col="gray")
                      })
  mtext(side=1,line=3,
        text="Variation in Yield (kg)\n(Add to grand mean of yield = 96.53)")
  par(las=2)
  axis(2, at=1:4, labels=levels(gps))
  par(las=0)
  if(device!="")dev.off()
}

g10.7 <-
function(device="", path="")
{
  if(device!="")hardcopy(width=4, height=4,
                         device=device, trellis=TRUE, pointsize=c(7, 5))
      if(!require(lattice))return("Package 'lattice' must be installed.")    
      if(!require(grid))return("Package 'grid' must be installed.")    
      if(!require(lme4))return("Package 'lme4' must be installed.")    
  trellis.par.set(layout.heights=list(key.top=0.5, axis.top=1.15,
                    bottom.padding=0.15, main.key.padding=2.5))
  if (!exists("kiwishade.lmer"))
    kiwishade.lmer <- lmer(yield ~ shade + (1|block/shade), data=kiwishade)
  pk1 <- xyplot(residuals(kiwishade.lmer) ~ fitted(kiwishade.lmer)|block,
                groups=shade, layout=c(3,1), par.strip.text=list(cex=1.0),
                pch=1:4, data=kiwishade, grid=TRUE,
                xlab="Fitted values (Treatment + block + plot effects)",
                ylab="Residuals",
                scales=list(x=list(alternating=FALSE), tck=0.35),
                legend=list(top=list(fun=textGrob,
                              args=list(label="A", x=0,
                                just="left"))),
                key=list(x=0.052, y=1.24, points=list(pch=1:4),
                  text=list(labels=levels(kiwishade$shade)),columns=4))
  print(pk1,  position=c(0,.52,1,1))
  ploteff <- ranef(kiwishade.lmer, drop=TRUE)[[1]]
  nam <- names(sort(ploteff, decreasing=TRUE)[1:4])
  trellis.par.set(layout.heights=list(key.top=0.5, axis.top=0.85,
                    bottom.padding=0.15, main.key.padding=1))
  pk2 <- qqmath(ploteff, ylab="Plot effect estimates",
                scales=list(tck=0.35),
                key=list(x=0, y=0.98, corner=c(0,1),
                  text=list(labels=nam, cex=0.75)),
                legend=list(top=list(fun=textGrob,
                              args=list(label="B", x=0,
                                just="left"))), aspect=1,
                xlab="Normal quantiles")
  print(pk2,  newpage=FALSE, position=c(0,0,.5,.5))
  simvals <- simulate(kiwishade.lmer, nsim=2)
  simeff <- apply(simvals, 2, function(y) scale(ranef(refit(kiwishade.lmer, y),
                                                      drop=TRUE)[[1]]))
  simeff <- data.frame(v1=simeff[,1], v2=simeff[,2])
  pk3 <- qqmath(~ v1+v2, data=simeff,
                ylab="Simulated plot effects\n(2 sets, standardized)",
                Qs=list(tck=0.35), aspect=1,
                legend=list(top=list(fun=textGrob,
                              args=list(label="C", x=0,
                                just="left"))),
                xlab="Normal quantiles")
  print(pk3,  newpage=FALSE, position=c(0.48,0,1,.5))
  if(device!="")dev.off()
  invisible()
}

g10.8 <-
function(device=""){
    if(!require(lattice))return("Package 'lattice' must be installed.")     
#    if(!require(nlme))return("Package 'nlme' must be installed.")   
    if(!require(lme4))return("Package 'lme4' must be installed.")     
    if(device!="")hardcopy(device=device, width=3.6, height=3.4,
                           color=FALSE, pointsize=9)
    oldpar <- par(mgp=c(2.0, .5,0), mar=c(4.1,4.1,1.6,1.1),
                  pty="s", cex=1, cex.lab=1)
    on.exit(par(oldpar))
    par(fig=c(0, 0.925, 0, 1))
    plot(o2 ~ wattsPerKg, data=humanpower1,
         pch=(1:5)[unlist(id)],
         col=palette()[unlist(id)],
         xlab="Watts per kilogram",
         ylab=expression("Oxygen intake ("*ml.min^{-1}*.kg^{-1}*")"))
    hp.lmList <- lmList(o2 ~ wattsPerKg|id, data=humanpower1)
    sapply(1:length(hp.lmList), function(i)abline(hp.lmList[[i]], col=i))
    mtext(side=3, line=0.25, "A", adj=-0.165)
    par(fig=c(0.575, 1, 0.025, 0.525), new=TRUE, cex=0.8, mgp=c(1.8,0.5,0))
    coefs <- data.frame(t(sapply(hp.lmList, coef)))
    names(coefs) <- c("Intercept", "Slope")
    par(pty="s", tck=-0.025, xpd=TRUE, cex.axis=0.8, cex.lab=0.75,
        mgp=c(1.25,0.25,0))
    plot(Slope ~ Intercept, data=coefs, pch=1:5, type="n")
    xy <- par()$usr
    chw <- par()$cxy[1]
    chh <- par()$cxy[2]
    rect(xy[1]-3.25*chw,xy[3]-2.75*chh,xy[2]+chw,xy[4]+chh, col="gray85")
    axis(1)
    axis(2)
    box()
    with(coefs, points(Slope ~ Intercept, cex=1, pch=1:5))
    par(xpd=F)
    abline(lm(Slope ~ Intercept, data=coefs)$coef)
    par(xpd=TRUE, cex=1)
    text(xy[1]-2.25*chw, xy[4]+0.25*chh, "B")
    mtext(side=1,line=1.5,"Intercept", cex=0.8)
    mtext(side=2,line=1.5,"Slope", cex=0.8)
    par(fig=c(0,1,0,1), cex=1, xpd=FALSE)
    if(device!="")dev.off()
  }

g10.9 <-
function(device=""){
    if(!require(lattice))return("Package 'lattice' must be installed.")    
    if(!require(lme4))return("Package 'lme4' must be installed.")          
    if(device!="")hardcopy(device=device, width=5.5, height=2.75,
                           trellis=TRUE, color=FALSE)
    trellis.par.set(layout.heights=list(key.top=0.25, axis.top=0.65))
    hp.lmer <- lmer(o2 ~ wattsPerKg + (wattsPerKg | id),
                    data=humanpower1)
    hat2 <- fitted(hp.lmer)
    hp1.plt1 <- xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1,
                       panel=function(x,y,subscripts,groups,hat2,...){
                         u <- lm(y ~ groups*x);
                         hat <- fitted(u)
                         panel.superpose(x,hat,subscripts,groups, type="l")
                         panel.superpose(x, y=hat2, subscripts, groups,
                                         type="l", lty=2)
                       },
                       scales=list(tck=0.5),
                       xlab="Watts per kilogram",
                       ylab=expression("Oxygen intake ("*ml.min^{-1}*.kg^{-1}*")"),
                       legend=list(top=list(fun=textGrob, args=list(label="A", x=0, gp=gpar(cex=1.15)))),
                       hat2=hat2, aspect=1)
    print(hp1.plt1, position=c(0,0,0.5,1))
    res <- resid(hp.lmer)
    hp1.plt2 <- xyplot(res ~ wattsPerKg, groups=id,
                       data=humanpower1, xlab="Watts per kilogram",
                       ylab="Residuals from random lines",
                       scales=list(tck=0.5), aspect=1,
                       legend=list(top=list(fun=textGrob, args=list(label="B", x=0, gp=gpar(cex=1.15)))), type="l")
    print(hp1.plt2, position=c(0.5,0,1,1), newpage=FALSE)
    if(device!="")dev.off()
  }

g10.10 <-
function(device="", log=2){
    if(!require(lattice))return("Package 'lattice' must be installed.")       
    if(device!="")hardcopy(device=device, width=5.25, height=3.25,
                           trellis=TRUE, color=FALSE, pointsize=c(7,5))
    plt <- xyplot(distance ~ age | Subject, groups=Sex,
                  data=Orthodont, type=c("p","r"),
                  par.strip.text=list(cex=0.75), scale=list(log=log),
                  layout=c(11,3))
    print(plt)
    if(device!="")dev.off()
  }

g10.11 <-
function(device=""){
    if(device!="")hardcopy(device=device, width=4.5, height=2.25,
                           color=TRUE, pointsize=8)
    oldpar <- par(mar=c(4.1,4.1,1.1,1.1), mgp=c(2.5,0.75,0), pty="s",
                  mfrow=c(1,2))
    on.exit(par(oldpar))
    if(!require(lme4))return("Package 'lme4' must be installed.")  
    if(!require(MEMSS))return("Package 'MEMSS' must be installed.")        
    Orthodont$age <- Orthodont$age-11
    ab <- coef(lmList(distance ~ age|Subject, data=Orthodont))
    ab$sex <- substring(rownames(ab),1,1)
    names(ab) <- c("a", "b", "sex")
    xlim <- range(ab[,1])
    xlim <- xlim+diff(xlim)*c(-.015, .185)
    ylim <- range(ab[,2])
    ylim <- ylim+diff(ylim)*c(-.015, .015)
    plot(ab[,1], ab[,2], col=c(F="gray40", M="black")[ab$sex],
         pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
         xlab="Intercept", ylab="Slope")
    use <- ab$a %in% range(ab$a) | ab$b %in% range(ab$b) |
    ab$b==min(ab$b[ab$sex=="M"])
    text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
    mtext(side=3, line=1,"A: Distances", adj=0)
    ##
    Orthodont$logdist <- log(Orthodont$distance)
    ab <- coef(lmList(logdist ~ age|Subject, data=Orthodont))
    ab$sex <- substring(rownames(ab),1,1)
    names(ab) <- c("a", "b", "sex")
    xlim <- range(ab[,1])
    xlim <- xlim+diff(xlim)*c(-.015, .185)
    ylim <- range(ab[,2])
    ylim <- ylim+diff(ylim)*c(-.015, .015)
    plot(ab[,1], ab[,2], col=c(F="gray40", M="black")[ab$sex],
         pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
         xlab="Intercept", ylab="Slope")
    text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
    mtext(side=3, line=1,"B: Logarithms of distances", adj=0)
    if(device!="")dev.off()
  }

g10.11Col <-
function(device=""){
    if(device!="")hardcopy(device=device, width=4.5, height=2.25,
                           color=TRUE, pointsize=8)
    oldpar <- par(mar=c(4.1,4.1,1.1,1.1), mgp=c(2.5,0.75,0), pty="s",
                  mfrow=c(1,2))
    on.exit(par(oldpar))
    if(!require(lme4))return("Package 'lme4' must be installed.")    
    if(!require(MEMSS))return("Package 'MEMSS' must be installed.")        
    Orthodont$age <- Orthodont$age-11
    ab <- coef(lmList(distance ~ age|Subject, data=Orthodont))
    ab$sex <- substring(rownames(ab),1,1)
    names(ab) <- c("a", "b", "sex")
    xlim <- range(ab[,1])
    xlim <- xlim+diff(xlim)*c(-.015, .185)
    ylim <- range(ab[,2])
    ylim <- ylim+diff(ylim)*c(-.015, .015)
    plot(ab[,1], ab[,2], col=c(F="red", M="blue")[ab$sex],
         pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
         xlab="Intercept", ylab="Slope")
    use <- ab$a %in% range(ab$a) | ab$b %in% range(ab$b) |
    ab$b==min(ab$b[ab$sex=="M"])
    text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
    mtext(side=3, line=1,"A: Distances", adj=0)
    ##
    Orthodont$logdist <- log(Orthodont$distance)
    ab <- coef(lmList(logdist ~ age|Subject, data=Orthodont))
    ab$sex <- substring(rownames(ab),1,1)
    names(ab) <- c("a", "b", "sex")
    xlim <- range(ab[,1])
    xlim <- xlim+diff(xlim)*c(-.015, .185)
    ylim <- range(ab[,2])
    ylim <- ylim+diff(ylim)*c(-.015, .015)
    plot(ab[,1], ab[,2], col=c(F="red", M="blue")[ab$sex],
         pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
         xlab="Intercept", ylab="Slope")
    text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
    mtext(side=3, line=1,"B: Logarithms of distances", adj=0)
    if(device!="")dev.off()
  }

pkgs <- c("DAAG", "lattice", "MEMSS", "lme4", "grid")
z <- sapply(pkgs, require, character.only=TRUE, warn.conflicts=FALSE, quietly=TRUE)
## Loading required package: Matrix
## Loading required package: Rcpp
if(any(!z)){
  notAvail <- paste(names(z)[!z], collapse=", ")
  warning(paste("The following packages should be installed:", notAvail))
}

g10.1()

plot of chunk unnamed-chunk-1

g10.2()

plot of chunk unnamed-chunk-1

g10.3()

plot of chunk unnamed-chunk-1

g10.4()

plot of chunk unnamed-chunk-1

g10.5()

plot of chunk unnamed-chunk-1

g10.6()