# Mark Liberman 12/30/2006
# This work is licensed under the 
#     Creative Commons Attribution-Noncommercial-Share Alike 2.5 License. 
# To view a copy of this license, 
# visit http://creativecommons.org/licenses/by-nc-sa/2.5/ 
# or send a letter to Creative Commons, 
# 543 Howard Street, 5th Floor, San Francisco, California, 94105, USA.

# Read in data from Fisher 2003
fish <- read.table("FisherStats")
colnames(fish) <- c("sex", "turns", "words", "time", "wpm", "id", "age", "edu")
#
# Set up some 
sex <- factor(fish[,"sex"])
words <- fish[,"words"]
wpm <- fish[,"wpm"]
age <- fish[,"age"]
edu <- fish[,"edu"]
#
males <- fish[,"sex"] == "m"
females <- fish[,"sex"] == "f"
# Fit a model and inspect the results...
M1.out <- summary(M1 <- lm(words ~ sex + age + edu + sex:age + sex:edu - 1))
#
# Boxplot of HS- vs. college-educated men & women more than 25 years old
#
# Regression aside, what are these effects like?
#
HS25 <- (fish[,"edu"] <= 12) & (fish[,"age"] >= 25)
COL25 <- (fish[,"edu"] >= 16) & (fish[,"age"] >= 25)
MHS25 <- HS25 & (fish[,"sex"] == "m")
FHS25 <- HS25 & (fish[,"sex"] == "f")
MCOL25 <- COL25 & (fish[,"sex"] == "m")
FCOL25 <- COL25 & (fish[,"sex"] == "f")

MHfish <- fish[MHS25,]
MCfish <- fish[MCOL25,]
FHfish <- fish[FHS25,]
FCfish <- fish[FCOL25,]

Pvals <- list(MHfish[,"words"], FHfish[,"words"], MCfish[,"words"], FCfish[,"words"])

png(filename="SexWords.png", width=400, height=400)

names(Pvals) <- c("M <=12e", "F <=12e", "M >=16e", "F >=16e")
boxplot(Pvals,
 col=c("blue","pink","blue","pink"),
 main="Data from Fisher 2003",
 ylab="Number of words spoken",
 notch=TRUE, boxwex=.5, outline=FALSE,
 ylim=c(0,1800))