# Sample R code # We need to modify this for it to work. ex1028 = read.table("sleuth/ex1028.csv",header=T,sep=",") attach(ex1028) africa = rep("A",ncol(ex1028)) africa[west.africa==0] = "dry" africa[west.africa==1] = "wet" x = data.frame(year,el.nino,africa,storms,hurricanes,index=storm.index) rm(africa) detach() attach(x) # postscript("ex1028.ps",paper="letter",horizontal=T) par(mfrow=c(1,2)) levelsElNino = levels(as.factor(el.nino)) levelsAfrica = levels(as.factor(africa)) plot(year,storms,type="n") for(i in 1:length(levelsElNino)) { set = el.nino==levelsElNino[i] points(year[set],storms[set],pch=i) } legend(1950,19.5,levelsElNino,pch=1:length(levelsElNino)) plot(year,storms,type="n") for(i in 1:length(levelsAfrica)) { set = africa==levelsAfrica[i] points(year[set],storms[set],pch=i) } legend(1950,19.5,levelsAfrica,pch=1:length(levelsAfrica)) par(mfrow=c(1,2)) levelsElNino = levels(as.factor(el.nino)) levelsAfrica = levels(as.factor(africa)) plot(year,log(storms),type="n") for(i in 1:length(levelsElNino)) { set = el.nino==levelsElNino[i] points(year[set],log(storms[set]),pch=i) } legend(1950,19.5,levelsElNino,pch=1:length(levelsElNino)) plot(year,log(storms),type="n") for(i in 1:length(levelsAfrica)) { set = africa==levelsAfrica[i] points(year[set],log(storms[set]),pch=i) } legend(1950,19.5,levelsAfrica,pch=1:length(levelsAfrica)) par(mfrow=c(1,1)) hist(storms) par(mfrow=c(1,2)) fit1 = lm(storms ~ year + el.nino + africa) summary(fit1) plot(fitted(fit1),residuals(fit1)) abline(h=0,lty=2) lines(lowess(fitted(fit1),residuals(fit1))) fit2 = lm(log(storms) ~ year + el.nino + africa) summary(fit2) plot(fitted(fit2),residuals(fit2)) abline(h=0,lty=2) lines(lowess(fitted(fit2),residuals(fit2))) par(mfrow=c(1,2)) qqnorm(residuals(fit1)) qqnorm(residuals(fit2)) fit3 = lm(storms ~ el.nino * africa) summary(fit3) plot(fitted(fit3),residuals(fit3)) abline(h=0,lty=2) lines(lowess(fitted(fit3),residuals(fit3))) fit4 = lm(log(storms) ~ el.nino * africa) summary(fit4) plot(fitted(fit4),residuals(fit4)) abline(h=0,lty=2) lines(lowess(fitted(fit4),residuals(fit4))) fit5 = lm(storms ~ el.nino + africa) summary(fit5) plot(fitted(fit5),residuals(fit5)) abline(h=0,lty=2) lines(lowess(fitted(fit5),residuals(fit5))) fit6 = lm(log(storms) ~ el.nino + africa) summary(fit6) plot(fitted(fit6),residuals(fit6)) abline(h=0,lty=2) lines(lowess(fitted(fit6),residuals(fit6))) fit7 = lm(log(storms) ~ el.nino) summary(fit7) plot(fitted(fit7),residuals(fit7)) abline(h=0,lty=2) lines(lowess(fitted(fit7),residuals(fit7))) warm = as.factor(el.nino=="warm") fit8 = lm(log(storms) ~ warm) summary(fit8) plot(fitted(fit8),residuals(fit8)) abline(h=0,lty=2) fit9 = lm(log(storms) ~ africa) summary(fit9) plot(fitted(fit9),residuals(fit9)) abline(h=0,lty=2) dev.off()