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

women <- subset(vitaminD, category == 2)
women$country <- relevel(women$country, ref = "6") ## poland as reference level
women$sunexp <- relevel(women$sunexp, ref = "2") ## "sometimes sun" as reference level
       
model <- lm(log10(vitd) ~ bmi + country:log10(vitdintake) + country:sunexp + country,
            data = women)

## pol, den, fin, irl
countrypch <- c(21, 4, 2, 19)

## normalised dfbeta for all coefficients in the model
dfb.scaled <- dfbetas(model)

## and for the 2 terms of interest
finland.intake <- dfb.scaled[,which(model$coef==model$coef[["country2:log10(vitdintake)"]])]
finland.avoid <- dfb.scaled[,which(model$coef==model$coef[["country2:sunexp1"]])]

par(mfrow = c(1,2))
plot(finland.intake ~ log10(women$vitdintake),
     xlab = expression(paste(log[10], "(vitamin D intake)")),
     ylab = "Influence (Finland, vitamin D intake)",
     pch = countrypch[women$country])

plot(finland.avoid ~ women$bmi,
     xlab = "BMI",
     ylab = "Influence (Finland, avoid sun)",
     pch = countrypch[women$country])