Che materia stai cercando?

Statistical learning Appunti scolastici Premium

Appunti di statistical learning su: Boosting, adaboost, ensemble learning, algoritmi e inferenza, test multipli, p-hacking, FDP, FWER, FDR, Bonferroni, Holm, BH, BY, stima Storey, variable selection (single splits, multiple splits, stability selection)

Esame di Statistical learning docente Prof. A. Solari

Anteprima

ESTRATTO DOCUMENTO

Histogram of sim[1, ]

15

Frequency 10

5

0 0 20 40 60 80 100

X best

= 'hat beta')

hist(sim[2,],xlab 68

Histogram of sim[2, ]

30

25

20

Frequency 15

10

5

0 −3 −2 −1 0 1 2 3 4

hat beta

= 'P-value')

hist(sim[3,],xlab 69

Histogram of sim[3, ]

60

50

40

Frequency 30

20

10

0 0.00 0.02 0.04 0.06 0.08

P−value

p-value',mean(sim[3,]))

paste('media

## [1] "media p-value 0.0109661718430299"

coverage',mean(sim[4,]))

paste('media

## [1] "media coverage 0.01"

Lezione 7

Inferenza dopo selezione delle variabili

Sfortunatamente, l’utilizzo degli stessi dati per due proposte:

Ŝ.

1) Selezionare i predittori

2) Effettuare inferenza (intervalli di confidenza, p-values ecc. . . ) rispetto ai parametri dei predittori selezionati.

Introduce una distorsione sostanziale e invalida le proprietà inferenziali.

β

Inferenza su j

Il problema è l’inferenza dopo la selezione delle variabili, poichè se seleziono le vere variabili, allora tutto va

bene, invece se seleziono quelle sbagliate, posso avere dei risultati che in apparenza sono buoni ma in realtà

hanno coefficiente associato al regressore pari a 0 mentre il valore che viene stimato risulta essere diverso da 0.

β β

Se un modello è corretto lo stimatore dei è concentrato su e ha distribuzione normale. Posso calcolare

intervalli di confidenza ecc.. 70

− α

Abbiamo una certa probabilità di copertura 1 e un certo p-value.

Stimatore: −1

T

2

β̂ N , σ X)

(β (X )

px1 p px1 −1

T

2

β̂ N , σ X)

(β (X )

j j j

−1 −1

T T

X) j) X)

