source('http://192.38.117.59/~linearpredictors/datafiles/readVitaminD.R')

irlpolwomen <- subset(vitaminD, ( country == 4 | country ==6 ) & category == 2)  
irlpolwomen$country <- factor(irlpolwomen$country, labels = c("Ireland", "Poland"))
irlpolwomen$country <- relevel(irlpolwomen$country, "Poland")

modelNoInteraction <- lm(log10(vitd) ~ country + I(bmi-25), data = irlpolwomen)
modelInteraction <- lm(log10(vitd) ~ country*I(bmi-25), data = irlpolwomen)

## difference in lines for Ireland and Poland
diff <- function(x){
 modelInteraction$coef[2] + modelInteraction$coef[4]*(x-25)
}

## variance for each point of the estimated line
varDiff <- function(x){
  sigma <- summary(modelInteraction)$sigma
  cov <- summary(modelInteraction)$cov.unscaled * sigma^2
  var.a <- cov[2,2]
  var.b <- cov[4,4]
  cov.ab <- cov[2,4]
 (var.a + (x-25)^2*var.b + 2*(x-25)*cov.ab)
}

l <- range(irlpolwomen$bmi)[1]
r <- range(irlpolwomen$bmi)[2]
curve(diff, l, r, 
      ylim = c(-.5,.5),
      xlab = "BMI",
      ylab = "Country effect")

## pointvise confidence bands
curve(diff(x) - 1.96*sqrt(varDiff(x)), l, r, 
      add = TRUE, lty = "11")
curve(diff(x) + 1.96*sqrt(varDiff(x)), l, r, 
      add=TRUE, lty = "11")

curve(modelNoInteraction$coef[2] + 0*x, l, r,
      add = TRUE, lty = "32")