#This is survey data collected by the US National Center for Health Statistics (NCHS) which has conducted a series of health and nutrition surveys since the early 1960's. Since 1999 approximately 5,000 individuals of all ages are interviewed in their homes every year and complete the health examination component of the survey. The health examination is conducted in a mobile examination centre (MEC).

library(NHANES)
library(lattice)

nh_young <- subset(NHANES, Age <= 18)

xyplot(Height~Age | Race1*Gender, data = nh_young)


## function to calculate a-hat and b-hat

s2 <- function(x) {
  m <- mean(x)
  mean((x - m)^2)
}

cxy <- function(x, y) {
  mx <- mean(x)
  my <- mean(y)
  mean((x - mx) * (y - my))
}

lsfit <- function(x, y) {
  if (length(x) != length(y)) stop("Lengths and x and y do not match")
  skip <- is.na(x) | is.na(y)
  x <- x[!skip]
  y <- y[!skip]
  mx <- mean(x)
  my <- mean(y)
  s2x <- s2(x)
  covxy <- cxy(x, y)
  b_hat <- covxy / s2x
  a_hat <- my - b_hat * mx
  c(a = a_hat, b = b_hat)
}

ab_young=with(nh_young,lsfit(Age,Height))
ab_young
coef(lm(Height ~ 1 + Age, data =nh_young ))

plot(Height~Age, data=nh_young)
abline(ab_young)

xyplot(Height ~ Age, data = nh_young, grid = TRUE, 
       groups = Gender, auto.key = TRUE, 
       type = c("p", "r"))

xyplot(Height ~ Age, data = nh_young, grid = TRUE, 
       groups = Race1, auto.key = TRUE, 
       type = c("p", "r"),
      )



nhSubYMW <- subset(NHANES, Race1 == "White" & Age <= 18 & Gender == "male")
abSubYMW=with(nhSubYMW,lsfit(Age,Height))
xyplot(Height ~ Age, data = nhSubYMW, abline=abSubYMW)
abSubYMW


lmSubYMW=lm(nhSubYMW$Height~nhSubYMW$Age)
summary(lmSubYMW)
anova(lmSubYMW)
par(mfrow=c(2,2))
plot(lmSubYMW)