R Notebook

Part 1

Lets first import the data:

bone <- read.csv("bone.csv", as.is=T)

** a) Begin with a plot of spnbmd against age, colorcoded by gender. **

plot (bone$age, bone$spnbmd,  col=ifelse(bone$gender=="male","blue","red"), pch = 8)

Does it appear that there are differences in bone mineral density (BMD) trajectories between males and females?

There does not appear to be an extreme difference between the BMD trajectories between males and females. But there is a bit more variablility of bone density accross the different ages for males compared to females. At a younger age the males seem to have lower BMD, and at around age 15 they seem to have a higer BMD compared to females. When becoming older though, at ages above 20, they appear to have the same trajectories.

** fit a series of models to spnbmd using a polynomial of age for each gender separately. Use knots at c1 and c2, where c1 and c2 are the 33rd and 67th percentiles of the data, respectively.**

install.packages("lsr")

Installing package into <U+393C><U+3E31>C:/Users/Nicole Salomons/Documents/R/win-library/3.3<U+393C><U+3E32>
(as <U+393C><U+3E31>lib<U+393C><U+3E32> is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.3/lsr_0.5.zip'
Content type 'application/zip' length 140820 bytes (137 KB)
downloaded 137 KB

package ‘lsr’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\Nicole Salomons\AppData\Local\Temp\RtmpUj5Bik\downloaded_packages

library(splines)
library(ISLR)
library(lsr)

bones2 = bone
bones2$agegrp <- quantileCut(bones2$age, 3)
labs <- levels(bones2$agegrp)
bds <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
      upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) ))
bds <- c(bds[,1], bds[nrow(bds),2])
bones.male = subset(bones2, gender == "male")
bones.female = subset(bones2, gender == "female")

** 1) Construct the model **

model.quadratic.male <- lm(spnbmd ~ agegrp*poly(age, 2, raw=T), data=bones.male)
model.quadratic.female <- lm(spnbmd ~ agegrp*poly(age, 2, raw=T), data=bones.female)

** 2) State the number of coeficients fitted, including the intercept **

length(coef(model.quadratic.male))

[1] 9

length(coef(model.quadratic.female))

[1] 9

The number of coeficients is always the same for both male and female. For a quadratic model there are 9 coeficients. A quadratic equation has 3 coeficentes. And each of the 3 age groups has its own quadratic model. So there is 3*3 = 9 coeficients.

** 3) Show a plot with the fitted models superimposed **

plot (bone$age, bone$spnbmd,  col=ifelse(bone$gender=="male","blue","red"), pch = 8)
abline(v=bds, lwd=2, lty=3, col="blue")

curve(predict(model.quadratic.male, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="blue", add=TRUE)
curve(predict(model.quadratic.female, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="red", add=TRUE)

The plot seems to show a lot of variabiity between the groups and jumps at the knots. This is not great because we would like to create a smoother model, which also does not have these extreme jumps at knots. Therefore the continuous models might be better options. The plot does show however that there is a clear difference between males and females before around the age of 18.

Continuous piecewise quadratic

** 1) Construct the model **

#bds
v <- seq(1, 10, by=.1)
#pmax(v-3, 0) 
model.quadratic.continuous.male <- lm(spnbmd ~ poly(age,2, raw=T) + 
            poly(pmax(I(age-bds[2]),0), 2, raw=T) +
            poly(pmax(I(age-bds[3]),0), 2, raw=T), data=bones.male)
model.quadratic.continuous.female <- lm(spnbmd ~ poly(age,3, raw=T) + 
            poly(pmax(I(age-bds[2]),0), 2, raw=T) +
            poly(pmax(I(age-bds[3]),0), 2, raw=T), data=bones.female)

** 2) State the number of coeficients fitted, including the intercept **

length(coef(model.quadratic.continuous.male))

[1] 7

Here the knots are continuous. So we lose 2 coeficients from the model, both at the knots. So instead of 9, now we have 7 coeficients.

