# Formation en biostatistiques (CIM 2007)
# Time-stamp: <2007-07-15 23:58:28 chl>

# définir PNG à TRUE permet de générer tous les graphiques lorsque l'on
# source le fichier 
PNG <- TRUE
set.seed(768)

##############################################################################
# test de Student -- illustration
x <- rnorm(20,mean=20,sd=2.5)
y <- rnorm(20,mean=22,sd=2.3)

if (PNG) png("ex1.png")
boxplot(x,y)
points(c(rep(1,20),rep(2,20)),c(x,y),col='gray50')
points(c(1,2),c(mean(x),mean(y)),pch='x',cex=2,col=c('blue','red'))
if (PNG) dev.off()

t.test(x,y,var.equal=T)  ## test clasique
t.test(x,y,var.equal=F)  ## test modifié de Welch

##############################################################################
# test de Student -- application
#
# Données : sommeil.txt
# 
# On a mesuré le temps de sommeil (en heure) de 10 sujets après 
# administration d'un certain médicament (numéroté 1 et 2). Chaque sujet 
# teste les 2 médicaments l'un après l'autre. On s'intèresse à l'effet du 
# type de médicament sur la durée de sommeil ; en particulier, on cherche à 
# savoir si les sujets dorment plus après avoir pris le médicament 2 que 
# lorsqu'ils prennent le médicament 1.
a <- read.table("sommeil.txt",header=T)
my.summary <- function(x,digits=2) {
  tmp <- c(summary(x),sd(x))
  names(tmp)[7] <- "Sd"
  return(round(tmp,digits))
}
attach(a)
tapply(temps,med,my.summary)

# visualisation des données
if (PNG) png("ex2.png")
par(mfrow=c(1,2))
boxplot(temps~med,col="gray",main="n=10",ylab="Temps (h)",xlab="M",names=c("1","2"))
plot(temps[med==1],temps[med==2],pch=16,main="n=10",xlim=c(2,12),ylim=c(2,12),xlab="M1",ylab="M2")
abline(0,1)
if (PNG) dev.off()

# test de Student pour échantillons appariés
t.test(temps[med==1],temps[med==2],var.equal=T, paired=T)
# à comparer au test erroné (sans prendre en compte l'appariement)
t.test(temps[med==1],temps[med==2],var.equal=T)

##############################################################################
# vérification des conditions d'application pour le test de Student
# code: C. Pallier (www.pallier.org)
nsuj <- 60
effectsize <- 0.5
subject <- gl(nsuj,1)
group <- gl(2,nsuj/2)
scores <- 1+rnorm(nsuj)+((0:(nsuj-1))%/%(nsuj/2))*effectsize
par(mfcol=c(1,3))
stripchart(scores~group,vert=T,method='jitter')
plot.design(scores~subject+group)
boxplot(scores~group)
by(scores,group,mean)
by(scores,group,summary)
t.test(scores~group)

# indépendance
if (PNG) png("ex3-1.png")
par(mfrow=c(1,2))
plot(scores[group==1]-mean(scores[group==1]))
abline(h=0,lty=2)
plot(scores[group==2]-mean(scores[group==2]))
abline(h=0,lty=2)
if (PNG) dev.off()

# homoscédasticité
if (PNG) png("ex3-2.png")
boxplot(scores~group)
if (PNG) dev.off()

# normalité
if (PNG) png("ex3-3.png")
par(mfrow=c(1,2))
qqnorm(scores[group==1])
qqline(scores[group==1])
qqnorm(scores[group==2])
qqline(scores[group==2])
if (PNG) dev.off()


##############################################################################
# anova 1 facteur -- illustration
# code: C. Pallier (www.pallier.org)
ns <- c(15,28,10,20,35)
n <- length(ns)
group <- factor(rep(1:n,ns),labels=paste("g",1:n,sep=""))
data <- rnorm(length(group),mean=100+(as.numeric(group)-2)^2)

if (PNG) png("ex4.png")
par(mfrow=c(1,3))
plot(data~group)
stripchart(data~group,method="jitter",vertical=T)
plot.design(data~group)
if (PNG) dev.off()

model2 <- aov(data~group)
model2
summary(model2)

if (PNG) png("ex5.png")
par(mfrow=c(2,2))
plot(model2)
if (PNG) dev.off()

##############################################################################
# anova 1 facteur -- application
#
# Données : pedago.txt
#
# Lors d'une expérimentation pédagogique, on désire comparer l'efficacité 
# de quatre méthodes d'enseignements.
# On dispose des notes obtenues à un examen par quatre groupes d'élèves 
# ayant chacun reçu un des 4 types d'enseignements.
a <- read.table('pedago.txt',header=T)[-1]
names(a) <- c("pedago","notes")
attach(a)

boxplot(notes~pedago)
#barplot(t(tapply(notes,pedago,mean)),ylim=c(0,20))
stripchart(notes~pedago,method='stack',vertical=T)
(moy <- as.numeric(tapply(notes,pedago,mean)))
points(1:4,moy,pch="X",col='blue',cex=2)
lines(1:4,moy,col='blue',lwd=2)

