dat <- read.table("http://www.stat.umn.edu/~wuxxx725/data/anorexia.txt", 
                  header = TRUE)
# boxplot
boxplot(change ~ therapy, data = dat)

# ANOVA format: aov(numerical var ~ group var)
mod <- aov(change ~ therapy, data = dat)
summary(mod)


mod <- aov(change ~ therapy, data = dat)
TukeyHSD(mod, "therapy", conf.level = 0.95)