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

## Irish and Polish women
irlpolwomen <- subset(vitaminD, ( country == 4 | country ==6 ) & category == 2)

## Linear model with bmi and country as covariates
model <- lm(vitd ~ bmi + as.factor(country),
            data = irlpolwomen)

intercept <- model$coef[1]
slope <- model$coef[2]
intercept.difference <- model$coef[3]

## The two parallel lines

curve(intercept + slope*x,
      min(irlpolwomen$bmi), max(irlpolwomen$bmi),
      xlab = "BMI",
      ylab = "Vitamin D",
      ylim = c(10,60)
      )

curve(intercept + intercept.difference + slope*x,
      min(irlpolwomen$bmi), max(irlpolwomen$bmi),
      add=TRUE,
      lty = "44"
      )