tapply(notes,pedago,summary)
tapply(notes,pedago,sd)

m <- aov(notes~pedago)
summary(m)

pairwise.t.test(notes,pedago)
TukeyHSD(m)
plot(TukeyHSD(m))

##############################################################################
# anova 2 facteurs -- illustration
x1 <- rnorm(25)+0.5*runif(25)
x2 <- rnorm(25)+0.25*runif(25)
x3 <- rnorm(25)+1.5*runif(25)
x4 <- rnorm(25)+2.5*runif(25)
x <- c(x1,x2,x3,x4)
a <- gl(2,50,100)
b <- gl(2,25,100)
levels(a) <- c("a1","a2")
levels(b) <- c("b1","b2")

plot(x~factor(a:b))
gmean <- tapply(x,list(a,b),mean)

model <- aov(x~a*b)
summary(model)

##############################################################################
# anova 2 facteurs -- application
#
# Données : tab13-2.dat
#
# Il s'agit de la même étude que celle présentée dans le chapitre 12. 
# Eysenck [Eyse1974] a en réalité mené une étude qui faisait varier aussi bien 
# l'âge que la condition de rétention. L'étude incluait 50 sujets dont l'âge se 
# situait entre 18 et 30 ans, ainsi que 50 sujets compris dans la tranche d'âge 
# 55-65 ans.
# Les données figurent dans le tableau ci-dessous. Ce ne sont pas les données 
# originales, mais elles ont été conçues pour avoir les mêmes moyennes et 
# écarts-types que les données enregistrées par Eysenck.
#
# Eysenck, M.W. (1974). Age differences in incidental learning. Developmental 
# Psychology, 10, 936-941.
a <- read.table("tab13-2.dat",header=F)
a$V1 <- factor(a$V1,labels=c("ages","jeunes"))
a$V2 <- factor(a$V2,labels=c("addition","rimes","adjectifs","images","controle"))
summary(a)

a.means <- with(a, tapply(V3,list(Age=V1,Groupe=V2),mean))
a.sd <- with(a, tapply(V3,list(Age=V1,Groupe=V2),sd))

my.barplot <- function(data, err, type, ...) {
  tmp <- barplot(data, ...)
  if (length(err) != length(data))
      stop("la longueur de 'data' et 'err' ne correspond pas")
  for (i in 1:length(err))
      arrows(tmp[i],data[i]-err[i],tmp[i],data[i]+err[i],code=3,angle=90,length=0.08)
}

my.barplot(a.means,a.sd,beside=T,ylab="Nb. mots",xlab="Condition",ylim=c(0,26),
                   names.arg=levels(a$V2),legend.text=levels(a$V1),cex.axis=.7,cex.names=.7,cex.lab=.7)

a.aov <- aov(V3 ~ V1 * V2,data=a)
summary(a.aov)
plot(a.aov)
model.tables(a.aov,se=T)

interaction.plot(a[,2],a[,1],a[,3])

##############################################################################
# anova "mesures répétées" -- illustration
suj <- as.factor(rep(seq(1,10),3))
c1 <- factor(c(rep(1,10),rep(2,10),rep(3,10)),labels=c("1jr","2jr","3jr"))
v <- rnorm(30,mean=20,sd=2)+as.numeric(c1)*rnorm(30,5)

boxplot(v~c1)
interaction.plot(c1,suj,v)

# modèle "classique" à mesures répétées
anova(lm(v~suj+c1))
summary(aov(v~suj+c1))
summary(aov(v~c1+Error(suj/c1)))

# approche modèle mixte
library(nlme)
summary(lme(v~c1,random=~1|suj))

##############################################################################
# régression linéaire simple -- illustration
a <- rnorm(100)
b <- 2*a+rnorm(100)

plot(b~a)

r <- lm(b~a)
anova(r)
summary(r)

abline(r)
plot(r)

##############################################################################
# régression linéaire simple -- application
#
# Données : moineaux.txt
#
# Les données recueillies indiquent la longueur des ailes de moineaux (en cm)
# en fonction de leur âge (en jours).
#
# Zar (1996). Biostatistical Analysis. Prentice Hall.
a <- read.table("moineaux.txt",header=T)
summary(a)
sd(a)

m <- lm(Y~X,data=a)
summary(m)

plot(a$Y~a$X,xlab="Age (jours)",ylab="Longueur (cm)")
abline(m)

plot(m)

anova(m)
library(car)  ## pour avoir un tableau d'ANOVA avec des SC de type III
Anova(m,type="III")

