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

## Women from Ireland and Poland
irlpolwomen <- subset(vitaminD, ( country == 4 | country == 6 ) & category == 2)

## Factor with levels corresponding to normal resp. overweight
## 1 = greater or equal to 25. Notice that NA's will still be NA using cut().
irlpolwomen$bmigroup <- cut(irlpolwomen$bmi, c(0, 25, Inf), right = F, labels = 0:1) 

model <- lm(log10(vitd) ~ bmigroup + country, data = irlpolwomen)

## Predicted Vitamin D values for the four groups.
levels <- expand.grid(bmigroup = factor(0:1), country = factor(c(4,6)))
predict(model, newdata = levels)