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

## Women from Ireland and Poland
irlpolwomen <- subset(vitaminD, ( country == 4 | country == 6 ) & category == 2)
irlpolwomen$overweight <- factor(ifelse(irlpolwomen$bmi < 25,0,1), labels=c("Normal weight", "Over weight"))

## Remove the excessive labels for country.
irlpolwomen$country <- factor(irlpolwomen$country, labels = c("Ireland", "Poland"))

## Table of counts
counts <- table(irlpolwomen$country, irlpolwomen$overweight, dnn = c("Country", "Weight"))
(counts <- addmargins(counts)) ## parenthesis here causes the call to be printed

## Table of total Vitamin D in each group
vitD.total <- tapply(irlpolwomen$vitd, list(irlpolwomen$country, irlpolwomen$overweight), sum)
vitD.total <- addmargins(vitD.total)

## Using that mean = total/counts
vitD.total/counts