######################################################
#the following packages are necessary for specific robust analyses
#they must be first installed
# for the robust library, use Packages > Install package(s) 
# for the WRS library, go to https://r-forge.r-project.org/R/?group_id=468 to download it
# or type: install.packages("WRS", repos="http://R-Forge.R-project.org")

library(robust)     # necessary for M-estimation, robust regression and ANOVA
library(WRS)        # necessary only for Winsorized mean and variance

## For reproduceable research, we present here all commands that led
## to the figures of the tutorial.
## As users, the following steps are NOT useful: (3), (8), (9), (10), (14).


######################################################
## section 1: Central tendency and dispersion

#(1) reading the data on cigarette dependence and time to first cigarette
cig = read.csv("cigArt.csv", header=T)

#(2) creates figure 1
windows()
boxplot(cig$ttf)

#(3) creates figure 2
windows()
par(mar=c(5, 10, 4, 4)+.1)
x = c(2, 3, 4, 4, 5, 5, 6, 6, 6, 20)         # simulated data with one clear outlier
plot(2:20, c(rep(c(0, 8), 9), 3), type="n", xlab="", ylab="", yaxt="n")
axis(2, at=c(7, 5.5, 4, 2.5, 1),c("Mean","Median","20% trimmed\n mean","20% Winsorized\n mean","M-estimator"), las=1, cex.axis=1.4)
text(c(2, 3, 3.9, 4.1, 4.9, 5.1, 5.8, 6, 6.2, 20), y=-0, "|", cex=2)
lines(c(2, 2, 20, 20), c(6, 7, 7, 6)+.5, lwd=2)                             #mean
polygon(x=c(mean(x)-.1, mean(x), mean(x)+.1), y=c(7.3, 7, 7.3)+.5, col="black")
lines(c(2, 4.8, 4.8, 5.2, 20), c(3, 3, 4, 4, 3)+2, type="S", lwd=2)         #median
polygon(x=c(median(x)-.1, median(x), median(x)+.1), y=c(4.3, 4, 4.3)+2, col="black")
lines(c(2, 4, 4, 6, 20), c(0, 0, 1, 1, 0)+3.5, type="S", lwd=2)             #20% trimmed mean
polygon(x=c(mean(x, tr=.2)-.1, mean(x, tr=.2), mean(x, tr=.2)+.1), y=c(1.3, 1, 1.3)+3.5, col="black")
lines(c(2, seq(2, 4, .1), seq(6, 20, .1), 20),
      c(0, 1/(5-seq(2, 4, .1)), 1/abs(5-seq(6, 20, .1)), 0)+2, lwd=2)       #20% winsorized mean
polygon(x=c(mean(x, tr=.2)-.1, mean(x, tr=.2), mean(x, tr=.2)+.1), y=c(1.3, 1, 1.3)+2, col="black")
lines(c(2, seq(2, 2*4.594, .1), 20),
      c(0, ((((seq(2, 2*4.594, .1)-4.594)/4.594)^2)-1)^2, 0)+.5, lwd=2)     #M-estimator of central tendency
polygon(x=c(lmrob(x~1)$coeff-.1, lmrob(x~1)$coeff, lmrob(x~1)$coeff+.1), y=c(1.3, 1, 1.3)+.5, col="black")

#(4) estimation of central tendency and dispersion of time to first cigarette

ttfNoNA = na.omit(cig$ttf)             # Winsorized mean does not deal correctly with missing values
                                       # so it is necessary to take out all missing values

mean(ttfNoNA)                          # mean
median(ttfNoNA)                        # median
mean(ttfNoNA, trim=.2)                 # trimmed mean
win(ttfNoNA)                           # Winsorized mean
summary(lmRob(cig$ttf~1, data=cig))    # M-estimator of central tendency = intercept

sd(ttfNoNA)                            # standard deviation
IQR(ttfNoNA)                           # IQR
mad(ttfNoNA)                           # MAD
ttfNoNA20 = ttfNoNA[ttfNoNA>quantile(ttfNoNA, prob=.2) & ttfNoNA<quantile(ttfNoNA, prob=.8)]
sd(ttfNoNA20)                          # both lines code for 20% trimmed standard deviation
sqrt(winvar(ttfNoNA))                  # Winsorized standard deviation
summary(lmRob(cig$ttf~1, data=cig))    # M-estimator of dispersion = standard deviation of residuals


