library(effects)
df <- read.csv("~/teach/data/442/braziltourism.csv")
row.names(df) <- paste("#",row.names(df))
df$sex <- factor(df$sex, levels=c(0,1),labels=c("M","F"))
df$favoursRoad <- factor(df$access.road, levels=c(0,1),labels=c("N","Y"))
df <- na.omit(df)
par(mfrow=c(2,2))
hist(trips,xlab="Number of trips")
prelim.plot(trips ~ sex + age + income, clear=FALSE)
sqrt.income <- sqrt(income)
clear(2,2)
pcts <- percents(table(sex,trips),denom=1,pretty=FALSE)
dotchart(t(pcts),xlab="Percents",ylab="Number of Trips", main="Trips by Sex")
prelim.plot(sqrt(trips) ~ age + income + logged.income,clear=FALSE)

par(mfrow=c(3,2))
par(cex=1.3,cex.main=1.1)
llFit <- glm( trips ~ sex + age + logged.income , family=poisson)
plot(llFit,cex=1.3,cex.main=1.2)
resids <- residuals(llFit)
plot(age,resids,main="Residuals vs. Age")
lines(supsmu(age, resids),col="red")
plot(logged.income,resids,main="Residuals vs. Logged Income")
lines(supsmu(logged.income, resids),col="red")

print(df["#412",c("trips","sex","age","income")])
detach()
attach(df[row.names(df)!="#412",])

llFit <- glm( trips ~ sex + age + logged.income , family=poisson,
   subset=(row.names(df) != "#412"))

par(mfrow=c(3,2))
par(cex=1.3,cex.main=1.1)
plot(llFit)
resids <- residuals(llFit)
plot(age,resids)
lines(supsmu(age, resids),col="red")
plot(logged.income,resids)
lines(supsmu(logged.income, resids),col="red")

print(summary(llFit))

llFitx <- glm( trips ~ sex *  (age + logged.income) , family=poisson)

par(mfrow=c(1,1))
plot(all.effects(llFitx),2)
print(cbind(anova(llFit,llFitx,test="Chi"),
      anova(llFit,llFitx,test="Cp")[,5,drop=FALSE]))