dove (X è il (j, elemento di (X Intervallo di confidenza:

j α

1− ·

± se(

β̂

β β̂ t

, ˆ )

] =

[β 2 j

j

j n−p

j

− α.

ha probabilità di copertura 1

P-value:

p per testare

j H j

H β p U nif orme(0,

: = 0 ha 1)

j j j

Simulazione n p

= 25, = 100 2

y X β ε , ε N , σ I

= + (0 )

nx1 nxp px1 nx1 nx1 nx1 nxn

X x

con Uniforme(0, 1)

nxp ij 2

β , σ

= 0 = 1

px1 px1

Il modello vero è: Y β N

= + (0, 1)

1 j

Quando vado a fare il modello fittato, non viene mai presa in considerazione l’intercetta pari a 1, quindi i

vengono presi casualmente. ⇒

Quando viene determinato a priori uniforme.

Quando seleziono non ho più un’uniforme.

Il modello fittato/stimato è: Y β β X N

= + + (0, 1)

j j

1

j

con fissato a priori o selezionato attraverso la forward stepwise selection, che è equivalente a:

| |

β̂ k

|ρ̂(X

j argmax , Y argmin ρ argmax

= )| = =

k k se(

β̂

ˆ )

k∈{1,...,p} k∈{1,...,p} k∈{1,...,p} k

71

n = 25

p = 100

set.seed(123)

X = ncol=p-1))

cbind(1,matrix(runif(n*(p-1)),

=

colnames(X) paste0("X",1:p)

beta = rep(0,p)

y = X %*% beta + rnorm(n)

#j è un qualsiasi valore basta che non sia pari a 1, dato che è l'intercetta.

#in if ho l'algoritmo di selezione forward.

sim <- function(X, beta, j=2, select = F){

n <- nrow(X) 72

p <- ncol(X)

y = X %*% beta + rnorm(n)

X = as.data.frame(X)

if (select==T) j = ) ) +1

which.max( apply(X[,-1],2,function(Xj) abs(cor(y,Xj))

fit = lm(y~X[,j],X)

hatbetaj = coef(summary(fit))[2,1]

pj = coef(summary(fit))[2,4]

confintj = confint(fit)[2,]

coveragej = ( confintj[1] <= beta[j] & confintj[2] >= beta[j] )

'hatbetaj'=hatbetaj , 'p-value j'=pj, 'copertura'=coveragej))

return( c('j'=j,

}

sim(X,beta,j=2)

## j hatbetaj p-value j copertura.2.5 %

## 2.0000000 0.4219562 0.5449314 1.0000000

B = 5

replicate(B, sim(X,beta,j=2))

## [,1] [,2] [,3] [,4] [,5]

## j 2.0000000 2.0000000 2.0000000 2.00000000 2.0000000

## hatbetaj -0.5292119 0.7010965 1.0298621 0.07346402 0.2126573

## p-value j 0.4205535 0.2933639 0.1772042 0.90343767 0.7889729

## copertura.2.5 % 1.0000000 1.0000000 1.0000000 1.00000000 1.0000000

select = T))

replicate(B, sim(X,beta,j=2,

## [,1] [,2] [,3] [,4]

## j.X35 35.000000000 70.000000000 6.000000e+01 95.000000000

## hatbetaj 1.433040679 1.486327839 2.202918e+00 1.712205958

## p-value j 0.009664025 0.005455845 2.787897e-04 0.003296043

## copertura.2.5 % 0.000000000 0.000000000 0.000000e+00 0.000000000

## [,5]

## j.X35 33.000000000

## hatbetaj -2.643281987

## p-value j 0.001275926

## copertura.2.5 % 0.000000000

#select=FALSE

B = 1000

j = 2

res = replicate(B, sim(X,beta,j=2))

# estimates

hatbetajs = res[2,] #prendo solo la seconda riga dell'output res, quindi hatbetaj

varhatbetajs = )[2,2]

solve(t(X[,c(1,j)])%*%X[,c(1,j)]

#calcolo la varianza tra i suoi hatbetaj

freq=F, 20, col="gray", border="white", las=1, main=paste("N(",beta[j], ",",

hist(hatbetajs, round(varh

sd=sqrt(varhatbetajs)), add=TRUE, lwd=2)

curve(dnorm(x,mean=beta[j], 73

N( 0 , 0.46 )

0.6

0.5

0.4

Density 0.3

0.2

0.1

0.0 −2 −1 0 1 2 3

hatbetajs

# p-values

pjs = res[3,]

freq=F, 20, col="gray", border="white", main="Uniform(0,1)", las=1)

hist(pjs, lwd=2)

abline(h=1, 74

Uniform(0,1)

1.2

1.0

0.8

Density 0.6

0.4

0.2

0.0 0.0 0.2 0.4 0.6 0.8 1.0

pjs

alpha = .05

typeIerrors = <= alpha)

mean(pjs

typeIerrors

## [1] 0.038

# coverage

coverages = res[4,]

mean(coverages)

## [1] 0.962

#select=TRUE

res = select=T))

replicate(B, sim(X,beta,j,

js = res[1,] xlab="selected j")

plot(table(js), 75

20

15

table(js) 10

5

0 2 7 13 20 27 34 41 48 55 62 69 76 83 90 97

selected j

# estimates

hatbetajs = res[2,]

freq=F, 20, col="gray", border="white", main="", las=1)

hist(hatbetajs,

varhatbetajs = )[2,2]

solve(t(X[,c(1,j)])%*%X[,c(1,j)]

sd=sqrt(varhatbetajs)), add=TRUE, lwd=2)

curve(dnorm(x,mean=beta[j], 76

0.4

0.3

Density 0.2

0.1

0.0 −3 −2 −1 0 1 2 3

hatbetajs

# p-values

pjs = res[3,]

freq=F, 20, col="gray", border="white", main="", las=1)

hist(pjs, lwd=2)

abline(h=1, 77

80

60

Density 40

20

0 0.00 0.01 0.02 0.03 0.04 0.05 0.06

pjs

typeIerrors = <= alpha)

mean(pjs

typeIerrors

## [1] 0.994

# coverage

coverages = res[4,]

mean(coverages)

## [1] 0.006

Lezione 8

Multiple testing

Testa o croce (introduzione al problema)

Supponiamo di lanciare una moneta per 5 volte di fila.

E

Sia l’evento “l’esito dei lanci è TTTTT oppure CCCCC”. E.

E

Se proviamo a lanciare una moneta per 5 volte di fila, otterremo oppure

E.

Provate a lanciare una moneta per 5 volte e ditemi se avete osservato

set.seed(1234567) = T,size = 5)

sample(c('T','C'),replace 78

## [1] "C" "C" "C" "T" "C"

“La moneta è bilanciata”.

Ipotesi nulla (H ):

0 “La moneta è truccata”.

Ipotesi alternativa (H): H E H

Rifiutiamo se otteniamo Rifiutare quando è vera.

Test statistico: Errore di primo tipo:

5

·

α = = 2 (0.5) = 6.25%

Probabilità errore di primo tipo: P(E)

Verifica di molte ipotesi

H “L’i-esima moneta è bilanciata”.

Ipotesi nulla i-esima :

i

m H , ..., H

Stiamo verificando ipotesi nulle m

1

m

S

F E E”

Sia = l’evento “almeno uno tra voi ha ottenuto

i

i=1

Se tutte le ipotesi nulle sono vere (ovvero le monete che avete utilizzato non sono truccate), allora

m m

− − −

) = 1 [1 = 1 (0.945)

P(F P(E)]

ovvero la probabilità di commettere almeno un errore di primo tipo.

test <- function(r){

res <- = T,size = r)

sample(c('T','C'),replace

e <- |

all(res==rep('T',r)) all(res==rep('C',r))

return(e)

}

test(5)

## [1] FALSE

mean(replicate(100,test(r=5)))

## [1] 0.06

Più test vengono effetuati e più avremo una probabilità di errore del primo tipo che tenderà ad 1. Diventa un

problema sopratuttto nel caso di molte variabili.

con il lancio di una moneta per 5 volte, su un campione di 20 persone che hanno testato, almeno 1 si prevede

che commetterà un errore di primo tipo, quindi almeno 1 persona avrà 5 teste o 5 croci.

79

Media e varianza del numero di errori

V = numero di errori del primo tipo.

m

P }

V = ha distribuzione Binomiale(m,

1{E P(E))

i

i=1 ·

m

) =

E(V P(E)

· · −

m

) = [1

Var(V P(E) P(E)]

pE <- 2*(0.5)^5 #Probabilità di $E$ che sia "TTTTT oppure CCCCC"

m <- 24 #numero di lanci

alpha <- 0.05 #livello di significatività

m*alpha #valore atteso

## [1] 1.2

m*alpha+(1-alpha) #varianza

## [1] 2.15 size=m, prob=pE)

qbinom(c(alpha/2,1-(alpha/2)),

## [1] 0 4

#intervallo di confidenza

Il problema H , ..., H .

Ipotesi nulle: m

1 ∈ T ⊆ {1,

H i ..., m}

è vera se

Ipotesi vere: i

T

L’insieme è incognito. ∈ R ⊆ {1,

H i ..., m}

è rifiutata se

Ipotesi rifiutate: i

Se l’esito dei test è basato sui è-values: R {i ≤

p α̃}

= : i

p α̃

dove è il p-value del test i-esimo e è un valore di soglia.

i R T ∩ R

Ottenere un insieme “grande” con “piccolo”.

Obiettivo:

∩ T

#(R ) = numero di errori del primo tipo. 80

Errore di primo tipo rifiuto di una ipotesi nulla = scoperta scientifica.

Scoperte scientifiche:

Errore di primo tipo vs errore di secondo tipo:

- Un errore di II tipo è un manacto passo in avanti.

- Un errore di I tipo è un passo nella direzione sbagliata.

Evitare gli errori di I tipo che sono più gravi degli errori di II tipo (ma non è sempre così,

Test multipli:

dipende dal contesto).

P-hacking

AIDSVAX - sperimentazione clinica

Il 24/2/2003 VaxGen ha annunciato i risultati della sperimentazione sull’efficacia di un vaccino per curare

l’AIDS:

Non si evidenziano differenze tra il vaccino e il gruppo placebo.

Osserviamo che i numeri sono piccoli, ma il risultato è statisticamente significativo. Si osservano incredibili

risultati.

AIDSvax <- 3330-191, 96, 191),nrow = 2,

matrix(c(1679-96,

dimnames = =

list(Group c("Placebo",

"Vaccine"),

Result = Infected", "Infected")))

c("Not

fisher.test(AIDSvax)$p.value

## [1] 1

White = 3003-179, 81, 179),2)

matrix(c(1508-81,

Black = 203-4, 9, 4),2)

matrix(c(111-9,

Asian = 53-4, 2, 4),2)

matrix(c(20-2,

Other = 71-6, 6, 6),2)

matrix(c(40-6, #non significativo

fisher.test(White)$p.value

## [1] 0.456536 81

#significativo

fisher.test(Black)$p.value

## [1] 0.01489004 #non significativo

fisher.test(Asian)$p.value

## [1] 0.662847 #non significativo

fisher.test(Other)$p.value

## [1] 0.3449772

fMRI (esperimento)

fMRI scan di un salmone morto (Bennet et al., 2009)

# numero di ipotesi (voxels)

m = 8064 #ipotesi nulle

# valore di soglia (euristico)

alphatilde = 0.001

# p-values

set.seed(123)

pvals = runif(m)

# ipotesi rifiutate (cioè tutte le ipotesi che sono maggiori del vaore soglia alpha tilde)

Rset = pvals <= alphatilde )

which(

# numero di rifiuti

R = length(Rset)

R

## [1] 7

#valore atteso di falsi positivi

m*alphatilde

## [1] 8.064

#il nostro errore è stato di 7, più o meno vicino al valore atteso

Che valore devo mettere alla soglia per avere 0.95 % di falsi positivi?

Abbiamo utilizzato un’uniforme perchè i p-value si distribuiscono più o meno in questo modo.

Tabella riassuntiva dell’esito di m test 82

X N (µ, 1)

i

H µ

: = 0

i i

H µ >

: 0

i i ∼

H X N

Quando è vera , (0, 1)

i i

set.seed(123)

m = 100

m0 = 50

setT = size=m0, replace=F)

sample(1:m,

stats <- rnorm(m)

stats[!setT] <- stats[!setT] + 1

setR = >

which(stats qnorm(.95))

rifiutate= 1:m %in% setR,

table(

vere = 1:m %in% setT)

## vere

## rifiutate FALSE TRUE

## FALSE 47 49

## TRUE 3 1

1

Errore I tipo: 47

Errore II tipo:

Guardo la diagonale principale. Commento in questo modo:

Le colonne corrispondono alle ipotesi nulle se sono vere oppure false (true or false). Le righe invece mi vanno

a dire quando sono rifiutate oppure no (TRUE=rifiutate, FALSE=non rifiutate).

Ho quindi la numerosità di errori di I tipo quando le ipotesi nulle sono vere e vengono rifiutate. La numerosità

degli errori di II tipo si osserva quando le ipotesi sono false e non vengono rifiutate.

set.seed(123)

m = 100

m0 = 50

setT = size=m0, replace=F)

sample(1:m,

stats <- rnorm(m)

stats[!setT] <- stats[!setT] + 1

pvals <- lower.tail = F)

pnorm(stats,

setR = <= 0.05)

which(pvals

rifiutate= 1:m %in% setR,

table(

vere = 1:m %in% setT) 83

## vere

## rifiutate FALSE TRUE

## FALSE 47 49

## TRUE 3 1

curve(dnorm,-2,2,ylim=c(0,0.4))

points(stats[1],0)

0.4

0.3

dnorm(x) 0.2

0.1

0.0 −2 −1 0 1 2

x

pvals[1]

## [1] 0.9541688

#guardando la cumulata

curve(pnorm,-2,2)

points(stats[1],0) 84

1.0

0.8

0.6

pnorm(x) 0.4

0.2

0.0 −2 −1 0 1 2

x

m

• è costante, ma potrebbe essere casuale quando seleziono un numero casuale di ipotesi in una variable

selection (è però un numero noto).

m m

• è un numero ignoto così come .

0 1

R

• è una variabile casuale quando dipende dai dati (in generale).

m R

• è casuale.

V U

• e sono variabili casuali.

FDP, FWER e FDR

False Discovery Proportion ( V se R > 0

R

F DP = altrimenti

0

Familywise Error Rate F W ER >

= 0)

P(V

False Discovery Rate F DR DP

= )

E(F

85

FWER e FDR (relazione)

FWER FDR ≥

F W ER > > DP F DR

= 0) = 0}) ) =

P(V E(1{V E(F

Se tutte le ipotesi nulle sono vere, FWER=FDR

F DP >

Se tutte le ipotesi nulle sono vere, = 0}

1{V

Nel caso di una sola ipotesi nullam FWER=FDR

F DP >

Nel caso di una sola ipotesi nulla, = 0}

1{V

sim <- function(m, m0, mu) {

stats <- rnorm(m)

stats[(m0+1):m] <- stats[(m0+1):m] + mu

setR <- >

which(stats qnorm(.95))

# numero di errori di I tipo

V <- setR <= m0 )

sum(

# family-wise error

FWE <- (V > 0)

#numero ipotesi rifiutate

R <- length(setR)

# false discovery proportion

FDP <- V/max(1,length(setR))

# true discovery proportion

TDP <- sum(setR>m0)/max(1,R)

FWE, R, FDP, TDP))

return(c(V,

}

res <- 50, 1))

replicate(1000, sim(100,

<- "FDP", "TDP")

row.names(res) c("V","I(V>0)","R",

plot(table(res[1,])) 86

250

200

])

table(res[1, 150

100

50

0 0 1 2 3 4 5 6 7 8

rowMeans(res[2:5,])

## I(V>0) R FDP TDP

## 0.9190000 15.5840000 0.1538233 0.8461767 N µ

Testo 100 ipotesi di cui 50 sono vere e le altre sono una realizzazzione di una (1, 1). Il valore di è dato

dalla media della normale.

FWER, su 1000, ottengo un errore del primo tipo 919 volte.

Il grafico mostra la distribuzione di V, cioè la frequenza di quante volte ho un errore del primo tipo.

Vedo cosa succede al variare segli input nella funzione.

sim <- function(m, m0, mu) {

stats <- rnorm(m)

stats[(m0+1):m] <- stats[(m0+1):m] + mu

setR <- >

which(stats qnorm(.95))

# numero di errori di I tipo

V <- setR <= m0 )

sum(

# family-wise error

FWE <- (V > 0)

#numero ipotesi rifiutate

R <- length(setR)

# false discovery proportion

FDP <- V/max(1,length(setR))

# true discovery proportion

TDP <- sum(setR>m0)/max(1,R)

FWE, R, FDP, TDP))

return(c(V, 87

}

res <- 175, 5))

replicate(1000, sim(200,

<- "FDP", "TDP")

row.names(res) c("V","I(V>0)","R",

plot(table(res[1,]))

120

])

table(res[1, 80

60

40

20

0 1 2 3 4 5 6 7 8 9 11 13 15 17 19

rowMeans(res[2:5,])

## I(V>0) R FDP TDP

## 1.0000000 33.8970000 0.2573181 0.7426819

sim <- function(m, m0, mu) {

stats <- rnorm(m)

stats[(m0+1):m] <- stats[(m0+1):m] + mu

setR <- >

which(stats qnorm(.95))

# numero di errori di I tipo

V <- setR <= m0 )

sum(

# family-wise error

FWE <- (V > 0)

#numero ipotesi rifiutate

R <- length(setR)

# false discovery proportion

FDP <- V/max(1,length(setR))

# true discovery proportion

TDP <- sum(setR>m0)/max(1,R)

FWE, R, FDP, TDP))

return(c(V,

}

res <- 3, 5))

replicate(1000, sim(200, 88

<- "FDP", "TDP")

row.names(res) c("V","I(V>0)","R",

plot(table(res[1,]))

800

600

])

table(res[1, 400

200

0 0 1 2

rowMeans(res[2:5,])

## I(V>0) R FDP TDP

## 1.480000e-01 1.970690e+02 7.728062e-04 9.992272e-01

Figure 1:

89

Lezione 9

P-values “probabilità di osservare un test statistico come estremo o più estremo della statistica test

Definizione:

osservata.”

Il p-value è una variabile casuale.

Un astatistica test standardizzata per ottenere una specifica distribuzione.

Distribuzione del p-value:

Un p-value valido ha distribuzione stocastica uguale o maggiore di un’uniforme

≤ ≤ ∀t ∈

t) t, (0, 1)

P(p

quando l’ipotesi nulla è vera, e stocasticamente più piccola di un’uniforme

≤ ∀t ∈

t) > t, (0, 1)

P(p

quando l’ipotesi nulla è falsa. ∼

p U nif orme(0,

• Se il test è esatto e la distribuzione della statistica test è continua, allora 1).

i

In pratica, i p-values sono spesso solo approssimati, come sono derivati attraverso argomenti asintotici

• o altre approssimazioni.

• I p-values asintotici possono essere poco accurati, specialmente per campioni di piccole dimensioni, e

che la loro relativa precisione diminuisce quando i p-values diventano piccoli.

Distribuzioni congiunte dei p-values: saldamente sotto controllo.

Distribuzioni marginali dei p-values:

potrebbe essere qualsiasi cosa. Tipicamente i p-values sono correlati.

Distribuzione congiunta: è il più grande problema pratico nei test multipli.

Distribuzione congiunta sconosciuta:

##FWER (procedure per il controllo)

Disuguaglianza di Boole

Per due eventi A e B ∪ − ∩

B) B)

= +

P(A P(A) P(B) P(A

quindi ∪ ≤

B) +

P(A P(A) P(B)

m E , ..., E

Per eventi m

1 m m

[ X

E )

P P(E

i i

i=1 i=1

abbiamo l’uguaglianza se gli eventi sono disgiunti. 90

sim <- function(n,p){

X <- cbind(1,matrix(runif(n*(p-1)),ncol=(p-1)))

XX <- as.data.frame(X)

beta <- c(1,0,0)

y <- X %*% beta + rnorm(n)

anova(lm(y~1),lm(y~.,XX))$"Pr(>F)"[2]

}

res <- replicate(1000,sim(n=10,p=3))

hist(res,freq=F)

abline(h=1,col=2) Histogram of res

1.2

1.0

0.8

Density 0.6

0.4

0.2

0.0 0.0 0.2 0.4 0.6 0.8 1.0

res

plot(ecdf(res))

curve(punif,0,1,add=T,col=2) 91

ecdf(res)

1.0

0.8

0.6

Fn(x) 0.4

0.2

0.0 0.0 0.2 0.4 0.6 0.8 1.0

x

Bonferroni α

α̃ =

Rifiutare tutte le ipotesi nulle con p-value inferiore a m

F W ER α

Garantisce senza nessuna assunzione.

Dimostrazione: [

α α α

X

≤ ≤ ≤ · ≤

≤ p m α

p P

P i i 0

m m m

i∈T i∈T

Tre disuguaglianze:

1) Usa la disuguaglianza di Boole. ≤ ≤ ∈ T

t) t i

2) Usa la proprietà dei p-values relativi a ipotesi nulle vere: per .

P(p

i

m m

3) Usa 0

P-values aggiustati: α.

Rifiutare tutte le ipotesi nulle con p-value aggiustato inferiore a

Il p-value aggiustato per la procedura Bonferroni è: ·

p̃ min(m p ,

= 1)

i i

92 α

α̃

Bonferroni è sempre valido e rifiuta tutte le ipotesi nulle con un p-value inferiore ad = .

m

p̃ corrisponde al p-value aggiustato e dalla sua forma risulta essere sempre superiore al p-value non aggiustato,

i ·

m p

poichè prende il minimo tra e 1.

i

Prostate cancer study dataset

n

• Studio effettuato su = 102 uomini, 50 controls contro 52 cases.

m

• = 6033 geni.

H i

• : attività del gene è lo stesso in casi e controlli.

i

• Dopo il codice per importare, possiamo osservare un boxplot dei livelli di espressione dell’i-esimo gene.

prostmat <- Learning/Dataset/prostmat.csv")

read.csv("F:/Statistical

dim(prostmat)

## [1] 101 102

X <- t(prostmat)

Holm = Bonferroni sequenziale

α

α̃

Iniziare con = m

Ripetere i passi 1) e 2) ≤ α̃

1) Rifiutare tutte le ipotesi con p-value

α

α̃

2) Aggiornare = m−r

r

dove è il numero di ipotesi nulle rifiutate finora.

HOLM rifiuta sempre di più che bonferroni, quindi è uniformemente più potente di Bonferroni.

93

###Holm (procedura step-down) 94

1) Ordinare i p-values dal più piccolo al più grande

≤ ≤ ≤

p p ... p

(1) (2) (m)

e le ipotesi corrispondenti ≤ ≤ ≤

H H ... H

(1) (2) (m)

{1,

j, ..., m}

2) Determinare il più piccolo indice tra 2, tale che

α

p >

(j) −

m j +1

H i < j

3) Rifiutare tutte le ipotesi tali che

(i)

###Hochberg (procedura step-up) 95

1) Ordinare i p-values dal più piccolo al più grande

≤ ≤ ≤

p p ... p

(1) (2) (m)

{1,

j, ..., m}

2) Determinare il più grande indice tra 2, tale che

α

p (j) −

m j +1

H i j

3) Rifiutare tutte le ipotesi tali che

(i)

La procedura di Hochberg rifiuta come o più di Holm, ma bisogna assumere che i p-values

{p ∈ T }

, i soddisfano la diseguaglianza di Simes (dipendenza positiva).

i vado verso destra: se sto sotto, continuo. Nel momento in cui vado sopra, mi fermo e rifiuto

Step-down

tutto quello che ho prima.

Nella invece vado da destra a sinistra e se sto sopra, allora continuo. Quando rifiuto, andrò sotto,

step-up

allora mi fermo e rifiuto tutto quello che ho a sinistra. Quindi nella secondo rifiuterò di più.

96

Relazioni logiche

ANOVA ad una via: 2 2 2

∼ ∼ ∼

Y N , σ Y N , σ Y N , σ

(µ ), (µ ), (µ )

1 1 2 2 3 3

Confronti a coppie: H µ µ H µ µ H µ µ

: = : = : =

12 1 2 13 1 3 23 2 3

Relazioni logiche tra ipotesi nulle:

H H H

Ad esempio, se è falsa, allora e non possono essere contemporaneamente vere.

12 13 23

Shaffer = Holm con relazioni logiche

α

α̃

Iniziare con = m

Ripetere i passi 1) e 2) ≤ α̃