######################################################
## section 2: Linear regression

#(5) reading the data on creativeness and originality
creativ = read.csv("creatArt.csv",header=T)

#(6) OLS regression
creativ.lm = lm(creativity~original, data=creativ)
summary(creativ.lm)

#(7) robust regression (biweight)
creativ.lmRob = lmRob(creativity~original, data=creativ)
summary(creativ.lmRob)

#(8) creates figure 3
windows()
plot(creativ[, c("original", "creativity")], xlab="originality", ylab="creativeness", xlim=c(-2.5,2.5))
abline(creativ.lm$coeff, lwd=2)
abline(creativ.lmRob$coeff, lty=2, lwd=2)
segments(-2.57, 1.6, -2.57, 3)
segments(-2.52, 1.6, -2.52, 3.47, lty=2)
segments(-2.52, 1.6, -2.32, 1.6)
text(-2.32, 1.6, "observed value", pos=4)
segments(-2.52, 3, -2.32, 3)
text(-2.32, 3, "predicted value (OLS)", pos=4)
segments(-2.52, 3.47, -2.32, 3.47)
text(-2.32, 3.47, "predicted value (robust)", pos=4)
segments(1.545, 3.265, 1.545, 3.93)
segments(1.595, 3.265, 1.595, 3.695, lty=2)
legend("bottomright", legend=c("OLS    (b=0.23, p=.02)","robust (b=0.06, p=.69)"), lty=c(1, 2), lwd=2)

#(9) creates figure 4
windows()
ct=1.2; xmin=-1.5; xmax=1.5
x  =  seq(xmin, xmax, length=100)
plot(x, x^2, xlab = "residual value", ylab = "penalty function", type = "l", lwd=4)
lines(x, abs(x), lwd = 4, lty=2)
xtrunc = pmax (-1, pmin (x/ct, 1))
lines(x, 1*(3 * (xtrunc/ct)^2 - 3 * (xtrunc/ct)^4 + (xtrunc/ct)^6), lwd=4, lty=3)
legend ("bottomright", legend=c("OLS", "L1", "Biweight"), lty=1:3, lwd=4, cex=1)


#(10) creates figure 5
windows()
par(mfrow=c(1, 2))
plot(creativ.lm, which=1, caption="OLS regression", sub.caption=NA)
plot.default(creativ.lmRob$fitted, creativ.lmRob$residuals, xlab="Fitted Values", ylab="Residuals")
mtext("Robust regression", 3, .2)

#(10b) Both panels of figure 5 can be obtained more easily (but will look less similar) with:
windows()
plot(creativ.lm, which=1, caption="OLS regression", sub.caption=NA) #OLS regression
windows()
plot(creativ.lmRob,  which=5)

#(11) creates figure 6
windows()
plot(creativ.lm, which=4, caption=NA, sub.caption=NA)


######################################################
## section 3: ANOVA

#(12) reading the data on cigarette dependence and time to first cigarette
cig = read.csv("cigArt.csv", header=T)
levels(cig$income)<-c("low","a bit low","average", "a bit high","high")  # To have them in the correct order

#(14) creates figure 7
windows()
par(mfrow=c(2, 1))
boxplot(cig$ttf[cig$sexe=="male"]~cig$income[cig$sexe=="male"], ylab="time to first cigarette", las=3)
legend("topright", legend="male")
boxplot(cig$ttf[cig$sexe=="female"]~cig$income[cig$sexe=="female"], ylab="time to first cigarette", las=3)
legend("topright", legend="female")

#(14b) easiest way to obtain the 10 boxplots, but on a row
windows()
boxplot(cig$ttf~cig$income*cig$sexe, ylab="time to first cigarette")


#(16) the following two lines create a new dataset "cig2" with listwise deletion of missing values
# this deletion is necessary for the function named "anova".
cig2 = cig[,c("ttf", "income", "sexe")]
cig2 = na.omit(cig2)

#(17) classical anova
ttf.lm = lm(ttf~income*sexe, data=cig2)
anova(ttf.lm)

#(18) creates figure 8
windows()
plot(ttf.lm, which=1)

#(19) robust anova
ttf.lmRob = lmRob(ttf~income*sexe, data=cig2)
anova(ttf.lmRob)
## plot(ttf.lmRob, which=7)