soluzione - aldo solari

Statistica Computazionale - Modulo R
Prova d’esame - 11 Febbraio 2014
COGNOME . . . . . . . . . . . . . . . . . . . . . . . . . . . NOME . . . . . . . . . . . . . . . . . . . . . MATRICOLA . . . . . . . . . . . . . . .
AVVERTENZE:
1. La prova dura 1 ora e 30 minuti. Non e’ consentito utilizzare il materiale delle lezioni.
2. Si risponda alle domande compilando negli spazi appositi le informazioni richieste. Il
testo va riconsegnato alla fine della prova.
3. Si riporti in un file di testo lo script R utilizzato per svolgere le analisi richieste (deve
riprodurre tutta l’analisi). Si converta questo file in formato .pdf, nominandolo con il
vostro numero di matricola, per poi salvarlo nell’apposita cartella CONSEGNA.
4. I datasets da utilizzare si trovano nella cartella TESTO.
1
nki70 data
Breve descrizione del dataset:
I dati nki70 sono stati ottenuti con la tecnologia microarray e rappresentano uno studio osservazionale
su 144 donne affette da tumore al seno alle quali e’ stata misurata l’espressione genica di 70 geni
ritenuti prognostici (si veda van de Vijver et al. 2002 NEJM). Sono inoltre presenti le seguenti
covariate:
• Diam: diametro del tumore;
• N: numero di linfonodi coinvolti;
• ER: status del recettore estrogeno;
• Grade: grado del tumore;
• Age: eta’ del paziente (in anni)
Obiettivo:
Si vogliono identificare i geni che risultano differenzialmente espressi nel confronto tra pazienti con
recettore estrogeno positivo e negativo.
1. Importare il file ALL.csv in R con il comando read.csv().
2. Si costruisca l’istogramma per la variabile Age.
1
3. Si ricodifichi la variabile Age nella variabile etaclassi con classi di intervallo (20,30], (30,45]
e (45,55].
4. Per il gene MELK, si calcolino i livelli di espressione genica medi per grado del tumore (comando
by()).
5. Si confronti il livello di espressione genica tra i gruppi ER==Positive e ER==Negative su tutti i
geni, utilizzando il test t di Welch, calcolando per ciascun gene il p-value (comando t.test()).
6. Si determini il numero di geni rifiutati ad α = 5%: con la procedura di Benjamini-Hochberg
per il controllo del False Discovery Rate.
7. Si identifichi il nome del gene piu’ significativo
>
>
>
>
nki70<-read.table(file="~/Documents/lectures/Statistica Computazionale/esame/nki70.csv"
hist(nki70$Age)
etaclassi<-cut(nki70$Age,breaks=c(20,30,45,55))
str(etaclassi)
Factor w/ 3 levels "(20,30]","(30,45]",..: 3 2 3 2 3 3 2 2 2 2 ...
> by(nki70$MELK,nki70$Grade,mean)
nki70$Grade: Intermediate
[1] -0.05704948
-----------------------------------------------------------nki70$Grade: Poorly diff
[1] 0.123264
-----------------------------------------------------------nki70$Grade: Well diff
[1] -0.1699878
>
>
>
>
>
+
+
>
risposte<-nki70[,-c(1:5)]
m<-ncol(risposte)
pvals<-vector("numeric", length=m)
names(pvals)<-names(risposte)
for (i in 1:m){
pvals[i]<-t.test(risposte[nki70$ER=="Positive",i],risposte[nki70$ER=="Negative",i])$p
}
sum(p.adjust(pvals,"BH")<=.05)
[1] 50
> names(pvals)[which.min(pvals)]
[1] "SCUBE2"
30
20
0
10
Frequency
40
50
Histogram of nki70$Age
25
30
35
40
45
50
55
nki70$Age
2
ASPIRIN data
Questi dati sono stati raccolti in uno studio di meta-analisi sull’efficacia dell’Aspirina (contro un
placebo) nel prevenire la morte dopo un infarto miocardico.
Studio
dp
tp
da
ta
Elwoodetal1974
67 624
49 615
Coronary1976
64 771
44 758
ElwoodSweetman1979 126 850 102 832
Breddinetal1979
38 309
32 317
Persantine1980
52 406
85 810
Aspirin1980
219 2257 246 2267
ISIS21988
1720 8600 1570 8587
• dp: numero di morti trattati con il placebo;
• tp: numero di soggetti trattati con il placebo;
• da: numero di morti trattati con l’Aspirina;
• ta: numero di soggetti trattati con l’Aspirina;
1. Si verifichi l’ipotesi sull’efficacia dell’Aspirina contro il Placebo, calcolando il p-value con il
test delle probabilita’ esatte di Fisher (comando fisher.test(), specificando come ipotesi
alternativa greater).
2. Si effettui una correzione per la molteplicita’ per i 7 test effettuati, calcolando i p-values
aggiustati con un metodo opportuno per il controllo del Familywise Error Rate.
>
>
>
>
>
>
>
+
+
+
+
+
>
require(HSAUR2)
data(aspirin)
M<-as.matrix(aspirin)
n<-nrow(M)
pvals<-vector("numeric", length=n)
names(pvals)<-rownames(aspirin)
for (i in 1:n){
m<-matrix(M[i,],nrow=2,byrow=T)
m[1,2]<-m[1,2]-m[1,1]
m[2,2]<-m[2,2]-m[2,1]
pvals[i]<-prop.test(m,alternative="greater")$p.value
}
pvals
\\cite{HSAuR:Elwoodetal1974}
0.057536698
\\cite{HSAuR:ElwoodSweetman1979}
0.071518114
\\cite{HSAuR:Persantine1980}
0.134056192
\\cite{HSAuR:ISIS21988}
0.002251286
\\cite{HSAuR:Coronary1976}
0.035545847
\\cite{HSAuR:Breddinetal1979}
0.227340887
\\cite{HSAuR:Aspirin1980}
0.889258463
> adj.pvals<-1 - (1 - pvals)^n
> adj.pvals
\\cite{HSAuR:Elwoodetal1974}
0.33953297
\\cite{HSAuR:ElwoodSweetman1979}
0.40514099
\\cite{HSAuR:Persantine1980}
0.63488643
\\cite{HSAuR:ISIS21988}
0.01565296
3
\\cite{HSAuR:Coronary1976}
0.22380451
\\cite{HSAuR:Breddinetal1979}
0.83559481
\\cite{HSAuR:Aspirin1980}
0.99999980
Domanda di Teoria
Si dimostri, esplicitando tutti i passaggi, che il metodo di Sidak controlla il familywise error rate.

Download Report