#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)

# Looking at only below age 18. Otherwise growth stops and we do not have an increasing function any more.
nh_young <- subset(NHANES, Age <= 18)

xyplot(Height~Age | Race1*Gender, data = nh_young)
#Females stop growing earlier than 18. For males, taking 18 as upper limit is okay.   
#There are more data points for whites than for other races.


## 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)
}
# All categories combined using our function
ab_young=with(nh_young,lsfit(Age,Height))
ab_young

# Compare results with built-in function lm of R
coef(lm(Height ~ 1 + Age, data =nh_young ))
plot(Height~Age, data=nh_young)
abline(ab_young)

#Genders by different color in same plot 
# "p" for points, "r" for regression line
# auto.key creates legend.
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"),
      )


# Concentrating on white males, as the assumption of linear model seems to hold and it is the largest group. 
nhSubYMW <- subset(NHANES, Race1 == "White" & Age <= 18 & Gender == "male")
abSubYMW=with(nhSubYMW,lsfit(Age,Height))
xyplot(Height ~ Age, data = nhSubYMW, abline=abSubYMW)
abSubYMW

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