> source("T:/Alumnes/Román/Clase_20304/bootstrap-t.r")

> source("T:/Alumnes/Román/Clase_20304/rbivnorm.r")

> law.school<-read.table("Law School.txt",header=T)

> n <- nrow(law.school) # n es el tamaño muestral

> indices <- sample(1:n, replace=T) # "indices" define una remuestra

> indices

[1] 75 25 19 42 60 43 15 10 72 6 40 19 51 46 3 16 54 77 47 26 73 56 32 5 57

[26] 77 56 37 54 44 55 67 52 63 68 65 35 4 77 53 26 26 31 66 58 23 7 26 38 59

[51] 68 10 6 32 55 79 73 61 49 53 2 32 9 41 1 5 50 12 20 59 36 14 76 62 52

[76] 81 53 1 79 40 56 35

> correl(law.school, indices) # correlacion de la remuestra definida por "indices", calculada de otra forma

[1] 0.7021442

> correl(law.school, 1:n)

[1] 0.7599979

#hay otra manera alternativa para el bootstrap no parametrico

#es emplear vectores no paramétricos

> presamp <-tabulate(indices)

> indices

[1] 75 25 19 42 60 43 15 10 72 6 40 19 51 46 3 16 54 77 47 26 73 56 32 5 57

[26] 77 56 37 54 44 55 67 52 63 68 65 35 4 77 53 26 26 31 66 58 23 7 26 38 59

[51] 68 10 6 32 55 79 73 61 49 53 2 32 9 41 1 5 50 12 20 59 36 14 76 62 52

[76] 81 53 1 79 40 56 35

> presamp

[1] 2 1 1 1 2 2 1 0 1 2 0 1 0 1 1 1 0 0 2 1 0 0 1 0 1 4 0 0 0 0 1 3 0 0 2 1 1 1

[39] 0 2 1 1 1 1 0 1 1 0 1 1 1 2 3 2 2 3 1 1 2 1 1 1 1 0 1 1 1 2 0 0 0 1 2 0 1 1

[77] 3 0 2 0 1

# Para ordenar los indices de menor a mayor

> sort(indices)

[1] 1 1 2 3 4 5 5 6 6 7 9 10 10 12 14 15 16 19 19 20 23 25 26 26 26

[26] 26 31 32 32 32 35 35 36 37 38 40 40 41 42 43 44 46 47 49 50 51 52 52 53 53

[51] 53 54 54 55 55 56 56 56 57 58 59 59 60 61 62 63 65 66 67 68 68 72 73 73 75

[76] 76 77 77 77 79 79 81

> presamp <- presamp / n

> presamp

[1] 0.02439024 0.01219512 0.01219512 0.01219512 0.02439024 0.02439024

[7] 0.01219512 0.00000000 0.01219512 0.02439024 0.00000000 0.01219512

[13] 0.00000000 0.01219512 0.01219512 0.01219512 0.00000000 0.00000000

[19] 0.02439024 0.01219512 0.00000000 0.00000000 0.01219512 0.00000000

[25] 0.01219512 0.04878049 0.00000000 0.00000000 0.00000000 0.00000000

[31] 0.01219512 0.03658537 0.00000000 0.00000000 0.02439024 0.01219512

[37] 0.01219512 0.01219512 0.00000000 0.02439024 0.01219512 0.01219512

[43] 0.01219512 0.01219512 0.00000000 0.01219512 0.01219512 0.00000000

[49] 0.01219512 0.01219512 0.01219512 0.02439024 0.03658537 0.02439024

[55] 0.02439024 0.03658537 0.01219512 0.01219512 0.02439024 0.01219512

[61] 0.01219512 0.01219512 0.01219512 0.00000000 0.01219512 0.01219512

[67] 0.01219512 0.02439024 0.00000000 0.00000000 0.00000000 0.01219512

[73] 0.02439024 0.00000000 0.01219512 0.01219512 0.03658537 0.00000000

[79] 0.02439024 0.00000000 0.01219512

#en un contexto de remuestreo no paramétrico tiene que ser lo mismo

#la función correl hace elmismo cálculo

> correl(law.school, p=presamp)

[1] 0.7356295

Warning messages:

1: longer object length

is not a multiple of shorter object length in: x * p

2: longer object length

is not a multiple of shorter object length in: y * p

3: longer object length

is not a multiple of shorter object length in: y * p

