# Ring et al. 2017 Height and weight spline fits and residuals

#### Caroline Ring

#### 2019-04-25

    library(httk)
    library(survey)
    library(data.table)
    library(ks)
    library(ggplot2)
    all.reth <- levels(nhanes_mec_svy$variables[, ridreth1])
    all.gendr <- levels(nhanes_mec_svy$variables[, riagendr])

Plot smoothing spline fits to height.

    par(mfcol=c(5,2),
        mar=c(2,2,2.5,2)+0.1,
        oma=c(4,4,4,0),
        mgp=c(1,1,0))

    for (gendr in all.gendr){
      for (reth in all.reth){
        grsub <- subset(nhanes_mec_svy, riagendr==gendr & ridreth1==reth)
        svyplot(logbmxhtlenavg~ridexagm, grsub, style="transparent",
                basecol="gray50", xlab=NA, ylab=NA, cex.axis=1.5)
        lines(sort(unique(grsub$variables[,ridexagm])), predict(spline_heightweight[g==gendr & r==reth,
                                                                   height_spline[[1]]],
                                      sort(unique(grsub$variables[,ridexagm])))$y)
        title(paste(gendr, reth), line=0.5, cex.main=1.5)
      }
    }
    title(xlab="Age, months",
          ylab="Log height, log cm",
          outer=TRUE,
          line=2, cex.lab=2)
    title(main="Height vs. age with spline fits",
          outer=TRUE, line=1, cex.main=2)

Plot log height residuals vs. age to check for heteroscedasticity.

    par(mfcol=c(5,2),
        mar=c(2,2,2.5,2)+0.1,
        oma=c(4,4,4,0),
        mgp=c(1,1,0))

    for (gendr in all.gendr){
      for (reth in all.reth){
        grsub <- subset(nhanes_mec_svy, riagendr==gendr & ridreth1==reth & !is.na(ridexagm))
        grsub <- update(grsub, logbmxhtlenavg_resid=logbmxhtlenavg-
                          predict(spline_heightweight[g==gendr &
                                                        r==reth,
                                                      height_spline[[1]]],
                                  grsub$variables[,ridexagm])$y)

        svyplot(logbmxhtlenavg_resid~ridexagm, grsub, style="transparent",
                xlab=NA, ylab=NA, cex.axis=1.5)
        title(paste(gendr, reth), line=0.5, cex.main=1.5)
        }
      }
    title(xlab="Age, months",
          ylab="Log height residual, log cm",
          outer=TRUE,
          line=2, cex.lab=2)
    title(main = "Height residuals vs. age",
          outer=TRUE,
          line=1,
          cex.main=2)

In general, it looks like residuals have fairly constant variance with age.

Plot smoothing spline fits to weight.

    par(mfcol=c(5,2),
        mar=c(2,2,2.5,2)+0.1,
        oma=c(4,4,4,0),
        mgp=c(1,1,0))

    for (gendr in all.gendr){
      for (reth in all.reth){
        grsub <- subset(nhanes_mec_svy, riagendr==gendr & ridreth1==reth & !is.na(ridexagm))
        svyplot(logbmxwt~ridexagm, grsub, style="transparent",
                basecol="gray50", xlab=NA, ylab=NA, cex.axis=1.5)
        lines(sort(unique(grsub$variables[,ridexagm])), predict(spline_heightweight[g==gendr & r==reth,
                                                                   weight_spline[[1]]],
                                      sort(unique(grsub$variables[,ridexagm])))$y)
        title(paste(gendr, reth), line=0.5, cex.main=1.5)
      }
    }
    title(xlab="Age, months",
          ylab="Log weight, log kg",
          outer=TRUE,
          line=2, cex.lab=2)
    title(main = "Weight vs. age with spline fits",
           outer=TRUE,
           line=1,
           cex.main=2)

Plot log weight residuals vs. age.

    par(mfcol=c(5,2),
        mar=c(2,2,2.5,2)+0.1,
        oma=c(4,4,4,0),
        mgp=c(1,1,0))

    for (gendr in all.gendr){
      for (reth in all.reth){
        grsub <- subset(nhanes_mec_svy, riagendr==gendr & ridreth1==reth & !is.na(ridexagm))
        grsub <- update(grsub, logbmxwt_resid=logbmxwt-
                          predict(spline_heightweight[g==gendr &
                                                        r==reth,
                                                      weight_spline[[1]]],
                                  grsub$variables[,ridexagm])$y)
        svyplot(logbmxwt_resid~ridexagm, grsub, style="transparent",
                xlab=NA, ylab=NA, cex.axis=1.5)
        title(paste(gendr, reth), line=0.5, cex.main=1.5)
        }
      }
    title(xlab="Age, months",
          ylab="Log weight residual, log kg",
          outer=TRUE,
          line=2, cex.lab=2)
    title(main = "Weight residuals vs. age",
          outer=TRUE, line=1, cex.main=2)

Again, in general, it looks like residuals have fairly constant variance with age.

We can expect log height residuals and log weight residuals to have some covariance. In other words, someone with above-average weight for their age is also more likely to have above-average height for their age. To capture this correlation in the virtual-individuals method, we fit a 2-dimensional KDE to the height and weight residuals. To visualize this, we plot log height residuals vs. log weight residuals, and overlay contours of the joint KDE.

    par(mfcol=c(5,2),
        mar=c(2,2,2.5,2)+0.1,
        oma=c(4,4,4,0),
        mgp=c(1,1,0))
    for (gendr in all.gendr){
      for (reth in all.reth){
        grsub <- subset(nhanes_mec_svy, riagendr==gendr & ridreth1==reth & !is.na(ridexagm))
        grsub <- update(grsub, logbmxwt_resid=logbmxwt-
                          predict(spline_heightweight[g==gendr &
                                                        r==reth,
                                                      weight_spline[[1]]],
                                  grsub$variables[,ridexagm])$y)
        grsub <- update(grsub, logbmxhtlenavg_resid=logbmxhtlenavg-
                          predict(spline_heightweight[g==gendr &
                                                        r==reth,
                                                      height_spline[[1]]],
                                  grsub$variables[,ridexagm])$y)
        hw_kde <- kde(x=as.matrix(grsub$variables[, .(logbmxwt_resid, logbmxhtlenavg_resid)]),
                                w=grsub$variables[, wtmec6yr])
         svyplot(logbmxhtlenavg_resid~logbmxwt_resid, grsub, style="transparent",
                 basecol="gray50",
                xlab=NA, ylab=NA, cex.axis=1.5)
        plot(hw_kde, cont=c(10,25,50,75,90,95), display="slice",
             col="black", labcex=1.1, vfont=c("sans serif", "bold"), add=TRUE)
        title(paste(gendr, reth), line=0.5, cex.main=1.5)
        }
      }
       title(xlab="Log weight residual, log kg",
          ylab="Log height residual, log cm",
          outer=TRUE,
          line=2, cex.lab=2)
    title(main = "Height resid vs. weight resid with KDE",
          outer=TRUE, line=1, cex.main=2)