** 3) Show a plot with the fitted models superimposed **

plot (bone$age, bone$spnbmd,  col=ifelse(bone$gender=="male","blue","red"), pch = 8)
abline(v=bds, lwd=2, lty=3, col="blue")

curve(predict(model.quadratic.continuous.male, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="blue", add=TRUE)
curve(predict(model.quadratic.continuous.female, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="red", add=TRUE)

This is a bit nicer, because it is continuous at the knots. But There are still sharp edges at the knots, which make it seem that at the ages of the knots drastic changes in spnbmd occur, which is not true. So the following models will take derivates, to have smoother knots.

Continuous piecewise quadratic with continuous first derivative

** 1) Construct the model **

model.quadratic.continuous.1derivative.male <- lm(spnbmd ~ poly(age,2, raw=T) + pmax(I(age-bds[2])^2,0) + pmax(I(age-bds[3])^2,0), data=bones.male)
model.quadratic.continuous.1derivative.female <- lm(spnbmd ~ poly(age,2, raw=T) + pmax(I(age-bds[2])^2,0) + pmax(I(age-bds[3])^2,0), data=bones.female)

** 2) State the number of coeficients fitted, including the intercept **

length(coef(model.quadratic.continuous.1derivative.female))

[1] 5

Taking the first derivate of a quadratic equation, means we lose the first terms for all but the first spline. So no we have 9-2 = 5 coeficients.

** 3) Show a plot with the fitted models superimposed **

plot (bone$age, bone$spnbmd,  col=ifelse(bone$gender=="male","blue","red"), pch = 8)
abline(v=bds, lwd=2, lty=3, col="blue")

curve(predict(model.quadratic.continuous.1derivative.male, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="blue", add=TRUE)

prediction from a rank-deficient fit may be misleading

curve(predict(model.quadratic.continuous.1derivative.female, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="red", add=TRUE)

prediction from a rank-deficient fit may be misleading

This model looks much nicer than the previous ones, it looks consistent accross the different groups. But maybe just a quadratic model is not detailed enough, so below there is a cubic continuous spline model.

Continuous piecewise cubic with continuous first and second derivatives (cubic spline)

** 1) Construct the model **

model.cubic.spline.male <- lm(spnbmd ~ poly(age,3, raw=T) + pmax(I(age-bds[2])^3,0) + 
            pmax(I(age-bds[3])^3,0), data=bones.male)
model.cubic.spline.female <- lm(spnbmd ~ poly(age,3, raw=T) + pmax(I(age-bds[2])^3,0) + 
            pmax(I(age-bds[3])^3,0), data=bones.female)

** 2) State the number of coeficients fitted, including the intercept **

length(coef(model.cubic.spline.male))

[1] 6

We fit cubic models for each of the groups, which has a total of 12 coeficients. But for all but the first group we only use the last cubic coefiecient. So there is a total of 4(first group) + 2 = 6.

** 3) Show a plot with the fitted models superimposed **

plot (bone$age, bone$spnbmd,  col=ifelse(bone$gender=="male","blue","red"), pch = 8)
abline(v=bds, lwd=2, lty=3, col="blue")