4: longer object length

is not a multiple of shorter object length in: x * x * p

> correl

function (datos, indices, p, ...) {

if (missing (indices)){

if (missing (p))

return (cor (datos[,1], datos[,2]))

else {

x <- datos[,1]

y <- datos[,2]

x <- x - sum (x * p)

y <- y - sum (y * p)

yp <- y*p

return ((sum(x*yp))/sqrt ((sum(x*x*p))*(sum(y*yp))))

}

}

else

return (cor (datos[indices,])[1,2])

}

# los puntos suspensivos quieren decir que significa que (ver después del # # return) llamaria a la función y lo pasaria donde los puntos suspensivos

# Definicion de se.correl, calculo del error estandar de la correlacion #muestral

> se.correl

function(datos, rho, se.fact, ...){

if (missing (rho))

rho <- correl(datos, ...)

if (missing (se.fact))

se.fact <- 1 / sqrt (nrow(datos) - 3)

return ((1 - rho * rho) * se.fact)

}

> se.r <- se.correl(law.school) # A partir de los datos originales

> se.correl(law.school, rho=r) # Sin necesidad de calcular de nuevo la correlacion

Error in se.correl(law.school, rho = r) : Object "r" not found

> se.const <- 1 / sqrt(n - 3)

> se.correl(rho=r, se.fact=se.const)

Error in se.correl(rho = r, se.fact = se.const) :

Object "r" not found

#hay una posibilidad de la que hablaremos de aquí a unos dias

#que es una manera reciente de hallar errores estandard y que

#es mucho mas cómodo plantearla empleando vectores de remuestreo.

#puede cortocircuitar, por razones de eficiencia, algunas cosas, si le damos

#rho ya no lo calcula y si le damos el denominador de la expresión, ya no

#lo calculose.correl

> se.correl

function(datos, rho, se.fact, ...){

if (missing (rho))

rho <- correl(datos, ...)

if (missing (se.fact))

se.fact <- 1 / sqrt (nrow(datos) - 3)

return ((1 - rho * rho) * se.fact)

}

#Finalmente tenemos otra función que calcula la correlación studentizada

> stud.correl

function(datos, rho.samp, rho, se, ...){

if (missing (rho.samp))

rho.samp <- correl(datos, ...)

if (missing (se))

se <- se.correl(datos, rho = rho.samp, ...)

return ((rho.samp - rho) / se)

}

> stud.boot.correl

function(datos, B = 999){

if (B < 0) stop("Numero negativo o nulo de replicas bootstrap")

tboot <- numeric(length=B)

n <- nrow(datos)

r <- correl(datos)

se.const <- 1 / sqrt (n - 3)

se.r <- se.correl(rho=r, se.fact=se.const)

indexos <- 1:n

for (b in 1:B) {

tboot[b] <- stud.correl(datos[sample(indexos, replace=T),],rho=r, se.fact=se.const)

}

attr(tboot, "rho") <- r

attr(tboot, "se.rho") <- se.r

return (tboot)

}

#la función attr permite definir atributos de cualquier objeto y de manera #dinámica crea o permite calcular atributos

> tboot <- stud.boot.correl(law.school)

attr(,"rho")

[1] 0.7599979

attr(,"se.rho")

[1] 0.04752408

#este tboot lo emplearemos como un vector, si decimos siete nos

#referimos a la posición siete del vector. Para recuperarlo llamamos

> attr(tboot, "rho")

[1] 0.7599979

> r <- attr(tboot, "rho")

> icboot.t.correl

function(datos, alfa=0.05, simetrizado=F, B=999){

remuestreo <- boot.t.correl( datos, B)

if (simetrizado) {

tU <- quantile (abs (remuestreo), probs = 1 - alfa)

tL <- -tU

}

else {

alfa2 <- 0.5*alfa

t.boot <- quantile (remuestreo, probs = c(alfa2,1-alfa2))

tL <- t.boot[1]

tU <- t.boot[2]

}

rho <- attr(remuestreo,"rho")

se.rho <- attr(remuestreo,"se.rho")

ic <- c(rho - tU*se.rho, rho - tL*se.rho)

names(ic) <- NULL

return(ic)

}

#finalmente esta función nos devuelve el intervalo de confianza bootstrap-t

> icboot.t.correl(law.school)

[1] 0.6272623 0.8409062