1) Rifiutare tutte le ipotesi con p-value

αs

α̃

2) Aggiornare =

s

dove è il numero massimo di ipotesi nulle che possono essere contemporaneamente vere, assumendo che le

ipotesi rifiutate finora siano ipotesi false.

Shaffer è uniformemente più potente di Holm.

Confronti a coppie:

H µ µ p µ µ p µ µ p

: = = 0.01H : = = 0.54H : = = 0.04

12 1 2 12 13 1 3 13 23 2 3 23

α̃

Valore critico

m α

1) = 3, = 5% 0.05 ≈

α̃ H

= 0.0167 e rifiuto 12

3 →

H H H s

2) Se è falsa, allora al massimo una tra e può essere vera = 1

12 13 23

0.05 H

α̃ e rifiuto

= 23

1 →

H H H s

3) Se e sono false, allora può essere vera = 1

12 23 13 0.05

α̃ = ma non rifiuto alcuna ipotesi

1

Ipotesi nulle dei confronti a coppie

H µ µ H µ µ H µ µ

: = : = : =

12 1 2 13 1 3 23 2 3

Ipotesi nulla globale ∩ ∩

H H H H µ µ µ

= : = =

123 12 13 23 1 2 3

Relazioni logiche

H H , H H

Se è falsa, al massimo una tra e può essere vera.

123 12 13 23

Vediamo graficamente ANOVA ad una via con 3 gruppi

97

Vado a fare test d’ipotesi sulla congiunta, ad esempio il test se rifiuto, allora vado ad analizzare le singole

F,

ipotesi e almeno una sarà vera.

Osserviamo un altro dataset:

prostz <- Learning/Dataset/prostz.txt", quote="\"", comment.char="")

read.table("F:/Statistical

pvals <- = F)*2

pnorm(abs(prostz$V1),lower.tail

# Histogram col="lightgray", main="", xlab="p-value", freq=F)

hist(pvals,20,

abline(h=1,lty=2) 98

1.5

1.0

Density 0.5

0.0 0.0 0.2 0.4 0.6 0.8 1.0

p−value

#grafico che ordina i p-values

Simulazione (FWER control)

m = length(pvals)

alpha = 0.05

<= alpha)

sum(pvals

## [1] 478

# Bonferroni <= alpha)

sum(p.adjust(pvals,"bonferroni")

## [1] 3

# Holm and Hochberg <= alpha)

sum(p.adjust(pvals,"holm")

## [1] 3 <= alpha)

sum(p.adjust(pvals,"hochberg")

## [1] 3 xlim=c(1,20), ylim = xlab="i", ylab="p-value")

plot(1:m, sort(pvals), c(0,.0001),

alpha/(m- 1:m +1))

lines(1:m, 99

8e−05

p−value 4e−05

0e+00 5 10 15 20

i

FDR (procedure per il controllo)

Benjamini & Hochberg

1) Ordinare i p-values dal più piccolo al più grande

≤ ≤ ≤

p p ... p

(1) (2) (m)

{1,

j, ..., m}

2) Determinare il più grande indice tra 2, che soddisfa

p (j) m

Diversamente da bonferroni ho un valore critico differente.

≤ α̃ p

3) Rifiutare tutte le ipotesi con p-value = (j)

• La procedura assume che i p-values sono indipendenti o più generalmente con dipendenza positiva

BH

- positive regression dipendence on the subset of p-values of true null hypotheses (PRDS).

m

F DR α

• Se è vera l’assunzione PRDS, garantisce

BH 0

m

• è una procedura

BH step-up. 100

101

#wilcoxon.test

#pvals2 <- apply(X,2, function(x) wilcox.test(x=x[1:50],y=x[-c(1:50)])$p.value)

#hist(pvals2,20)

#sum(pvals2 <= alpha)

#sum(p.adjust(pvals2,"BH") <= alpha)

In questo caso rifiutiamo di più, quindi controllo la percentuale di falsi positivi in valore atteso. Mi aspetterò

·

circa 3 falsi positivi, dato dal seguente calcolo: 61 0.05.

Benjamini & Yekutieli

1) E’ una variante di valida per qualsiasi forma di dipendenza tra i p-values.

BH iα , i ..., m

2) Rende più stringenti i valori critici = 1, moltiplicandoli per

m 1

1 m→∞

=

m 1

P log(m)

i=1 i

102

#wilcoxon.test

#pvals2 <- apply(X,2, function(x) wilcox.test(x=x[1:50],y=x[-c(1:50)])$p.value)

#hist(pvals2,20)

#sum(pvals2 <= alpha)

#sum(p.adjust(pvals2,"BY") <= alpha)

Qui non funziona molto bene, poichè il metodo è troppo stringente.

Simulazione (FDR control)

<= alpha)

sum(p.adjust(pvals,"BH")

## [1] 21 xlim=c(1,50), ylim = xlab="i", ylab="p-value")

plot(1:m, sort(pvals), c(0,.002),

alpha/(m- 1:m +1))

lines(1:m, alpha*(1:m)/m, col=2)

lines(1:m,

0.0020

p−value 0.0010

0.0000 0 10 20 30 40 50

i

#Benjamini & Yekutieli <= alpha)

sum(p.adjust(pvals,"BY")

## [1] 2 103

Lezione 10

FDP (stima e confidenza)

Stima FDP di Storey

Impostazione rifiutata:

R {H ≤

p t} R

Supponiamo di rifiutare = : con = #R.

i i 104

Intuizione: m t

Dall’uniformità dei p-values sotto FDP nulla 0

R

m

Stima di 0 ≥ −

> λ}) m λ),

Idea: (1 così

E(#{p i 0 > λ}

#{p + 1

i

m̂ =

0 − λ

1

Stima di FDP (q-value) m̂ t t > λ}

#{p + 1

i

0

ˆ

F DP (t) = = − ≤

R λ t}

1 #{p

i

Figure 2:

Storey e Benjamini & Hochberg

Relazioni vicine:

Una via alternativa di costruzione di impostazione rifiutata:

BH

1) Stimare = 1 invece di una stima di Storey.

0 ˆ ≤

t F DP α.

2) Prendere il valore più grande come (t)

Un look alternativo di Storey:

Metodo Storey = controllo adattato FDR.

Un look alternativo di BH:

Stima conservativa di FDP.

Storey e dipendenza

Metodo delle stime dei momenti:

Solo dipendenzain media Non affetta da una struttura correlata.

Però:

La variabilità della stima può essere grande se i p-values sono correlati.

105

Standard errors non disponibili:

Disponibili solo per p-values indipendenti.

Simulazione (Prostate data)

#dataset già importato in precedenza

# Sorted p-values plot type="l", xlab="sorted p-value")

plot(1:m, sort(pvals),

lines(1:m,(1:m)/m,lty=2)

1.0

0.8

sort(pvals) 0.6

0.4

0.2

0.0 0 1000 2000 3000 4000 5000 6000

sorted p−value

Stima FDP: Storey

(m = length(pvals))

## [1] 6033

# Stima del numero di ipotesi vere

lambda = 0.5

(hatm0 <- (sum(pvals > lambda) + 1)/(1-lambda))

## [1] 5582

# Stima della proporzione delle ipotesi vere

(hatpi0 = hatm0/m)

## [1] 0.9252445 106

alpha = 0.001

(R = <= alpha))

sum(pvals

## [1] 60

(hatfdp = (hatm0*alpha)/R)

## [1] 0.09303333 ≤

α 5%

Quali hanno FDP

alphas = seq(0,0.05,length=5000)

Rs = function(alpha) <= alpha))

sapply(alphas, sum(pvals

plot(alphas,Rs,type="l")

400

300

Rs 200

100

0 0.00 0.01 0.02 0.03 0.04 0.05

alphas

hatfdps = (hatm0*alphas)/Rs ylim=c(0,1))

plot(alphas,hatfdps,type="l",

abline(h=.05,lty=2) 107

1.0

0.8

0.6

hatfdps 0.4

0.2

0.0 0.00 0.01 0.02 0.03 0.04 0.05

alphas

(alphatilde = alphas[which.min(hatfdps<=0.05)])

## [1] 0.000190038

(Rs[which.min(hatfdps<=0.05)])

## [1] 21 m

Vediamo molte stime 0

lambdas = seq(0.05,0.5,length=1000)

hatm0s <- function(lambda) (sum(pvals > lambda) + 1)/(1-lambda) )

sapply(lambdas,

hatm0s, type="l", ylim=c(0,m))

plot(lambdas,

abline(h=m,lty=2) 108

5000

hatm0s 3000

1000

0 0.1 0.2 0.3 0.4 0.5

lambdas

Test dello scimpanzè (Ignorance Project)

E’ stato svolto un questionario da parte di 88 studenti riguardo ad argomenti di cultura generale, ma

abbastanza specifici.

Domande che ci poniamo prima di effettuare l’analisi

Numero risposte al questionario = 88.

Numero domande in ogni questionario = 12.

Ipotesi:

Non casuale: 1 1

6

H π H π

Nulla: : = Alternativa: : =

j i

i i 3 3

Meglio del caso: 1 1

H π H π >

Nulla: : Alternativa: :

i i j i

3 3

Peggio del caso: 1 1

H π H π <

Nulla: : Alternativa: :

i i j i

3 3

Confronti a coppie:

1) Confronto maschi contro femmine (12 test).

2) Oroscopo si contro no (12 test).

3) Amici su Facebook (1 test di correlazione).

library(readr)

Rosling <- Learning/Dataset/Rosling.csv")

read_csv("F:/Statistical 109

## Parsed with column specification:

## cols(

## .default = col_character(),

## `4. L' aspettativa di vita alla nascita nel Regno Unito è di 81 anni. Qual è l'aspettativa di vita

## `11. Nel 1965, il numero di bambini nati per ogni donna nel mondo, in media, era cinque. Quanto val

## `H. Scegli un numero tra 0 e 9` = col_integer(),

## `I. Qual è il tuo peso? (solo il numero e.g. 70.4)` = col_number()

## )

## See spec(...) for full column specifications.

quest <- as.data.frame(Rosling[-52,2:13])

<- sep='')

