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

Caroline Ring

2018-01-23

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)