> icboot.t.correl(law.school, alfa=0.01)

[1] 0.5902067 0.8628937

> icboot.t.correl(law.school,simetrizado=T)

[1] 0.6572470 0.8627487

> icboot.t.correl(law.school,sim=T)

[1] 0.6463448 0.8736509

#están hechas para la correlación pero a poco que se utilicen pueden

#emplearse para cualquier estadístico, con ciertas limitaciones, etc.

#Para el caso concreto de datos normal bivariante y el coeficiente de correlación de Pearson tenemos el #cambio de Fisher que si tengo rho puedo cambiar de escala mediante la siguiente transformación:

{ {

> #finalmente esta función nos devuelve el intervalo de confianza bootstrap-t

> icboot.t.correl(law.school)

[1] 0.6272623 0.8409062

> icboot.t.correl(law.school, alfa=0.01)

[1] 0.5902067 0.8628937

> icboot.t.correl(law.school,simetrizado=T)

[1] 0.6572470 0.8627487

> icboot.t.correl(law.school,sim=T)

[1] 0.6463448 0.8736509

#están hechas para la correlación pero a poco que se utilicen pueden

#emplearse para cualquier estadístico, con ciertas limitaciones, etc.

> icnorm.correl(law.school)

[1] 0.6502299 0.8386849

#lo da mucho más ràpidamente, como ha podido verse.

#Para saber cual de las anteriores funciones funciona mejor hay que hacer #simulación, que es lo que intentan hacer las simulaciones que siguen

> sim.boot.t.correl

function(nsim = 1000, n=10, m1=0, m2=0, s1=1, s2=1, rho=0.5, alfa=0.05, simetrizado=F, B = 999){

if (nsim <= 0) stop("Numero negativo o nulo de replicas de simulacion")

if (n <= 0) stop("Tamaño muestral negativo o nulo")

if (s1 <= 0 || s2 <= 0) stop("La desviacion estandar no puede ser negativa o nula")

if (rho < -1.0 || rho > 1.0) stop("La correlacion tiene que estar entre -1 i +1")

datos.sim <- rbivnorm(nsim*n, m1, m2, s1, s2, rho)

sim.results <- matrix (NA, nrow=3, ncol=nsim)

for (isim in 1:nsim) {

sim.results[,isim] <- ic.sim (isim, rho, n, datos.sim, alfa, simetrizado, B)

}

return (t (sim.results))

# Equivalente, supuestamente mas eficiente:

# matrix (unlist (lapply (1:nsim, ic.sim, rho, n, datos.sim, alfa, simetrizado, B)), ncol=3, byrow=T)

}

# Solo valen para el caso de querer generar normales bivariantes.

>Sim.norm.correl

>m1=mean (law.school [ ,1 ] )

>m2=mean (law.school [ ,2 ] )

>s1=mean (law.school [ ,1 ] )

>s2=mean (law.school [ ,2 ] )

# realizamos la simulación

>sim.results.norm < - sim.nor.correl (n = n, m1 = m1, m2 = m2, s1 = s1, s2 = s2, rho = 0.8 )

>sim.results.norm

>mean (sim.results.norm ) [ , 3 ] )

>mean (sim.results.norm ) [ , 2 ] )

>mean (sim.results.norm ) [ , 2 ] ) – sim.results.norm [ , 1 ] )

COMPARACIÓN DE DOS GRUPOS (medias) CON EL EJEMPLO DE CHACALES

El contraste de permutaciones era adecuado para este estúdio y veremos como el método Boostrap tiene unas ventajas o desventajas.

Ventajas.

Caso permutaciones: No tenemos porqué tener una muestra aleatoria de la población. Es un test condicional pero no presupone aleatorización.

Desventajas

Boostrap: Presupone que detrás d elos datos hay un modelo probabilístico, en definitiva, un muestreo aleatório simple muy laxo ya que no hay que fijar la distribución pero existe.

#Ejecutaremos el fichero Contrastes Boostrap.r paso a paso.

#Generamos los datos de los chacales del ejemplo

sexo <- factor(c(rep(1,10),rep(2,10)),labels=c("macho","hembra"))

sexo

mandibulas <- data.frame(sexo=sexo,longitud=c(120,107,110,116,114,111,113,117,114,112,110,111,107,108,110,105,107,106,111,111))

mandibulas