##############################################################################
# régression linéaire multiple -- application
#
# Données : rli_voix.csv
#
# Chez 101 patients déprimés hospitalisés en psychiatrie, on désire tester
# à l'entrée à l'hôpital, l'existence d'une association linéaire entre la 
# monotonie de la voix (mesurée par l'Af0s, paramètre quantitatif évaluant
# la variabilité de la hauteur de la voix) et l'intensité de la dépression, 
# mesurée par l'échelle de dépression de Hamilton (HDRS). La voix étant 
# sensiblement modifiée par les médicaments psychotropes, il est jugé nécessaire
# d'ajuster sur les variables binaires : antidépresseur tricyclique,
# antidépresseur sérotoninergique, neuroleptique, benzodiazépine (oui=1, non=0).
#
# Falissard (2005). Comprendre et utiliser les statistiques dans les sciences
# de la vie. Masson
a <- read.csv("rli_voix.csv",header=T,sep=";",dec=",")
str(a)
head(a)
attach(a)

par(mfrow=c(1,2))
hist(ham,prob=T)
lines(density(ham,na.rm=T),col='red')
hist(af0s,prob=T)
lines(density(af0s,na.rm=T),col='red')

pairs(a)

summary(lm(ham~af0s))
summary(lm(ham~af0s+tricyc+serot+neurol+bz))
summary(lm(ham~af0s+tricyc+serot+neurol+bz+tyrer))

##############################################################################
# sélection du meilleur modèle prédictif -- illustration
# exemple tiré de Azaïs & Bardet (2005). Le modèle linéaire par l'exemple.
# (Dunod)
# `seed` testé : 38, 1, 25
i <- 1:100
Z1 <- i
Z2 <- i^2
Z3 <- i^3
Z4 <- i^4
Z5 <- sqrt(i)
Z6 <- 1/i
Z7 <- log(i)
eps <- 20*rnorm(100,0)

# réalisations aléatoires tirées du vrai modèle tronquées sur une fenêtre
# [1,100]
Y <- 5-0.03*Z2+0.0002*Z3-3*Z7+eps

# le vrai modèle
j <- 1:150
YY <- 5-0.03*(j^2)+0.0002*(j^3)-3*log(j)

if (PNG) png("sel-model-4.png")
plot(j,YY,type="l",xlim=c(-10,160),ylim=c(-150,50))
points(i,Y)

# régression pas-à-pas descendante
# Note:
# les résultats diffèrent par rapport à l'ouvrage de Azaïs et Bardet
# puisque les observations sont issues de simulations numériques
# avec des "graines" variables
require(car)
y.lm7 <- lm(Y~Z1+Z2+Z3+Z4+Z5+Z6+Z7)
new <- data.frame(Z1=j,Z2=j^2,Z3=j^3,Z4=j^4,Z5=sqrt(j),Z6=1/j,Z7=log(j))
y.pred7 <- predict(y.lm7,new)
lines(j,y.pred7,lty="dashed")

# ex. de procédure de sélection "à la main" : tests de Fisher emboîtés
# Note:
# les valeurs de p reportées correspondent à 1 première simulation ;
# il ne faut pas en tenir compte et refaire la sélection des variables
# à chaque fois (mais on arrivera toujours à ne retenir que Z2 et Z3)
Anova(y.lm7)                    
Anova(lm(Y~Z1+Z2+Z3+Z4+Z5+Z7))  
Anova(lm(Y~Z1+Z2+Z3+Z4+Z5))     
Anova(lm(Y~Z1+Z2+Z3+Z4))        
Anova(lm(Y~Z2+Z3+Z4))           
Anova(lm(Y~Z2+Z3))              
# on retient finalement le modèle Y~Z2+Z3
y.lm2 <- lm(Y~Z2+Z3)
y.pred2 <- predict(y.lm2,new)
lines(j,y.pred2,lty="dashed",col="red")

# procédure de sélection à l'aide du critère Cp de Mallows
library(leaps)
Z <- matrix(c(Z1,Z2,Z3,Z4,Z5,Z6,Z7),ncol=7)
colnames(Z) <- c("Z1","Z2","Z3","Z4","Z5","Z6","Z7")
r <- leaps(Z,Y)
r$whi
r$Cp
t <- (r$Cp==min(r$Cp))
# variable(s) à retenir
colnames(Z)[r$whi[t]]          
y.pred.cp <- predict(lm(Y~Z1+Z4),new)
lines(j,y.pred.cp,lty="dashed",col="blue")

# procédure de sélection à l'aide des critères AIC/BIC
require(MASS)
ZZ <- as.data.frame(Z)
y.lm <- lm(Y~.,data=ZZ)
# pour utiliser BIC, on choisit k=log(n)
y.bic <- stepAIC(y.lm,k=log(100))
y.aic <- stepAIC(y.lm,k=2,trace=FALSE)
# variable(s) à retenir
y.bic                           
y.aic                           
y.pred.aic <- predict(lm(Y~Z1+Z2+Z3+Z5+Z7),new)
lines(j,y.pred.aic,lty="dashed",col="green")
y.pred.bic <- predict(lm(Y~Z2+Z3),new)
points(j,y.pred.bic,col="cyan")

legend("bottomleft",c("True Model","Complete Model","Stepwise desc.","Cp Mallows","AIC","BIC"),
       lty=c(1,rep(2,4)),col=c(rep("black",2),"red","blue","green","cyan"),bty="n")
if (PNG) dev.off()