curve(predict(model.cubic.spline.male, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="blue", add=TRUE)
curve(predict(model.cubic.spline.female, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="red", add=TRUE)

This model is more complex than the previous models because it is cubic. But it looks like it represents the data well, and it is continuous at the knots.

Part 3

** Describe which model in Part 2 appears most suited to the data, based on visual inspection. Are there aspects of this fitted model that seem concerning? **

Both models for the Continuous piecewise quadratic with continuous first derivative(c) and the cubic spline(d) seem like good models. If a more detailed model is needed maybe d might be a better choice. If a simpler and easier to interpret model is needed, then c might be better.

Personally I think d is a better model. It seems to capture the data well, and create good models for both male and female data points. It also clearly shows the differences of models that occurr do to the different regions of age.

An aspect that is concerning about this model is the start and end ages, where it seems to quicly go up or down, which would not be realistic with the data. For example if we tried to predict a new female data point with age 8. It probably would be around -0.2, which is probalby very incorrect, because it cant be negative. For both male and female, predicitions of lower ages than 10, would probalby not yield good results.

** For this best approach, implement 5-fold cross-validation to select between using 2, 3, 4, 5, 6, or 7 knots, spaced out at uniform percentiles of the data. **

for ( j in 2:7)
{
  mse.male = 0
  mse.female = 0
  bones3 = bone
  bones3$agegrp <- quantileCut(bones3$age, j+1)
  bones.male = subset(bones3, gender == "male")
  bones.female = subset(bones3, gender == "female")
  labs <- levels(bones3$agegrp)
  bds <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
      upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) ))
  bds <- c(bds[,1], bds[nrow(bds),2])
  folds.female <- as.numeric(cut(1:nrow(bones.female), 5))
  folds.male <- as.numeric(cut(1:nrow(bones.male), 5))
  
  for (i in 1:5)
  {
    train.female <- bones.female[folds.female !=i,]
    val.female <- bones.female[folds.female == i,]
    train.male <- bones.male[folds.male !=i,]
    val.male <- bones.male[folds.male == i,]
    model.female=lm(spnbmd ~ bs(age ,knots = bds ),data=train.female) 
    model.male=lm(spnbmd ~ bs(age ,knots = bds ),data=train.male) 
    pred.female = predict(model.female)
    pred.male = predict(model.male)
    mse.female = mse.female + mean((pred.female - train.female$spnbmd)^2)
    mse.male = mse.male + mean((pred.male - train.male$spnbmd)^2)
  }
  cat(sprintf ("female - %s - s%f\n", j,mse.female))
  cat(sprintf ("male - %s - s%f\n", j,mse.male))
}

female - 2 - s0.006190
male - 2 - s0.008183
female - 3 - s0.005904
male - 3 - s0.008170
female - 4 - s0.005860
male - 4 - s0.008153
female - 5 - s0.005869
male - 5 - s0.008111
female - 6 - s0.005845
male - 6 - s0.008103
female - 7 - s0.005831
male - 7 - s0.008049

** What choice of knots works best? Refit a model on all of the data with this number of knots and plot the model. **

The best model used 4 knots (5 different groups), where the error was 0.005860. I chose this model because it has a low error rate, and it does not have too many knots, separating the data too much, and creating too many different groups.

Bellow is the plot of the final model, separated by gender.

bones2 = bone
bones2$agegrp <- quantileCut(bones2$age, 5)
labs <- levels(bones2$agegrp)
bds <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
      upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) ))
bds <- c(bds[,1], bds[nrow(bds),2])
bones.male = subset(bones2, gender == "male")
bones.female = subset(bones2, gender == "female")
final.male <- lm(spnbmd ~ bs(age ,knots = bds ),data=train.male) 
final.female <- lm(spnbmd ~ bs(age ,knots = bds ),data=train.female) 
plot (bone$age, bone$spnbmd,  col=ifelse(bone$gender=="male","blue","red"), pch = 8)
abline(v=bds, lwd=2, lty=3, col="blue")

curve(predict(final.male, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="blue", add=TRUE)

some 'x' values beyond boundary knots may cause ill-conditioned basesprediction from a rank-deficient fit may be misleading

curve(predict(final.female, data.frame(agegrp=cut(x, breaks=bds),
                              age=x)), lwd=2, col="red", add=TRUE)

some 'x' values beyond boundary knots may cause ill-conditioned basesprediction from a rank-deficient fit may be misleading

This final model with 4 knots captures the data very well. It shows the curve for women at a low age, were their spnbmd values are higher, then drops and stays relatively consistent after age 15. The peak for men is a bit later, at around age 15, and at around age 20 it stabilizes, being the same as the womans spnbmd.