names(quest)[1:12] paste('Q',1:12,

p = ncol(quest)

head(quest)

## Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11

## 1 Quattro milardi C 60% 60 5 anni 50% 10% 50% B Raddoppiata 2.5

## 2 Quattro milardi B 40% 60 3 anni 20% 1% 50% A Rimasta la stessa 2.5

## 3 Quattro milardi B 60% 60 3 anni 20% 5% 50% C Raddoppiata 2.5

## 4 Quattro milardi C 40% 50 5 anni 20% 5% 50% A Rimasta la stessa 3.5

## 5 Quattro milardi B 60% 60 5 anni 50% 5% 50% A Rimasta la stessa 3.5

## 6 Quattro milardi B 80% 60 5 anni 50% 5% 50% B Dimezzata 2.5

## Q12

## 1 Rimasto lo stesso

## 2 Dimezzato

## 3 Raddoppiato

## 4 Raddoppiato

## 5 Rimasto lo stesso

## 6 Dimezzato

X <- as.data.frame(Rosling[-52,c(14,20,23)])

<-

names(X) paste('X',1:3,sep='')

head(X)

## X1 X2 X3

## 1 Femmina No 800

## 2 Maschio No 300

## 3 Femmina Si 1300

## 4 Maschio Si 300

## 5 Maschio No 550

## 6 Maschio No 650

#Risposte corrette

rights <- miliardi","A","80%",70,"7 anni","80%","1%","80%","B","Dimezzata",2.5,"Dimezzato")

c("Due

correct <- function(i) quest[i]==rights[i] )

sapply(1:p,

Scores

score <- rowSums(correct)

plot(table(score)) 110

30

25

table(score) 20

15

10

5

0 0 1 2 3 4 5 6 7 8

score

#Posso consiederare un'eliminazione degli outlier nel seguente modo

#which(score==12) fatto all'inizio, quando ho importato il dataset (tolto l'osservazione 52)

qscore = colMeans(correct)

barplot(rbind(qscore,1-qscore),

names.arg = names(quest))

abline(h=1/3) 111

1.0

0.8

0.6

0.4

0.2

0.0 Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q12

Intervalli di confidenza

alpha = 0.05/37 #diviso il numero di test che effettuiamo (penso sia sbagliato, perchè ci sono 25 test c

alt <- 'less' #saperne meno degli chimpanzee

(cis <- function(i) alternative=alt,

sapply(1:p, binom.test(sum(correct[,i]),length(correct[,i]),p=1/3,

## [,1] [,2] [,3] [,4] [,5] [,6] [,7]

## [1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

## [2,] 0.2541249 0.2266291 0.3412211 0.2720533 0.4072026 0.1376718 0.5767683

## [,8] [,9] [,10] [,11] [,12]

## [1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

## [2,] 0.3909613 0.3579807 0.3070986 0.6276499 0.3827801

Bonferroni aggiustato:

# Bonferroni adjusted

(cis.bonf = function(i) alternative=a

sapply(1:p, binom.test(sum(correct[,i]),length(correct[,i]),p=1/3,

## [,1] [,2] [,3] [,4] [,5] [,6] [,7]

## [1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

## [2,] 0.2937949 0.2652976 0.3829495 0.3122787 0.4495237 0.1715805 0.6172811

## [,8] [,9] [,10] [,11] [,12]

## [1,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000

## [2,] 0.4332073 0.3999331 0.3482065 0.6666943 0.4249712 1

Vediamo graficamente gli intervalli di confidenza e l’abline che determina il caso 3

cis[1,], ylim=c(0,1), pch="-", xlab="Question number", ylab="Probability of correct answer")

plot(1:p, cis[2,], pch="-")

points(1:p, 112

col=2)

abline(h=1/3,

for (i in 1:p) cis[1,i], y1 = cis[2,i])

segments(x0=i,

1.0

answer 0.8

correct −

0.6 −

of 0.4 − − −

Probability − −

− −

0.2 −

0.0 − − − − − − − − − − − −

2 4 6 8 10 12

Question number

aggiustiamo ora il p.value

alpha = 0.05

alt <- 'less' #saperne meno degli chimpanzee

ps = function(i) alternative=alt, con

sapply(1:p, binom.test(sum(correct[,i]),length(correct[,i]),p=1/3,

#momento in cui aggiustiamo il p-value

round(p.adjust(ps,'holm'),4)

## [1] 0.0000 0.0000 0.0171 0.0000 0.3095 0.0000 1.0000 0.2002 0.0467 0.0012

## [11] 1.0000 0.1658

Maschi contro Femmine

Osservo che considero in questo caso prop.test e non più binom.test. E’ evidente che sto applicando a qualcosa

di più pratico. Non vado più ad effettuare test che riguardano il caso.

alt = "two.sided"

ps2 <- function(j)

sapply(1:p, ) ,

prop.test(x=c(sum(correct[X$X1=="Maschio",j]),sum(correct[X$X1!="Maschio",j])

n=c(sum(X$X1=="Maschio"),sum(X$X3!="Maschio")), alternative=alt)$p.value

)

## Warning in prop.test(x = c(sum(correct[X$X1 == "Maschio", j]),

## sum(correct[X$X1 != : Chi-squared approximation may be incorrect

round(ps2,4)

## [1] 0.2406 0.0032 0.0004 0.0286 0.0013 0.1376 0.0000 0.0964 0.7310 0.0018

113

## [11] 0.0806 0.0028

qscoreA = colMeans(correct[X$X1=="Maschio",])

qscoreB = colMeans(correct[X$X1!="Maschio",])

beside=T, names.arg =

barplot(rbind(qscoreA,qscoreB), names(quest))

0.6

0.5

0.4

0.3

0.2

0.1

0.0 Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q12

oroscopo

alt = "two.sided"

ps3 <- function(j)

sapply(1:p, ) ,

prop.test(x=c(sum(correct[X$X2=="Si",j]),sum(correct[X$X2!="Si",j])

n=c(sum(X$X2=="Si"),sum(X$X2!="Si")), alternative=alt)$p.value

)

## Warning in prop.test(x = c(sum(correct[X$X2 == "Si", j]), sum(correct[X

## $X2 != : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X2 == "Si", j]), sum(correct[X

## $X2 != : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X2 == "Si", j]), sum(correct[X

## $X2 != : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X2 == "Si", j]), sum(correct[X

## $X2 != : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X2 == "Si", j]), sum(correct[X

## $X2 != : Chi-squared approximation may be incorrect

114

round(ps3,4)

## [1] 0.6067 0.8508 0.0510 1.0000 0.9586 0.4657 0.0300 0.8158 0.9399 0.6290

## [11] 0.1570 0.4683

qscoreA = colMeans(correct[X$X2=="Si",])

qscoreB = colMeans(correct[X$X2!="Si",])

beside=T, names.arg =

barplot(rbind(qscoreA,qscoreB), names(quest))

0.6

0.5

0.4

0.3

0.2

0.1

0.0 Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q12

Facebook

alt = "two.sided"

ps4 <- function(j)

sapply(1:p, ) ,

prop.test(x=c(sum(correct[X$X3>100,j]),sum(correct[X$X3<=100,j])

n=c(sum(X$X3>100),sum(X$X3<=100)), alternative=alt)$p.value

)

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

115

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

## Warning in prop.test(x = c(sum(correct[X$X3 > 100, j]), sum(correct[X$X3

## <= : Chi-squared approximation may be incorrect

round(ps4,4)

## [1] 0.4265 0.7165 0.6179 0.1959 0.4762 1.0000 0.6163 1.0000 0.7764 0.6439

## [11] 0.3146 0.2718

qscoreA = colMeans(correct[X$X3>100,])

qscoreB = colMeans(correct[X$X3<=100,])

beside=T, names.arg =

barplot(rbind(qscoreA,qscoreB), names(quest))

0.5

0.4

0.3

0.2

0.1

0.0 Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q12

116

round(p.adjust(c(ps,ps2,ps3,ps4),'holm'),4)

## [1] 0.0000 0.0000 0.0955 0.0002 1.0000 0.0000 1.0000 1.0000 0.2802 0.0063

## [11] 1.0000 1.0000 1.0000 0.1186 0.0161 1.0000 0.0539 1.0000 0.0000 1.0000

## [21] 1.0000 0.0710 1.0000 0.1073 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

## [31] 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

## [41] 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

Lezione 11

Sample-Splitting Inference

In questa lezione ci concentriamo sull’utilizzo del sample-splitting come metodo per effettuare inferenza per

modelli di grandi dimensioni (high-dimensional).

- Questi metodi tendono ad essere computazionalmente intensi, in quanto possono comportare stime a modelli

high-dimensional centinaia o migliaia di volte.

- Questo non è necessariamnete proibitivo dal punto di vista dell’analisi di un singolo dataset, ma presenta

barriere, ostacoli per studi di simulazione completi.

Model

Linear model: y X β ε

= +

nx1 nxp px1 nx1

2

Error ε N σ I

(0, )

n

M atrice di disegno X ixed or ranodm)

(f 0

P arametri β , ..., β

= (β )

px1 p

1

≤ −

Basso n) o alto > n) impostazione della dimensione

(low) (p (high) (p

Inferenza (inferenze)

H β

P-values per : = 0

j j β

Intervalli di confidenza per j

{j 6

S β j ..., p}

= : = 0, = 1,

j

Simulated setting and example 2

σ

Considero la numerosità delle osservazioni pari a 100, con il numero di predittori pari a 60 e = 1.

6

β

p = 6 corrisponde al predittore di tipo con = 0. Sono i predittori con coefficienti associati non nulli e

A j

1

risultano essere quelli fortemente diversi da 0.

Tra questi 6 coefficienti possiamo distinguere:

±

β

- 2 forti predittori con = 1.

j ±

β

- 4 predittori deboli con = 0.5.

j

Ognuno di questi 6 predittori di tipo è correlato (ρ = 0.5) con 2 altri tipi di predittori di tipo Questi

A B.

β

predittori correlati con hanno coefficienti = 0.

A j β

I restanti predittori (42) sono di tipo e sono solo puro rumore, con = 0 e sono indipendenti da tutti gli

C j

altri predittori.

Casual diagram 117

Un esempio potrebbe essere: gialle, zodiacale.

Y=cancro, A=fumatore, B=dita C=segno

#imposto un seme al fine di avere sempre lo stesso risultato

set.seed(1)

n <- 100 #numerosità osservazioni

p1 <- 6 #predittori di tipo A con coefficiente diverso da 0

N <- 100

p <- 60 #predittori totali

K <- 3

bb <- numeric(6*K)

bb[(0:5)*K+1] <- c(1,-1,0.5,0.5,-0.5,-0.5)

require(Matrix)

A <- 3, 3) + 0.5*diag(3)

matrix(0.5,

#matrice di correlazione di A

B <- diag(p-K*p1)

Sigma <- A, A, A, A, A, B)

bdiag(A,

R <- chol(Sigma)

X <- * p), n, p) %*% R)

as.matrix(matrix(rnorm(n

<- 1:p)

colnames(X) paste0("V",

y <- X[,1:(p1*K)] %*% bb + rnorm(n)

varType <- p)

vector("character",

varType[(0:5)*K+1] <- "A"

varType[c((0:5)*K+2, (0:5)*K+3)] <- "B"

varType[(p1*K+1):p] <- "C"

varType

## [1] "A" "B" "B" "A" "B" "B" "A" "B" "B" "A" "B" "B" "A" "B" "B" "A" "B"

## [18] "B" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C"

## [35] "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C"

## [52] "C" "C" "C" "C" "C" "C" "C" "C" "C"

#Vediamo che è presente una leggera correlazione tra A e B, con un andamento crescente

plot(X[,2],X[,1]) 118

2

1

1]

X[, 0

−1

−2 −2 −1 0 1 2

X[, 2]

#Vediamo attraverso una step forward se riesco a selezionare i predittori con coefficiente diverso da 0

yX <- x=X) #unisco i due dataframe (risposta e predittori)

data.frame(y=y,

full <- lm(y~.,yX)

null <- lm(y~1,yX)

fit <- scope=list(upper=full), direction='forward') #sembrerebbe che la forward mi prende tut

step(null,

## Start: AIC=142.64

## y ~ 1

##

## Df Sum of Sq RSS AIC

## + x.V1 1 138.538 269.59 103.17

## + x.V4 1 117.494 290.64 110.69

## + x.V16 1 47.575 360.56 132.25

## + x.V2 1 44.763 363.37 133.03

## + x.V7 1 42.929 365.20 133.53

## + x.V6 1 41.063 367.07 134.04

## + x.V5 1 29.587 378.55 137.12

## + x.V47 1 29.156 378.98 137.23

## + x.V22 1 22.691 385.44 138.92

## + x.V3 1 17.781 390.35 140.19

## + x.V36 1 14.808 393.32 140.95

## + x.V37 1 11.270 396.86 141.84

## + x.V39 1 10.592 397.54 142.01

## + x.V51 1 9.988 398.15 142.16

## + x.V30 1 9.986 398.15 142.16 119

## + x.V53 1 8.854 399.28 142.45

## <none> 408.13 142.64

## + x.V20 1 8.035 400.10 142.65

## + x.V55 1 7.922 400.21 142.68

## + x.V10 1 6.906 401.23 142.94

## + x.V33 1 6.773 401.36 142.97

## + x.V56 1 5.963 402.17 143.17

## + x.V8 1 5.724 402.41 143.23

## + x.V42 1 5.377 402.76 143.32

## + x.V50 1 5.292 402.84 143.34

## + x.V9 1 5.167 402.97 143.37

## + x.V52 1 5.104 403.03 143.38

## + x.V44 1 4.399 403.73 143.56

## + x.V15 1 4.385 403.75 143.56

## + x.V14 1 3.856 404.28 143.69

## + x.V13 1 3.299 404.83 143.83

## + x.V29 1 3.218 404.91 143.85

## + x.V34 1 3.058 405.08 143.89

## + x.V18 1 2.693 405.44 143.98

## + x.V25 1 2.329 405.80 144.07

## + x.V58 1 2.064 406.07 144.13

## + x.V46 1 2.041 406.09 144.14

## + x.V26 1 2.012 406.12 144.15

## + x.V38 1 1.853 406.28 144.19

## + x.V11 1 1.613 406.52 144.25

## + x.V57 1 1.583 406.55 144.25

## + x.V21 1 1.519 406.61 144.27

## + x.V35 1 1.467 406.67 144.28

## + x.V31 1 1.453 406.68 144.29

## + x.V59 1 1.442 406.69 144.29

## + x.V28 1 1.410 406.72 144.30

## + x.V43 1 1.290 406.84 144.33

## + x.V32 1 1.034 407.10 144.39

## + x.V19 1 1.003 407.13 144.40

## + x.V41 1 0.997 407.14 144.40

## + x.V27 1 0.820 407.31 144.44

## + x.V48 1 0.624 407.51 144.49

## + x.V49 1 0.512 407.62 144.52

## + x.V17 1 0.326 407.81 144.56

## + x.V60 1 0.230 407.90 144.59

## + x.V24 1 0.148 407.98 144.61

## + x.V40 1 0.087 408.05 144.62

## + x.V54 1 0.042 408.09 144.63

## + x.V12 1 0.032 408.10 144.63

## + x.V23 1 0.016 408.12 144.64

## + x.V45 1 0.000 408.13 144.64

##

## Step: AIC=103.17

## y ~ x.V1

##

## Df Sum of Sq RSS AIC

## + x.V4 1 105.244 164.35 55.683

## + x.V5 1 45.713 223.88 86.595

## + x.V16 1 40.871 228.72 88.734 120

## + x.V6 1 31.756 237.84 92.642

## + x.V7 1 30.882 238.71 93.009

## + x.V22 1 12.654 256.94 100.368

## + x.V37 1 11.434 258.16 100.841

## + x.V36 1 10.606 258.99 101.162

## + x.V30 1 9.844 259.75 101.455

## + x.V47 1 9.078 260.52 101.749

## + x.V39 1 8.420 261.18 102.002

## + x.V31 1 6.697 262.90 102.659

## <none> 269.60 103.175

## + x.V11 1 4.535 265.06 103.478

## + x.V51 1 4.506 265.09 103.489

## + x.V55 1 4.425 265.17 103.520

## + x.V56 1 4.098 265.50 103.643

## + x.V59 1 3.534 266.06 103.855

## + x.V58 1 3.054 266.54 104.036

## + x.V15 1 2.932 266.66 104.081

## + x.V33 1 2.900 266.69 104.093

## + x.V13 1 2.789 266.81 104.135

## + x.V18 1 2.680 266.91 104.176

## + x.V49 1 2.669 266.93 104.180

## + x.V57 1 2.608 266.99 104.203

## + x.V10 1 2.561 267.03 104.221

## + x.V26 1 2.458 267.14 104.259

## + x.V53 1 2.444 267.15 104.264

## + x.V8 1 2.441 267.15 104.265

## + x.V28 1 2.419 267.18 104.273

## + x.V44 1 2.355 267.24 104.298

## + x.V52 1 2.335 267.26 104.305

## + x.V23 1 2.168 267.43 104.368

## + x.V3 1 2.162 267.43 104.370

## + x.V29 1 1.703 267.89 104.541

## + x.V2 1 1.549 268.05 104.599

## + x.V43 1 1.277 268.32 104.700

## + x.V14 1 1.264 268.33 104.705

## + x.V46 1 1.202 268.39 104.728

## + x.V17 1 1.186 268.41 104.734

## + x.V20 1 1.157 268.44 104.745

## + x.V54 1 1.144 268.45 104.750

## + x.V9 1 1.016 268.58 104.797

## + x.V38 1 1.004 268.59 104.802

## + x.V25 1 0.683 268.91 104.921

## + x.V12 1 0.647 268.95 104.935

## + x.V24 1 0.234 269.36 105.088

## + x.V34 1 0.226 269.37 105.091

## + x.V48 1 0.141 269.45 105.122

## + x.V45 1 0.136 269.46 105.124

## + x.V50 1 0.125 269.47 105.128

## + x.V42 1 0.115 269.48 105.132

## + x.V32 1 0.060 269.53 105.153

## + x.V19 1 0.046 269.55 105.158

## + x.V35 1 0.046 269.55 105.158

## + x.V60 1 0.045 269.55 105.158

## + x.V40 1 0.034 269.56 105.162 121

## + x.V21 1 0.025 269.57 105.165

## + x.V41 1 0.001 269.59 105.175

## + x.V27 1 0.000 269.59 105.175

##

## Step: AIC=55.68

## y ~ x.V1 + x.V4

##

## Df Sum of Sq RSS AIC

## + x.V7 1 18.6154 145.74 45.662

## + x.V16 1 14.9280 149.42 48.161

## + x.V13 1 10.8423 153.51 50.858

## + x.V14 1 10.7762 153.57 50.901

## + x.V52 1 8.7691 155.58 52.200

## + x.V10 1 8.3951 155.96 52.440

## + x.V47 1 8.0933 156.26 52.633

## + x.V37 1 7.5277 156.82 52.994

## + x.V11 1 7.4693 156.88 53.032

## + x.V30 1 6.8159 157.53 53.447

## + x.V9 1 5.3892 158.96 54.349

## + x.V5 1 5.0252 159.32 54.578

## + x.V15 1 4.7064 159.64 54.777

## + x.V50 1 4.6747 159.68 54.797

## + x.V31 1 4.4263 159.92 54.953

## <none> 164.35 55.683

## + x.V26 1 3.1828 161.17 55.727

## + x.V39 1 3.1157 161.23 55.769

## + x.V59 1 3.0572 161.29 55.805

## + x.V58 1 3.0240 161.33 55.826

## + x.V22 1 2.4600 161.89 56.175

## + x.V17 1 2.3603 161.99 56.236

## + x.V19 1 2.1324 162.22 56.377

## + x.V12 1 2.1043 162.25 56.394

## + x.V60 1 1.9227 162.43 56.506

## + x.V51 1 1.6024 162.75 56.703

## + x.V46 1 1.5329 162.82 56.746

## + x.V8 1 1.3375 163.01 56.866

## + x.V36 1 1.2959 163.05 56.891

## + x.V6 1 1.1217 163.23 56.998

## + x.V33 1 1.0242 163.33 57.058

## + x.V42 1 0.9925 163.36 57.077

## + x.V55 1 0.9556 163.40 57.100

## + x.V21 1 0.9094 163.44 57.128

## + x.V35 1 0.6127 163.74 57.309

## + x.V45 1 0.5663 163.78 57.338

## + x.V29 1 0.5552 163.79 57.345

## + x.V43 1 0.5389 163.81 57.355

## + x.V18 1 0.4839 163.87 57.388

## + x.V49 1 0.4837 163.87 57.388

## + x.V34 1 0.4596 163.89 57.403

## + x.V2 1 0.4159 163.93 57.430

## + x.V32 1 0.3604 163.99 57.463

## + x.V44 1 0.3510 164.00 57.469

## + x.V3 1 0.3085 164.04 57.495

## + x.V27 1 0.2264 164.12 57.545 122

## + x.V23 1 0.1792 164.17 57.574

## + x.V56 1 0.1508 164.20 57.591

## + x.V48 1 0.1433 164.21 57.596

## + x.V28 1 0.0693 164.28 57.641

## + x.V54 1 0.0684 164.28 57.641

## + x.V57 1 0.0433 164.31 57.657

## + x.V41 1 0.0412 164.31 57.658

## + x.V25 1 0.0405 164.31 57.658

## + x.V38 1 0.0242 164.33 57.668

## + x.V24 1 0.0158 164.33 57.673

## + x.V53 1 0.0154 164.34 57.674

## + x.V20 1 0.0016 164.35 57.682

## + x.V40 1 0.0001 164.35 57.683

##

## Step: AIC=45.66

## y ~ x.V1 + x.V4 + x.V7

##

## Df Sum of Sq RSS AIC

## + x.V13 1 14.4801 131.25 37.197

## + x.V16 1 11.2837 134.45 39.603

## + x.V14 1 10.5476 135.19 40.149

## + x.V10 1 10.4938 135.24 40.189

## + x.V52 1 8.1933 137.54 41.876

## + x.V15 1 7.9952 137.74 42.019

## + x.V11 1 7.7167 138.02 42.221

## + x.V5 1 6.2000 139.53 43.314

## + x.V47 1 5.0686 140.67 44.122

## + x.V30 1 4.3980 141.34 44.598

## + x.V50 1 4.2232 141.51 44.721

## + x.V31 1 3.5623 142.17 45.187

## + x.V39 1 3.0608 142.67 45.539

## + x.V37 1 2.9795 142.75 45.596

## + x.V59 1 2.9386 142.80 45.625

## <none> 145.74 45.662

## + x.V12 1 2.6638 143.07 45.817

## + x.V51 1 2.2666 143.47 46.094

## + x.V22 1 2.2166 143.52 46.129

## + x.V26 1 2.2147 143.52 46.130

## + x.V8 1 2.2019 143.53 46.139

## + x.V36 1 2.1439 143.59 46.180

## + x.V58 1 2.0360 143.70 46.255

## + x.V19 1 1.7844 143.95 46.430

## + x.V21 1 1.6082 144.13 46.552

## + x.V17 1 1.4209 144.31 46.682

## + x.V6 1 1.3850 144.35 46.707

## + x.V29 1 1.1981 144.54 46.836

## + x.V2 1 1.1502 144.59 46.869

## + x.V46 1 1.1081 144.63 46.899

## + x.V35 1 1.1052 144.63 46.901

## + x.V34 1 0.7888 144.95 47.119

## + x.V49 1 0.7467 144.99 47.148

## + x.V45 1 0.6419 145.09 47.220

## + x.V55 1 0.4720 145.26 47.337

## + x.V48 1 0.3281 145.41 47.436 123

## + x.V44 1 0.2925 145.44 47.461

## + x.V43 1 0.2748 145.46 47.473

## + x.V18 1 0.2738 145.46 47.474

## + x.V54 1 0.2638 145.47 47.481

## + x.V42 1 0.2382 145.50 47.498

## + x.V40 1 0.2267 145.51 47.506

## + x.V33 1 0.2095 145.53 47.518

## + x.V38 1 0.1873 145.55 47.533

## + x.V28 1 0.1607 145.57 47.551

## + x.V20 1 0.1320 145.60 47.571

## + x.V57 1 0.1279 145.61 47.574

## + x.V60 1 0.1260 145.61 47.575

## + x.V3 1 0.1069 145.63 47.588

## + x.V25 1 0.0777 145.66 47.608

## + x.V24 1 0.0321 145.70 47.640

## + x.V56 1 0.0234 145.71 47.646

## + x.V23 1 0.0197 145.72 47.648

## + x.V53 1 0.0140 145.72 47.652

## + x.V9 1 0.0010 145.73 47.661

## + x.V27 1 0.0003 145.73 47.662

## + x.V32 1 0.0003 145.73 47.662

## + x.V41 1 0.0002 145.74 47.662

##

## Step: AIC=37.2

## y ~ x.V1 + x.V4 + x.V7 + x.V13

##

## Df Sum of Sq RSS AIC

## + x.V10 1 17.2820 113.97 25.079

## + x.V11 1 11.5548 119.70 29.982

## + x.V16 1 11.4258 119.83 30.089

## + x.V52 1 6.4791 124.78 34.135

## + x.V5 1 5.2660 125.99 35.102

## + x.V22 1 4.9204 126.33 35.376

## + x.V12 1 4.0127 127.24 36.092

## + x.V50 1 3.8453 127.41 36.224

## + x.V30 1 3.5983 127.66 36.417

## + x.V59 1 3.2059 128.05 36.724

## + x.V31 1 2.8580 128.40 36.995

## + x.V58 1 2.8055 128.45 37.036

## + x.V47 1 2.7778 128.48 37.058

## + x.V2 1 2.6469 128.61 37.160

## <none> 131.25 37.197

## + x.V18 1 2.4460 128.81 37.316

## + x.V8 1 2.4345 128.82 37.325

## + x.V36 1 2.2747 128.98 37.449

## + x.V19 1 2.1439 129.11 37.550

## + x.V39 1 1.9209 129.33 37.723

## + x.V40 1 1.6401 129.62 37.940

## + x.V34 1 1.5214 129.73 38.031

## + x.V35 1 1.3861 129.87 38.135

## + x.V44 1 1.3212 129.93 38.185

## + x.V15 1 1.2171 130.04 38.265

## + x.V17 1 1.1452 130.11 38.321

## + x.V26 1 1.1309 130.12 38.332 124

## + x.V38 1 1.0479 130.21 38.395

## + x.V29 1 1.0232 130.23 38.414

## + x.V37 1 0.9553 130.30 38.466

## + x.V14 1 0.9034 130.35 38.506

## + x.V51 1 0.8655 130.39 38.535

## + x.V54 1 0.7462 130.51 38.627

## + x.V6 1 0.6930 130.56 38.668

## + x.V21 1 0.6555 130.60 38.696

## + x.V46 1 0.5691 130.69 38.762

## + x.V20 1 0.4942 130.76 38.820

## + x.V48 1 0.4883 130.77 38.824

## + x.V49 1 0.3327 130.92 38.943

## + x.V42 1 0.3103 130.94 38.960

## + x.V55 1 0.2942 130.96 38.973

## + x.V45 1 0.2834 130.97 38.981

## + x.V60 1 0.2218 131.03 39.028

## + x.V24 1 0.1971 131.06 39.047

## + x.V41 1 0.1960 131.06 39.048

## + x.V56 1 0.1577 131.10 39.077

## + x.V57 1 0.1544 131.10 39.079

## + x.V32 1 0.1392 131.12 39.091

## + x.V33 1 0.1159 131.14 39.109

## + x.V27 1 0.0960 131.16 39.124

## + x.V28 1 0.0519 131.20 39.157

## + x.V43 1 0.0416 131.21 39.165

## + x.V25 1 0.0287 131.23 39.175

## + x.V53 1 0.0123 131.24 39.188

## + x.V23 1 0.0103 131.24 39.189

## + x.V3 1 0.0020 131.25 39.195

## + x.V9 1 0.0000 131.25 39.197

##

## Step: AIC=25.08

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10

##

## Df Sum of Sq RSS AIC

## + x.V16 1 18.3083 95.664 9.5675

## + x.V22 1 8.7087 105.264 19.1301

## + x.V19 1 4.5015 109.471 23.0490

## + x.V18 1 4.3372 109.635 23.1991

## + x.V52 1 4.1449 109.828 23.3743

## + x.V34 1 3.9414 110.031 23.5594

## + x.V8 1 3.6794 110.293 23.7972

## + x.V2 1 3.3725 110.600 24.0751

## + x.V5 1 3.2300 110.743 24.2039

## + x.V50 1 2.8117 111.161 24.5809

## + x.V30 1 2.6596 111.313 24.7176

## + x.V58 1 2.5842 111.388 24.7853

## <none> 113.973 25.0788

## + x.V31 1 2.1929 111.780 25.1361

## + x.V59 1 2.1105 111.862 25.2097

## + x.V11 1 2.0724 111.900 25.2437

## + x.V39 1 2.0483 111.924 25.2653

## + x.V24 1 1.8193 112.153 25.4697

## + x.V47 1 1.7797 112.193 25.5050 125

## + x.V35 1 1.5888 112.384 25.6750

## + x.V44 1 1.5512 112.421 25.7085

## + x.V29 1 1.3881 112.585 25.8534

## + x.V55 1 1.3451 112.628 25.8916

## + x.V46 1 1.3415 112.631 25.8948

## + x.V37 1 1.3170 112.656 25.9166

## + x.V36 1 1.1395 112.833 26.0740

## + x.V3 1 1.0685 112.904 26.1369

## + x.V26 1 0.8593 113.113 26.3220

## + x.V28 1 0.7018 113.271 26.4612

## + x.V57 1 0.7012 113.271 26.4617

## + x.V9 1 0.6946 113.278 26.4676

## + x.V14 1 0.6938 113.279 26.4682

## + x.V20 1 0.5385 113.434 26.6052

## + x.V56 1 0.5106 113.462 26.6298

## + x.V21 1 0.3670 113.606 26.7563

## + x.V33 1 0.3570 113.616 26.7651

## + x.V60 1 0.3376 113.635 26.7822

## + x.V51 1 0.3368 113.636 26.7829

## + x.V32 1 0.3148 113.658 26.8023

## + x.V23 1 0.2809 113.692 26.8321

## + x.V41 1 0.2276 113.745 26.8789

## + x.V49 1 0.1866 113.786 26.9150

## + x.V15 1 0.1767 113.796 26.9237

## + x.V54 1 0.1739 113.799 26.9262

## + x.V6 1 0.1545 113.818 26.9432

## + x.V53 1 0.0901 113.883 26.9998

## + x.V43 1 0.0868 113.886 27.0026

## + x.V42 1 0.0742 113.898 27.0137

## + x.V48 1 0.0394 113.933 27.0443

## + x.V27 1 0.0268 113.946 27.0553

## + x.V38 1 0.0175 113.955 27.0635

## + x.V12 1 0.0151 113.958 27.0656

## + x.V40 1 0.0149 113.958 27.0657

## + x.V25 1 0.0072 113.965 27.0725

## + x.V45 1 0.0072 113.965 27.0725

## + x.V17 1 0.0034 113.969 27.0759

##

## Step: AIC=9.57

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16

##

## Df Sum of Sq RSS AIC

## + x.V22 1 6.8005 88.864 4.1935

## + x.V17 1 4.1955 91.469 7.0828

## + x.V34 1 4.1640 91.500 7.1172

## + x.V2 1 4.0381 91.626 7.2547

## + x.V58 1 3.9362 91.728 7.3659

## + x.V19 1 3.3561 92.308 7.9963

## + x.V46 1 3.1847 92.480 8.1818

## + x.V29 1 3.1430 92.521 8.2268

## + x.V11 1 2.9561 92.708 8.4287

## + x.V5 1 2.9369 92.727 8.4494

## + x.V59 1 2.8468 92.818 8.5465

## + x.V28 1 2.0397 93.625 9.4123 126

## + x.V39 1 2.0187 93.646 9.4348

## + x.V50 1 1.9819 93.682 9.4741

## <none> 95.664 9.5675

## + x.V44 1 1.7498 93.914 9.7214

## + x.V30 1 1.6931 93.971 9.7818

## + x.V55 1 1.4774 94.187 10.0111

## + x.V52 1 1.3301 94.334 10.1674

## + x.V57 1 1.1941 94.470 10.3115

## + x.V8 1 1.0922 94.572 10.4193

## + x.V47 1 1.0413 94.623 10.4730

## + x.V9 1 1.0204 94.644 10.4951

## + x.V33 1 0.9867 94.678 10.5308

## + x.V24 1 0.9490 94.715 10.5705

## + x.V51 1 0.8895 94.775 10.6333

## + x.V14 1 0.7893 94.875 10.7391

## + x.V36 1 0.6806 94.984 10.8536

## + x.V27 1 0.5207 95.144 11.0217

## + x.V37 1 0.5153 95.149 11.0275

## + x.V20 1 0.4900 95.174 11.0540

## + x.V35 1 0.4326 95.232 11.1143

## + x.V41 1 0.4313 95.233 11.1157

## + x.V38 1 0.4177 95.247 11.1300

## + x.V31 1 0.4101 95.254 11.1379

## + x.V21 1 0.2986 95.366 11.2549

## + x.V42 1 0.2911 95.373 11.2628

## + x.V6 1 0.2384 95.426 11.3180

## + x.V3 1 0.2124 95.452 11.3453

## + x.V25 1 0.2105 95.454 11.3473

## + x.V23 1 0.2035 95.461 11.3546

## + x.V53 1 0.1720 95.492 11.3875

## + x.V40 1 0.1529 95.511 11.4076

## + x.V15 1 0.1024 95.562 11.4604

## + x.V18 1 0.1000 95.564 11.4629

## + x.V49 1 0.0868 95.577 11.4767

## + x.V32 1 0.0628 95.602 11.5018

## + x.V12 1 0.0558 95.608 11.5092

## + x.V60 1 0.0423 95.622 11.5233

## + x.V43 1 0.0325 95.632 11.5335

## + x.V56 1 0.0302 95.634 11.5360

## + x.V26 1 0.0251 95.639 11.5413

## + x.V54 1 0.0185 95.646 11.5482

## + x.V45 1 0.0176 95.647 11.5491

## + x.V48 1 0.0051 95.659 11.5622

##

## Step: AIC=4.19

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22

##

## Df Sum of Sq RSS AIC

## + x.V2 1 4.3901 84.474 1.1270

## + x.V58 1 3.4032 85.461 2.2885

## + x.V19 1 3.1842 85.680 2.5445

## + x.V34 1 3.0998 85.764 2.6430

## + x.V11 1 3.0981 85.766 2.6450

## + x.V17 1 3.0348 85.829 2.7187 127

## + x.V59 1 2.7843 86.080 3.0102

## + x.V5 1 2.6148 86.249 3.2068

## + x.V46 1 2.6032 86.261 3.2203

## + x.V39 1 2.0582 86.806 3.8501

## + x.V29 1 2.0369 86.827 3.8747

## + x.V28 1 1.9421 86.922 3.9838

## <none> 88.864 4.1935

## + x.V50 1 1.4290 87.435 4.5723

## + x.V30 1 1.3819 87.482 4.6262

## + x.V36 1 1.3633 87.501 4.6475

## + x.V8 1 1.3591 87.505 4.6522

## + x.V35 1 1.2189 87.645 4.8124

## + x.V20 1 1.1957 87.668 4.8388

## + x.V33 1 1.1755 87.688 4.8619

## + x.V44 1 1.0685 87.795 4.9838

## + x.V57 1 1.0040 87.860 5.0573

## + x.V23 1 0.9258 87.938 5.1462

## + x.V47 1 0.8846 87.979 5.1931

## + x.V3 1 0.8558 88.008 5.2258

## + x.V38 1 0.8350 88.029 5.2494

## + x.V52 1 0.7842 88.080 5.3071

## + x.V27 1 0.6444 88.219 5.4657

## + x.V55 1 0.5566 88.307 5.5651

## + x.V14 1 0.5447 88.319 5.5786

## + x.V42 1 0.5358 88.328 5.5887

## + x.V51 1 0.4972 88.367 5.6324

## + x.V15 1 0.4777 88.386 5.6545

## + x.V37 1 0.4596 88.404 5.6749

## + x.V6 1 0.4582 88.406 5.6766

## + x.V24 1 0.4556 88.408 5.6795

## + x.V54 1 0.4524 88.411 5.6831

## + x.V18 1 0.3952 88.469 5.7478

## + x.V9 1 0.3750 88.489 5.7706

## + x.V31 1 0.2906 88.573 5.8660

## + x.V60 1 0.2677 88.596 5.8917

## + x.V40 1 0.2532 88.611 5.9082

## + x.V53 1 0.2227 88.641 5.9425

## + x.V21 1 0.1884 88.675 5.9813

## + x.V43 1 0.0841 88.780 6.0988

## + x.V12 1 0.0696 88.794 6.1151

## + x.V32 1 0.0610 88.803 6.1248

## + x.V49 1 0.0500 88.814 6.1372

## + x.V25 1 0.0437 88.820 6.1443

## + x.V41 1 0.0425 88.821 6.1457

## + x.V56 1 0.0347 88.829 6.1544

## + x.V45 1 0.0115 88.852 6.1806

## + x.V26 1 0.0030 88.861 6.1901

## + x.V48 1 0.0003 88.864 6.1931

##

## Step: AIC=1.13

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2

##

## Df Sum of Sq RSS AIC

## + x.V11 1 4.0893 80.384 -1.83495 128

## + x.V19 1 3.9795 80.494 -1.69842

## + x.V5 1 3.6588 80.815 -1.30081

## + x.V46 1 2.7058 81.768 -0.12852

## + x.V35 1 2.6430 81.831 -0.05168

## + x.V29 1 2.5343 81.939 0.08099

## + x.V58 1 2.5108 81.963 0.10970

## + x.V39 1 2.3948 82.079 0.25109

## + x.V17 1 2.2280 82.246 0.45414

## + x.V59 1 2.1896 82.284 0.50084

## + x.V34 1 2.0543 82.419 0.66507

## + x.V8 1 1.8058 82.668 0.96620

## <none> 84.474 1.12705

## + x.V23 1 1.5681 82.906 1.25334

## + x.V28 1 1.5171 82.957 1.31479

## + x.V33 1 1.4265 83.047 1.42396

## + x.V38 1 1.3219 83.152 1.54978

## + x.V30 1 1.2944 83.179 1.58289

## + x.V50 1 1.2916 83.182 1.58621

## + x.V57 1 1.1866 83.287 1.71233

## + x.V52 1 1.1813 83.292 1.71881

## + x.V20 1 1.1597 83.314 1.74472

## + x.V42 1 0.7941 83.680 2.18258

## + x.V36 1 0.7298 83.744 2.25930

## + x.V51 1 0.7012 83.773 2.29347

## + x.V44 1 0.6663 83.807 2.33515

## + x.V47 1 0.6241 83.850 2.38550

## + x.V31 1 0.6165 83.857 2.39453

## + x.V24 1 0.5220 83.952 2.50720

## + x.V9 1 0.5153 83.958 2.51522

## + x.V6 1 0.5149 83.959 2.51563

## + x.V43 1 0.4795 83.994 2.55781

## + x.V55 1 0.3676 84.106 2.69090

## + x.V27 1 0.3540 84.120 2.70715

## + x.V14 1 0.3286 84.145 2.73731

## + x.V40 1 0.3200 84.154 2.74746

## + x.V15 1 0.2054 84.268 2.88356

## + x.V18 1 0.1990 84.275 2.89124

## + x.V12 1 0.1542 84.320 2.94434

## + x.V3 1 0.1491 84.325 2.95040

## + x.V54 1 0.1182 84.356 2.98700

## + x.V37 1 0.1072 84.366 3.00002

## + x.V32 1 0.0986 84.375 3.01030

## + x.V60 1 0.0866 84.387 3.02451

## + x.V41 1 0.0533 84.420 3.06394

## + x.V56 1 0.0438 84.430 3.07517

## + x.V25 1 0.0355 84.438 3.08495

## + x.V45 1 0.0302 84.444 3.09130

## + x.V48 1 0.0188 84.455 3.10474

## + x.V49 1 0.0116 84.462 3.11331

## + x.V21 1 0.0079 84.466 3.11772

## + x.V53 1 0.0058 84.468 3.12021

## + x.V26 1 0.0009 84.473 3.12601

##

## Step: AIC=-1.83 129

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11

##

## Df Sum of Sq RSS AIC

## + x.V5 1 3.9643 76.420 -4.8924

## + x.V19 1 3.7380 76.646 -4.5967

## + x.V35 1 3.1867 77.198 -3.8800

## + x.V58 1 3.0217 77.363 -3.6664

## + x.V29 1 3.0056 77.379 -3.6457

## + x.V39 1 2.6517 77.733 -3.1894

## + x.V46 1 2.4862 77.898 -2.9766

## + x.V17 1 2.4019 77.983 -2.8685

## + x.V8 1 2.2225 78.162 -2.6388

## + x.V23 1 2.1810 78.203 -2.5857

## + x.V50 1 1.7573 78.627 -2.0453

## <none> 80.384 -1.8349

## + x.V20 1 1.4951 78.889 -1.7125

## + x.V59 1 1.4632 78.921 -1.6720

## + x.V34 1 1.1935 79.191 -1.3308

## + x.V52 1 1.1557 79.229 -1.2832

## + x.V28 1 1.1269 79.258 -1.2467

## + x.V38 1 1.0815 79.303 -1.1894

## + x.V44 1 1.0576 79.327 -1.1594

## + x.V27 1 0.8726 79.512 -0.9264

## + x.V31 1 0.7384 79.646 -0.7578

## + x.V43 1 0.6596 79.725 -0.6588

## + x.V36 1 0.6392 79.745 -0.6333

## + x.V42 1 0.6320 79.752 -0.6243

## + x.V57 1 0.6194 79.765 -0.6085

## + x.V51 1 0.5762 79.808 -0.5543

## + x.V47 1 0.5552 79.829 -0.5280

## + x.V14 1 0.5407 79.844 -0.5098

## + x.V30 1 0.4812 79.903 -0.4354

## + x.V33 1 0.4750 79.909 -0.4276

## + x.V15 1 0.4572 79.927 -0.4053

## + x.V32 1 0.4095 79.975 -0.3457

## + x.V6 1 0.3698 80.015 -0.2961

## + x.V24 1 0.3174 80.067 -0.2306

## + x.V9 1 0.2749 80.109 -0.1776

## + x.V12 1 0.2717 80.113 -0.1736

## + x.V55 1 0.2054 80.179 -0.0908

## + x.V3 1 0.1954 80.189 -0.0784

## + x.V54 1 0.1932 80.191 -0.0756

## + x.V41 1 0.1351 80.249 -0.0031

## + x.V53 1 0.1280 80.256 0.0057

## + x.V48 1 0.0968 80.288 0.0446

## + x.V40 1 0.0947 80.290 0.0471

## + x.V45 1 0.0806 80.304 0.0647

## + x.V37 1 0.0534 80.331 0.0986

## + x.V60 1 0.0502 80.334 0.1026

## + x.V18 1 0.0219 80.363 0.1378

## + x.V49 1 0.0201 80.364 0.1401

## + x.V56 1 0.0151 80.369 0.1463

## + x.V26 1 0.0140 80.370 0.1477 130

## + x.V25 1 0.0019 80.383 0.1627

## + x.V21 1 0.0011 80.383 0.1636

##

## Step: AIC=-4.89

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11 + x.V5

##

## Df Sum of Sq RSS AIC

## + x.V19 1 3.3811 73.039 -7.4176

## + x.V29 1 3.2743 73.146 -7.2715

## + x.V35 1 2.8431 73.577 -6.6837

## + x.V58 1 2.5326 73.888 -6.2626

## + x.V23 1 2.5280 73.892 -6.2564

## + x.V17 1 2.4624 73.958 -6.1676

## + x.V46 1 2.0554 74.365 -5.6188

## + x.V52 1 1.8675 74.553 -5.3665

## + x.V59 1 1.8295 74.591 -5.3154

## + x.V8 1 1.7270 74.693 -5.1781

## <none> 76.420 -4.8924

## + x.V39 1 1.3368 75.083 -4.6571

## + x.V50 1 1.2961 75.124 -4.6030

## + x.V44 1 1.2152 75.205 -4.4953

## + x.V28 1 1.1434 75.277 -4.3999

## + x.V20 1 1.0716 75.349 -4.3045

## + x.V34 1 0.9267 75.493 -4.1124

## + x.V47 1 0.9215 75.499 -4.1055

## + x.V27 1 0.8436 75.577 -4.0025

## + x.V42 1 0.8204 75.600 -3.9717

## + x.V36 1 0.8163 75.604 -3.9663

## + x.V51 1 0.7428 75.677 -3.8692

## + x.V33 1 0.7083 75.712 -3.8236

## + x.V38 1 0.5923 75.828 -3.6705

## + x.V31 1 0.4723 75.948 -3.5123

## + x.V32 1 0.4528 75.967 -3.4867

## + x.V57 1 0.4377 75.982 -3.4668

## + x.V12 1 0.3901 76.030 -3.4042

## + x.V53 1 0.3512 76.069 -3.3530

## + x.V14 1 0.3445 76.076 -3.3442

## + x.V37 1 0.2925 76.128 -3.2758

## + x.V24 1 0.2505 76.170 -3.2207

## + x.V30 1 0.2483 76.172 -3.2178

## + x.V15 1 0.2351 76.185 -3.2004

## + x.V43 1 0.2251 76.195 -3.1873

## + x.V3 1 0.1982 76.222 -3.1520

## + x.V40 1 0.1842 76.236 -3.1337

## + x.V55 1 0.1512 76.269 -3.0904

## + x.V48 1 0.1487 76.271 -3.0872

## + x.V41 1 0.1343 76.286 -3.0683

## + x.V26 1 0.1165 76.304 -3.0450

## + x.V54 1 0.1056 76.315 -3.0306

## + x.V60 1 0.0864 76.334 -3.0056

## + x.V9 1 0.0855 76.335 -3.0043

## + x.V45 1 0.0728 76.347 -2.9877

## + x.V25 1 0.0339 76.386 -2.9368 131

## + x.V56 1 0.0285 76.392 -2.9296

## + x.V49 1 0.0214 76.399 -2.9204

## + x.V18 1 0.0093 76.411 -2.9046

## + x.V6 1 0.0027 76.417 -2.8960

## + x.V21 1 0.0023 76.418 -2.8953

##

## Step: AIC=-7.42

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11 + x.V5 + x.V19

##

## Df Sum of Sq RSS AIC

## + x.V23 1 3.5522 69.487 -10.4033

## + x.V29 1 3.3273 69.712 -10.0802

## + x.V35 1 2.6289 70.410 -9.0834

## + x.V17 1 2.5377 70.501 -8.9538

## + x.V52 1 2.1231 70.916 -8.3675

## + x.V58 1 2.0783 70.961 -8.3043

## + x.V59 1 1.8488 71.190 -7.9815

## + x.V8 1 1.5944 71.445 -7.6248

## + x.V46 1 1.5173 71.522 -7.5169

## <none> 73.039 -7.4176

## + x.V50 1 1.3887 71.650 -7.3372

## + x.V39 1 1.3234 71.716 -7.2461

## + x.V44 1 1.2542 71.785 -7.1497

## + x.V28 1 1.2286 71.810 -7.1141

## + x.V47 1 1.1882 71.851 -7.0578

## + x.V42 1 1.0537 71.985 -6.8708

## + x.V34 1 1.0126 72.026 -6.8137

## + x.V9 1 0.6983 72.341 -6.3783

## + x.V51 1 0.6726 72.366 -6.3428

## + x.V20 1 0.6651 72.374 -6.3324

## + x.V38 1 0.6636 72.375 -6.3303

## + x.V31 1 0.6273 72.412 -6.2802

## + x.V53 1 0.6031 72.436 -6.2467

## + x.V32 1 0.5849 72.454 -6.2216

## + x.V37 1 0.4672 72.572 -6.0594

## + x.V15 1 0.4247 72.614 -6.0008

## + x.V12 1 0.3890 72.650 -5.9516

## + x.V27 1 0.3712 72.668 -5.9271

## + x.V14 1 0.3642 72.675 -5.9175

## + x.V40 1 0.3416 72.697 -5.8864

## + x.V30 1 0.3021 72.737 -5.8321

## + x.V33 1 0.2728 72.766 -5.7918

## + x.V24 1 0.2515 72.788 -5.7625

## + x.V60 1 0.2492 72.790 -5.7594

## + x.V43 1 0.2158 72.823 -5.7135

## + x.V55 1 0.2091 72.830 -5.7044

## + x.V18 1 0.2044 72.835 -5.6979

## + x.V26 1 0.1820 72.857 -5.6671

## + x.V54 1 0.1388 72.900 -5.6078

## + x.V45 1 0.1120 72.927 -5.5711

## + x.V36 1 0.1119 72.927 -5.5709

## + x.V48 1 0.1069 72.932 -5.5642

## + x.V57 1 0.0896 72.949 -5.5404 132

## + x.V3 1 0.0611 72.978 -5.5013

## + x.V21 1 0.0483 72.991 -5.4838

## + x.V56 1 0.0126 73.026 -5.4349

## + x.V6 1 0.0064 73.033 -5.4264

## + x.V41 1 0.0059 73.033 -5.4257

## + x.V49 1 0.0056 73.033 -5.4253

## + x.V25 1 0.0020 73.037 -5.4204

##

## Step: AIC=-10.4

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11 + x.V5 + x.V19 + x.V23

##

## Df Sum of Sq RSS AIC

## + x.V17 1 3.1872 66.300 -13.0986

## + x.V35 1 2.7333 66.754 -12.4162

## + x.V29 1 2.3373 67.150 -11.8247

## + x.V59 1 2.2095 67.277 -11.6347

## + x.V58 1 1.9137 67.573 -11.1959

## + x.V8 1 1.8774 67.609 -11.1423

## + x.V52 1 1.6686 67.818 -10.8339

## + x.V47 1 1.6558 67.831 -10.8151

## + x.V34 1 1.4178 68.069 -10.4647

## <none> 69.487 -10.4033

## + x.V39 1 1.3126 68.174 -10.3104

## + x.V46 1 1.2923 68.195 -10.2806

## + x.V9 1 1.0871 68.400 -9.9801

## + x.V51 1 0.9850 68.502 -9.8310

## + x.V20 1 0.9167 68.570 -9.7314

## + x.V28 1 0.8356 68.651 -9.6131

## + x.V12 1 0.8349 68.652 -9.6121

## + x.V31 1 0.6805 68.806 -9.3875

## + x.V44 1 0.6738 68.813 -9.3777

## + x.V14 1 0.5447 68.942 -9.1902

## + x.V32 1 0.5375 68.949 -9.1798

## + x.V33 1 0.5165 68.970 -9.1493

## + x.V50 1 0.4336 69.053 -9.0292

## + x.V42 1 0.4133 69.074 -8.9998

## + x.V53 1 0.3798 69.107 -8.9514

## + x.V15 1 0.3108 69.176 -8.8516

## + x.V27 1 0.2812 69.206 -8.8088

## + x.V40 1 0.2797 69.207 -8.8067

## + x.V54 1 0.2770 69.210 -8.8027

## + x.V25 1 0.2590 69.228 -8.7767

## + x.V38 1 0.2548 69.232 -8.7706

## + x.V55 1 0.2506 69.236 -8.7646

## + x.V37 1 0.2372 69.250 -8.7452

## + x.V30 1 0.2226 69.264 -8.7241

## + x.V43 1 0.2152 69.272 -8.7134

## + x.V24 1 0.2042 69.283 -8.6976

## + x.V56 1 0.1721 69.315 -8.6512

## + x.V3 1 0.1707 69.316 -8.6492

## + x.V60 1 0.1470 69.340 -8.6150

## + x.V49 1 0.1423 69.345 -8.6083

## + x.V36 1 0.1248 69.362 -8.5831 133

## + x.V45 1 0.1122 69.375 -8.5648

## + x.V26 1 0.0974 69.389 -8.5436

## + x.V6 1 0.0363 69.451 -8.4555

## + x.V18 1 0.0316 69.455 -8.4487

## + x.V48 1 0.0226 69.464 -8.4358

## + x.V57 1 0.0099 69.477 -8.4175

## + x.V41 1 0.0048 69.482 -8.4102

## + x.V21 1 0.0042 69.483 -8.4093

##

## Step: AIC=-13.1

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11 + x.V5 + x.V19 + x.V23 + x.V17

##

## Df Sum of Sq RSS AIC

## + x.V35 1 2.80761 63.492 -15.426

## + x.V59 1 2.70452 63.595 -15.263

## + x.V29 1 1.99076 64.309 -14.147

## + x.V34 1 1.42486 64.875 -13.271

## + x.V39 1 1.41168 64.888 -13.251

## + x.V12 1 1.38127 64.918 -13.204

## + x.V52 1 1.37041 64.929 -13.187

## + x.V47 1 1.36344 64.936 -13.177

## <none> 66.300 -13.099

## + x.V31 1 1.28612 65.014 -13.057

## + x.V58 1 1.28091 65.019 -13.050

## + x.V51 1 1.20765 65.092 -12.937

## + x.V18 1 1.01110 65.289 -12.635

## + x.V8 1 0.95927 65.340 -12.556

## + x.V14 1 0.94456 65.355 -12.534

## + x.V44 1 0.92546 65.374 -12.504

## + x.V28 1 0.83810 65.462 -12.371

## + x.V9 1 0.79876 65.501 -12.311

## + x.V32 1 0.75431 65.545 -12.243

## + x.V33 1 0.68318 65.616 -12.134

## + x.V37 1 0.64796 65.652 -12.081

## + x.V20 1 0.41992 65.880 -11.734

## + x.V46 1 0.39896 65.901 -11.702

## + x.V27 1 0.37263 65.927 -11.662

## + x.V25 1 0.37137 65.928 -11.660

## + x.V24 1 0.29400 66.006 -11.543

## + x.V43 1 0.29015 66.009 -11.537

## + x.V55 1 0.28434 66.015 -11.528

## + x.V49 1 0.25211 66.048 -11.480

## + x.V50 1 0.22942 66.070 -11.445

## + x.V42 1 0.22776 66.072 -11.443

## + x.V40 1 0.21158 66.088 -11.418

## + x.V53 1 0.20373 66.096 -11.406

## + x.V30 1 0.19914 66.100 -11.399

## + x.V60 1 0.16364 66.136 -11.346

## + x.V15 1 0.14823 66.151 -11.322

## + x.V54 1 0.12432 66.175 -11.286

## + x.V57 1 0.09599 66.204 -11.243

## + x.V36 1 0.09067 66.209 -11.235

## + x.V41 1 0.08835 66.211 -11.232 134

## + x.V38 1 0.06808 66.232 -11.201

## + x.V45 1 0.04147 66.258 -11.161

## + x.V3 1 0.03872 66.261 -11.157

## + x.V26 1 0.02263 66.277 -11.133

## + x.V6 1 0.01814 66.281 -11.126

## + x.V48 1 0.01669 66.283 -11.124

## + x.V56 1 0.00991 66.290 -11.114

## + x.V21 1 0.00009 66.300 -11.099

##

## Step: AIC=-15.43

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11 + x.V5 + x.V19 + x.V23 + x.V17 + x.V35

##

## Df Sum of Sq RSS AIC

## + x.V59 1 3.6328 59.859 -19.317

## + x.V47 1 1.8540 61.638 -16.389

## + x.V39 1 1.7559 61.736 -16.230

## + x.V58 1 1.4896 62.002 -15.800

## + x.V18 1 1.4304 62.062 -15.704

## + x.V51 1 1.3918 62.100 -15.642

## + x.V12 1 1.3865 62.106 -15.633

## + x.V52 1 1.3177 62.174 -15.523

## + x.V29 1 1.2954 62.197 -15.487

## <none> 63.492 -15.426

## + x.V32 1 1.1960 62.296 -15.327

## + x.V44 1 0.9531 62.539 -14.938

## + x.V9 1 0.9125 62.579 -14.873

## + x.V34 1 0.9062 62.586 -14.863

## + x.V8 1 0.8541 62.638 -14.780

## + x.V55 1 0.8380 62.654 -14.754

## + x.V28 1 0.6253 62.867 -14.415

## + x.V37 1 0.5646 62.927 -14.319

## + x.V33 1 0.5415 62.951 -14.282

## + x.V31 1 0.5320 62.960 -14.267

## + x.V20 1 0.4071 63.085 -14.069

## + x.V24 1 0.3748 63.117 -14.018

## + x.V46 1 0.3533 63.139 -13.983

## + x.V14 1 0.3456 63.146 -13.971

## + x.V25 1 0.3419 63.150 -13.966

## + x.V49 1 0.2926 63.199 -13.887

## + x.V50 1 0.2627 63.229 -13.840

## + x.V60 1 0.2464 63.246 -13.814

## + x.V42 1 0.2256 63.266 -13.782

## + x.V43 1 0.2135 63.279 -13.762

## + x.V48 1 0.2085 63.283 -13.755

## + x.V53 1 0.1571 63.335 -13.673

## + x.V27 1 0.1476 63.344 -13.658

## + x.V36 1 0.1295 63.362 -13.630

## + x.V54 1 0.1250 63.367 -13.623

## + x.V30 1 0.1160 63.376 -13.608

## + x.V21 1 0.0634 63.429 -13.525

## + x.V40 1 0.0542 63.438 -13.511

## + x.V38 1 0.0531 63.439 -13.509

## + x.V15 1 0.0520 63.440 -13.508 135

## + x.V26 1 0.0310 63.461 -13.474

## + x.V41 1 0.0208 63.471 -13.458

## + x.V6 1 0.0144 63.478 -13.448

## + x.V56 1 0.0100 63.482 -13.441

## + x.V3 1 0.0100 63.482 -13.441

## + x.V57 1 0.0042 63.488 -13.432

## + x.V45 1 0.0001 63.492 -13.426

##

## Step: AIC=-19.32

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11 + x.V5 + x.V19 + x.V23 + x.V17 + x.V35 + x.V59

##

## Df Sum of Sq RSS AIC

## + x.V39 1 1.93693 57.922 -20.607

## + x.V47 1 1.37635 58.483 -19.644

## + x.V18 1 1.30893 58.550 -19.529

## + x.V29 1 1.27040 58.589 -19.463

## + x.V12 1 1.24447 58.615 -19.418

## <none> 59.859 -19.317

## + x.V28 1 1.17123 58.688 -19.294

## + x.V58 1 1.09566 58.764 -19.165

## + x.V55 1 1.08083 58.778 -19.140

## + x.V52 1 1.04578 58.813 -19.080

## + x.V51 1 0.98220 58.877 -18.972

## + x.V32 1 0.87018 58.989 -18.782

## + x.V25 1 0.80453 59.055 -18.671

## + x.V46 1 0.77729 59.082 -18.625

## + x.V34 1 0.70312 59.156 -18.499

## + x.V44 1 0.67936 59.180 -18.459

## + x.V27 1 0.66943 59.190 -18.442

## + x.V37 1 0.61838 59.241 -18.356

## + x.V31 1 0.58776 59.271 -18.304

## + x.V33 1 0.53651 59.323 -18.218

## + x.V9 1 0.50059 59.359 -18.157

## + x.V49 1 0.39186 59.467 -17.974

## + x.V24 1 0.37717 59.482 -17.950

## + x.V48 1 0.37521 59.484 -17.946

## + x.V43 1 0.31017 59.549 -17.837

## + x.V8 1 0.29249 59.567 -17.807

## + x.V50 1 0.28110 59.578 -17.788

## + x.V40 1 0.26585 59.593 -17.763

## + x.V42 1 0.25471 59.604 -17.744

## + x.V54 1 0.15284 59.706 -17.573

## + x.V20 1 0.11808 59.741 -17.515

## + x.V56 1 0.11194 59.747 -17.505

## + x.V21 1 0.10987 59.749 -17.501

## + x.V60 1 0.10893 59.750 -17.500

## + x.V14 1 0.09172 59.767 -17.471

## + x.V53 1 0.07547 59.784 -17.444

## + x.V41 1 0.07284 59.786 -17.439

## + x.V36 1 0.05982 59.799 -17.418

## + x.V15 1 0.04774 59.811 -17.397

## + x.V6 1 0.04331 59.816 -17.390

## + x.V26 1 0.02312 59.836 -17.356 136

## + x.V57 1 0.02188 59.837 -17.354

## + x.V38 1 0.01461 59.845 -17.342

## + x.V30 1 0.01395 59.845 -17.341

## + x.V3 1 0.00293 59.856 -17.322

## + x.V45 1 0.00045 59.859 -17.318

##

## Step: AIC=-20.61

## y ~ x.V1 + x.V4 + x.V7 + x.V13 + x.V10 + x.V16 + x.V22 + x.V2 +

## x.V11 + x.V5 + x.V19 + x.V23 + x.V17 + x.V35 + x.V59 + x.V39

##

## Df Sum of Sq RSS AIC

## + x.V47 1 1.30156 56.621 -20.880

## <none> 57.922 -20.607

## + x.V29 1 1.10627 56.816 -20.535

## + x.V55 1 1.02988 56.892 -20.401

## + x.V18 1 1.00878 56.913 -20.364

## + x.V12 1 0.97812 56.944 -20.310

## + x.V46 1 0.90031 57.022 -20.173

## + x.V34 1 0.87286 57.049 -20.125

## + x.V28 1 0.82712 57.095 -20.045

## + x.V25 1 0.72402 57.198 -19.865

## + x.V51 1 0.71808 57.204 -19.854

## + x.V42 1 0.63368 57.289 -19.707

## + x.V52 1 0.59891 57.323 -19.646

## + x.V32 1 0.58139 57.341 -19.616

## + x.V58 1 0.54252 57.380 -19.548

## + x.V44 1 0.52777 57.394 -19.522

## + x.V50 1 0.52404 57.398 -19.516

## + x.V9 1 0.52166 57.401 -19.512

## + x.V37 1 0.47392 57.448 -19.428

## + x.V27 1 0.46773 57.455 -19.418

## + x.V24 1 0.45916 57.463 -19.403

## + x.V33 1 0.39550 57.527 -19.292

## + x.V40 1 0.37340 57.549 -19.254

## + x.V43 1 0.30143 57.621 -19.129

## + x.V21 1 0.26859 57.654 -19.072

## + x.V48 1 0.21498 57.707 -18.979

## + x.V31 1 0.19758 57.725 -18.948

## + x.V60 1 0.16558 57.757 -18.893

## + x.V49 1 0.15235 57.770 -18.870

## + x.V20 1 0.14993 57.772 -18.866

## + x.V36 1 0.11986 57.802 -18.814

## + x.V8 1 0.11790 57.804 -18.811

## + x.V54 1 0.09775 57.825 -18.776

## + x.V6 1 0.07695 57.845 -18.740

## + x.V14 1 0.06807 57.854 -18.724

## + x.V53 1 0.06518 57.857 -18.719

## + x.V26 1 0.05950 57.863 -18.710

## + x.V38 1 0.05562 57.867 -18.703

## + x.V45 1 0.04262 57.880 -18.680

## + x.V30 1 0.03901 57.883 -18.674

## + x.V41 1 0.03584 57.886 -18.669

## + x.V3 1 0.01754 57.905 -18.637

## + x.V56 1 0.01605 57.906 -18.635 137


PAGINE

157

PESO

5.41 MB

AUTORE

Pagani21

PUBBLICATO

+1 anno fa


DETTAGLI
Corso di laurea: Corso di laurea magistrale in scienze statistiche ed economiche
SSD:
Docente: Solari Aldo
A.A.: 2017-2018

I contenuti di questa pagina costituiscono rielaborazioni personali del Publisher Pagani21 di informazioni apprese con la frequenza delle lezioni di Statistical learning e studio autonomo di eventuali libri di riferimento in preparazione dell'esame finale o della tesi. Non devono intendersi come materiale ufficiale dell'università Milano Bicocca - Unimib o del prof Solari Aldo.

Acquista con carta o conto PayPal

Scarica il file tutte le volte che vuoi

Paga con un conto PayPal per usufruire della garanzia Soddisfatto o rimborsato

Recensioni
Ti è piaciuto questo appunto? Valutalo!

Altri appunti di Corso di laurea magistrale in scienze statistiche ed economiche

Statistica Aziendale M
Appunto
Previsione serie storiche applicando la foresta casuale
Tesi
Economia Applicata M
Esercitazione
Schema Statistica economica M
Appunto