mandi.machos <- mandibulas[,2][mandibulas[,1]=="macho"]

mandi.hembras <- mandibulas[,2][mandibulas[,1]=="hembra"]

Para realizar las comparaciones de médias nos basaremos en el tes estadístico t-Student.

t.test(mandi.machos,mandi.hembras)

Welch Two Sample t-test

data: mandi.machos and mandi.hembras

t = 3.4843, df = 14.894, p-value = 0.00336

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

1.861895 7.738105

sample estimates:

mean of x mean of y

Ejemplos del test t-student donde si no hay igualdad de varianzas el mismo realiza el test adecuado (Welch) o si solo le pasas un vector entonces hace test

#Le pasamos una fórmula=>Sexo es factor expicativo de longitud sería como un modelo linea

t.test(formula=longitud~sexo, data=mandibulas)

#test unilateral

t.test(formula=longitud~sexo, data=mandibulas, alternative="greater")

Welch Two Sample t-test

data: longitud by sexo

t = 3.4843, df = 14.894, p-value = 0.00168

alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:

2.383870 Inf

sample estimates:

mean in group macho mean in group hembra

113.4 108.6

#Fijamos la igualdad de varianzas=> t-Student

t.test(formula=longitud~sexo, data=mandibulas, var.equal=T)

Two Sample t-test

data: longitud by sexo

t = 3.4843, df = 18, p-value = 0.002647

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

1.905773 7.694227

sample estimates:

mean in group macho mean in group hembra

t.test(mandibulas[,"longitud"][mandibulas[,"sexo"]=="macho"],mandibulas[,"longitud"][mandibulas[,"sexo"]=="hembra"])

t.test(formula=longitud~sexo, data=mandibulas, alternative="greater")$statistic

t

3.484324

El objeto que devuelve este objeto es htest donde podemos encontrar el valor del estadístico o del p-valor. Se crea una función muy sencilla para crear este test que llamamos tStat.

tStat <- function(...){

return (t.test (...)$statistic)

}

Posible manera de implementar el test Boostrap respecto a la media mediante el tStat, en definitiva, t.test

tStat.boot <- function(x, y, nboot=999, ...){

tObs =tStat(x, y, ...)

mpool = mean(c(x,y))

x1 = x - mean(x) + mpool

x2 = y - mean(y) + mpool

count = 0

for (iboot in 1:nboot) {

if (tStat (sample (x1,replace=T), sample(x2,replace=T),...) >= tObs)

count = count + 1

}

return ((count+1)/(nboot+1))

}

Bajo la hipótesis nula, las restricciones deben cumplir que las dos medias (x,y) tienen que ajustarse al hecho que las medias son iguales.

H0: m1=m2------à media muestral(X1)= media muestral(X2)

Para ello, creamos una media común y entonces reescalamos los valores de cada x e y , de tal forma que la media tanto de x como de y será la común.

x1 = x - mean(x) + mpool

x2 = y - mean(y) + mpool

Hacemos la estimación mediante boostrap para ver como es el p-valor.

tStat.boot(mandi.machos,mandi.hembras)

[1] 0.003

Veamos ahora el test boostrap para comparación de medias de datos apareados. Hemos de dejar la estructura fija.

tStat.boot.paired <- function(x, y, nboot=999, ...){

tObs = tStat(x, y, ...)

mpool = mean(c(x,y))

x1 = x - mean(x) + mpool

x2 = y - mean(y) + mpool

count = 0

n=length(x1)

for (iboot in 1:nboot) {

indices<-sample(1:n,replace=T)

if (tStat (x1[indices], x2[indices],...) >= tObs)

count = count + 1

}

return ((count+1)/(nboot+1))

}

fertilized.cross<-c(92,0,72,80,57,76,81,67,50,77,90,72,81,88,0)

fertilized.self<-c(43,67,64,64,51,53,53,26,36,48,34,48,6,28,48)

tStat.boot.paired(fertilized.cross, fertilized.self , paired=TRUE)

Los datos son del tamaño de la mandibula de chacales, 10 para machos y 10 para hembras.

Ponemos en un vector los datos los 10 primeros machos y depués las hembras.

canis<-c(120,107,110,116,114,111,113,117,114,112, 110,111,107,108,110,105,107,106,111,111)

Para generar permutaciones de los datos simplemente le pedimos una muestra que al no pedirlo con reemplazamiento la muestra permuta simplemente.

