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)