sample(canis)

[1] 110 106 110 117 120 113 111 108 112 114 107 107 111 114 111 116 107 111 110

[20] 105

Para comparar las medias simplemente deberemos calcular el estadístico t(ver página 5), pero es demasiado complicado teniendo en cuenta que el denominador es común para todos los cálculos. Luego la primera transformación sería.

t’=media muestral1 –media muestral 2

Pero realmente, simplemente si tomamos una suma de uno de los grupos ya nos dará lo que buscamos(al ser balanceado).

Intenemos generar las perturbaciones pero con unos índices, ya que será útil para el cálculo posterior.

indices<-1:length(canis)

canis[sample(indices)]

igualmente podemos hacer

canis[sample(1:length(canis))]

Creamos una función que crea la diferencia de medias pero con el test bilateral lo que obliga a poner en valor absoluto.

dif.medias<-function(datos,n1){

return(abs(mean(datos[1:n1])- mean(datos[(n1+1): length(datos)])))

}

Versión mejorada:

dif.medias<-function(datos,n1){

n<- length(datos)

if(n1>( n-1))stop("el grupo 1 no puede superar el total")

return(abs(mean(datos[1:n1])- mean(datos[(n1+1):n])))

}

Lo probamos (estimación ^t para toda la muestra)

dif.medias(canis,10)

[1] 4.8

Generamos ahora datos aleatorios que calcularan el estadístico de las permutaciones (nperm=999)

perm.medias<-function(datos, n1,nperm=999){

t0<- dif.medias(datos,n1)

n<- length(datos)

count<-0

for( i in 1:nperm)

{

count<- count +(dif.medias(sample(datos),n1)>= t0)

}

return(t0,count, (count+1)/(nperm+1))

}

Versión 2

perm.medias<-function(datos, n1,nperm=999)

{

t0<- dif.medias(datos,n1)

count<-0

for( i in 1:nperm)

{

count<- count +ifelse(dif.medias(sample(datos),n1)>= t0,1,0)

}

return(t0,count, (count+1)/(nperm+1))

}

Debemos recordar que (count+1)/(nperm+1) nos da una estimación insesgada del estimador pero que nos es más útil.

perm.medias(canis,10)

$t0

[1] 4.8

$count

[1] 3

[[3]]

[1] 0.004

Como en el ejemplo se obtiene un p-valor=0.0013 creemos que esto debe ser debido a que estudia de forma unilateral el test.

dif.medias<-function(datos,n1){

return(mean(datos[1:n1])- mean(datos[(n1+1): length(datos)]))

}

> perm.medias(canis,10)

$t0

[1] 4.8

$count

[1] 1

[[3]]

[1] 0.002

-CASO 2-

Estudiemos ahora para datos apareados. El ejemplo del punto 4.1. [Deberes] Es muy parecido al anterior pero para datos apareados.

-CASO 3-

Veamos como funciona con la regresión, de la página 93 (tabla 6.1). Pensemos que nuestra variable independiente es 1/altitude y la variable dependiente es Hk 1.00 frequency y evaluamos con el coeficiente de correlación de pearson. En realidad nos bastaría la suma de los productos de xy.

Para permutar la solución es dejar una de las variables fijas y permutar la columna de valores de la otra.

colonia<-matrix(c(2,1.25,1.75,1.82,2.63,1.08,2.08,1.59,0.67,0.57,0.5,0.24,0.4,0.5,0.15,0.13,0.11,0.10,98,36,72,67,82,72,65,1,40,39,9,19,42,37,16,4,1,4),ncol=2)

#Hacemos permutaciones de una columna prefijada por nosotros

perm <-function(datos,j)

{

x<- datos[,-j]

y<-datos[,j]

return(sum(x* sample(y)))

}

#Generamos el estadístico por permutación de la suma de x por y.

perm.correl<-function(datos,j=2,nperm=999)

{

t0<-sum(datos[,-j]* datos[,j])

count<-0

for(i in 1:nperm)

{

count<- count +ifelse(perm(datos,j)>= t0,1,0)

}

return(t0,count, (count+1)/(nperm+1))

}

> perm.correl(colonia,2)

$t0

[1] 1015.97

$count

[1] 0

[[3]]

[1] 0.001

-CASO 4-

A Métodos de Montecarlo

A Bioestadística

Volver al índice