## martes, 30 de abril de 2013

### Knitr en R

A Minimal Example for Markdown <!-- Styles for R syntax highlighter

# A Minimal Example for Markdown

This is a minimal example of using knitr to produce an HTML page from Markdown.

## R code chunks

# set global chunk options: images will be 7x5 inches
opts_chunk$set(fig.width = 7, fig.height = 5)  Now we write some code chunks in this markdown file: x <- 1 + 1 # a simple calculator set.seed(123) rnorm(5) # boring random numbers  ## [1] -0.56048 -0.23018 1.55871 0.07051 0.12929  We can also produce plots: par(mar = c(4, 4, 0.1, 0.1)) with(mtcars, { plot(mpg ~ hp, pch = 20, col = "darkgray") lines(lowess(hp, mpg)) })  ## Inline code Inline R code is also supported, e.g. the value of x is 2, and 2 × π = 6.2832. ## Math LaTeX math as usual: $$f(\alpha, \beta) \propto x^{\alpha-1}(1-x)^{\beta-1}$$. ## Misc You can indent code chunks so they can nest within other environments such as lists. 1. the area of a circle with radius x pi * x^2  ## [1] 12.57  2. OK, that is great To compile me, use library(knitr) knit("knitr-minimal.Rmd")  ## Conclusion Markdown is super easy to write. Go to knitr homepage for details . --> Read more... ## martes, 21 de agosto de 2012 ### Clase 6: Índices de biodiversidad en R-project ## Índices de biodiversidad Continuando con este breve curso sobre bioestadística en R (Clase 1, Clase 2, Clase 3, Clase 4 y Clase 5), en esta oportunidad intentaremos abordar el problema de los índices de biodiversidad utilizando R. Comenzaremos con algunos conceptos previos. #### ¿Qué queremos decir por biodiversidad? La biodiversidad o diversidad biológica se define como “la variabilidad entre los organismos vivientes de todas las fuentes, incluyendo, entre otros, los organismos terrestres, marinos y de otros ecosistemas acuáticos, así como los complejos ecológicos de los que forman parte; esto incluye diversidad dentro de las especies, entre especies y de ecosistemas” (UNEP, 1992). El término comprende, por tanto, diferentes escalas biológicas: desde la variabilidad en el contenido genético de los individuos y las poblaciones, el conjunto de especies que integran grupos funcionales y comunidades completas, hasta el conjunto de comunidades de un paisaje o región (Solbrig, 1991; Halffter y Ezcurra, 1992; Heywood, 1994; UNEP, 1992; Harper y Hawksworth, 1994). Es decir, existen distintas escalas de diversidad: • Diversidad a escala genética. • Diversidad de especies: la diversidad alfa es la riqueza de especies de una comunidad particular a la que consideramos homogénea. • Diversidad a nivel de comunidades: la diversidad beta es el grado de cambio o reemplazo en la composición de especies entre diferentes comunidades en un paisaje, y la diversidad gamma es la riqueza de especies del conjunto de comunidades que integran un paisaje, resultante tanto de las diversidades alfa como de las diversidades beta (Whittaker, 1972). ### Diversidad de especies. La diversidad alfa tiene dos componentes fundamentales: • riqueza específica: número de especies que tiene un ecosistema • estructura (equitabilidad): mide la distribución de la abundancia de las especies, es decir, cómo de uniforme es un ecosistema. #### ¿Cómo podemos traducir un concepto intuitivo en algo que podamos medir? Fuente: http://entomologia.rediris.es/sea/manytes/metodos.pdf #### Métodos paramétricos: índices La medida de ocurrencia Jaccard o Coeficiente de comunidad es una medida muy simple e intuitiva (Shi 1993; Magurran 2004). La medida de Sorenson o índice de Dice, Czekanowski or Coincidencem hace más énfasis a las especies presentes compartidas más que a las no compartidas, y su cálculo es relativamente simple e intuitivo (Wolda 1981; Hub?lek 1982). Otros dos índices similares sue se utilizan ocasionalmente son las medidas de Ochiai y Kulczynski. Sin embargo, se eligen las medidas de Jaccard o Soreson porque son las más utilizadas actualmente. El índice de Simpson intenta solucionar uno de los problemas más típicos, el de diferentes tamaños muestrales entre localidades. Comprar sitios con intensidades desiguales puede generar visiones sesgadas. Teniendo esto en cuenta, Simpson (1960) desarrolló una medida para tamaños muestrales distintos. Su fórmula reescala el valor por el número de especies del último sitio muestreo, por lo que una submuestra tendrá similitud con la muestra original. Por tanto, este índice es utilizado con datos altamente variables en intensidad de muestreo, por ejemplo los datos de fósiles. La medida de Bray-Curtis es posiblemente la más utilizada, debido a su gran relación con las distancias ecológicas bajo condiciones variadas (Bray and Curtis 1957; Faith etal. 1987; Minchin 1987; Clarke 1993). Esta medida es equivalente a la de Soreson, cuando la utilizamos como medida de similitud con datos de ocurrencia. A pesar de ser menos conocido, el índice de Mosisita-Horn es recomendable debido a su relativa independencia del tamaño muestral y diversidad (Wolda 1981; Magurran 2004). Existen muchas variantes de este índice. #### Métodos no paramétricos y de rarefacción. Uno de los métodos típicos para solucionar los problemas de muestreo desigual en intensidad, es el índice de rareza (o también conocido con el peculiar nombre de "rarefacción" o interpolación de los datos; Sanders 1968). Los índices de rareza permiten comparar entre distintas comunidades donde cada una es "enrarecida" a un número igual de ejemplares muestreados (Heck et al. 1975; Foote 1992; Colwell and Coddington 1994). Los índices de rareza interpolan nuevos datos, mientras que los estimadores no paramétricos de especies extrapolan los datos para encontrar el "verdadero" número de especies (Colwell and Coddington 1994). Los estimadores no paramétricos utilizan el número de especies raras para encontrar la manera de calcular qué tan probable es que hayan más especies sin encontrar en la muestra. Por ejemplo, el estimador de Chao (Chao 1984; Colwell and Coddington 1994) calcula la diversidad de especies verdadera estimada de la muestra mediante los datos de abundancia. En caso de que no contemos con datos de abundancia sino de ocurrencia, podremos utilizer el estimator Chao 2 (Chao 1987; Colwell and Coddington 1994), que estima la diversidad de las especies a partir de los datos de ocurrencia de múltiples muestras. Otro estimador es el de Jackknife (Burnham & Overton, 1978, 1979) originalmente utilizados para dates de capture y recaptura. Dentro del paquete fossil existe un método de rarefacción conocido como la curva de Coleman(Coleman 1981; Coleman et al. 1982), que se desarrolla como un método de remuestreo más que mediante una fórmula de rareza. Por último está el estimator de riqueza boostrap (Smith & van Belle, 1984). #### Aplicación en R: paquetes fossil y vegan   ########################################################################### # Análisis de datos con R # # 2012 # # Msc. Rosana Ferrero # # http://statisticalecology.blogspot.com/ # ###########################################################################  ###########################################################################  # 6. ANÁLISIS DE BIODIVERSIDAD # # librerías fossil y vegan ###########################################################################  ########################################################################### ############ con la librería fossil ################# ############ 1. Manipulación de los datos ############ 2.  ########################################################################### library(fossil) data(fdata.list) head(fdata.list)  ############ 1. Manipulación de los datos #las funciones estimadoras de especies pueden agruparse en 2 categorías: las que usan los datos de ocurrencia y aquellas que usan los datos de abundancia. #matriz de ocurrencias ocurrencias=create.matrix(fdata.list,tax.name="species",locality="locality")  #toma la lista de las especies y su ocurrencias, y las convierte en una matriz de especies (filas) y localidades (columnas).   #matriz de abundancia abundancias=create.matrix(fdata.list,tax.name="species",locality="locality", abund=TRUE,abund.col="abundance")  #crea una matriz de abundancias, se trata de la misma función anterior pero ahora incluye la opción abund=TRUE y el argumento abund.col que corresponde a la columna de los valores de abundancia. #para este ejemplo, el resultado es el mismo que data(fdaa.mat)  #matriz de coordenadas coordenadas=create.lats(fdata.list,loc="locality",long="longitude", lat="latitude")  #extrae las coordenadas de los sitios, eliminando aquellas entradas que se repiten.  #La función dino.dist() toma una matrix de ocurrencia de especies vs. localidad (o grupos) y devuelve una matriz de distancia pareadas. ##matriz de distancia entre cualquier número de localidades #ecol.dist(x, method = sorenson, type = "dis") #dino.dist(x, method = sorenson, type = "dis") data(fdata.mat) ecol.dist(fdata.mat) ecol.dist(fdata.mat,simpson,"sim")  #gráfico data(fdata.mat) fdata.dist=dino.dist(fdata.mat) fdata.mst=dino.mst(fdata.dist) data(fdata.lats) library(maps) map("state") mstlines(fdata.mst,coordinates(fdata.lats)) points(coordinates(fdata.lats),pch=16,col="white",cex=3) points(coordinates(fdata.lats),pch=1,cex=3) text(coordinates(fdata.lats),labels=LETTERS[1:12])   ############ 2. Índices de distancia/similitud/diversidad beta  #Coefficient Name Formulae Alternate Name Function Call #Jaccard a/(a+b+c) Coefficient of Community jaccard() #Sorenson 2a/(2a+b+c) Dice, Czekanowski, Coincidence Inde sorenson() #Simpson a/(a+min(b,c)) - simpson() #Braun-Blanquet a/(a+max(b,c)) - braun.blanquet() #Ochiai a/sqrt{(a+b)(a+c)} Coefficient of Closeness ochiai() #Kulczynski [a/(a+b)+a/(a+c)]/2 - kulczynski() #otros: braun.blanquet(x, y) bray.curtis(x, y) euclidean(x,y) kulczynski(x,y) jaccard(x, y) manhattan(x, y) morisita.horn(x, y) ochiai(x, y) simpson(x, y) sorenson(x, y) #Todas las funciones de similitud se utilizan de la misma manera. Las funciones necesitan dos argumentos que representen las dos muestras. Es importante que las ocurrencias de especies estén colocadas de misma manera para cada sitio, y que cada especie ausente se represente por el cero.  sampleA=c(1,1,0,1,1,1,1) sampleB=c(0,1,1,0,0,1,1) sorenson(sampleA,sampleB) ##Conteo de especies de 2 localizaciones diferentes. a = c(1,0,4,3,5,0,0,7) b = c(2,1,3,0,0,1,0,6) bray.curtis(a,b) jaccard(a,b) simpson(a,b) sorenson(a,b) morisita.horn(a,b)    ##### 3. Métodos con muestras complejas (incompletas) o pequeñas: aleatorización, remuestreo, boostrap.  data(fdata.mat) head(fdata.mat) #es nuestra matriz de "abundancias" chao1(fdata.mat) jack1(fdata.mat) ## con aleatorización #spp.est(). La función tiene varias opciones: el número de aleatorizaciones y si utilizar o no datos de abundancias. Calcula la curva de rareza de Chao, Coverage Estimators y Jacknife, así como las desviaciones estándar para todas las estimaciones. Por defecto la función corre 10 aleatorizaciones de los datos, sin embargo si queremos mejores estimaciones tendremos que aumentar dicho número. e run. #ejemplo con datos de abundancia data(fdata.mat) spp.est(fdata.mat, abund = TRUE, counter = FALSE) #ejemplo con datos de ocurrencia data(fdata.mat) spp.est(fdata.mat, abund = FALSE, counter = FALSE) ## los estimators de Chao-Jaccard y Chao-Sorenson como estimadores de las especies cuando tenemos datos incompletos  #chao.sorenson(x, y) chao.jaccard(x, y) #se deben datos dos vectores separados, con las especies dispuestas con el mismo orden, desde el área A a la B. Si la especie está presente en uno de los sitios pero no en el otro, debe ser identificada en ambos sitio; en el sitio donde no se encontró se codificará con un cero. Especies que no estén presentes en ninguno de los dos sitios serán ignoradas.  ##conteo de especies en 2 localizaciones diferentes  a = c(1,0,4,3,5,0,0,7) b = c(2,1,3,0,0,1,0,6) chao.sorenson(a,b) chao.jaccard(a,b) ## Método de Bootstrap en estimadores de riqueza de especies, para datos de abundancia o para datos de presencia-ausencia.  ## vector muestra a=c(0,5,1,1,2,0,0,1,0,0,8,45) bootstrap(a,samples=45) ## formato matricial a=matrix(c(0,5,1,1,2,0,0,1,0,0,8,45),4,3) bootstrap(a) bootstrap(a,,FALSE) ## matriz de presencia-ausencia a=matrix(c(0,1,1,1,1,0,0,1,0,0,1,1),4,3) bootstrap(a,,FALSE) ## matiz de abundancia bootstrap(fdata.mat,,TRUE,samples=1000)        ########################################################################### ############## con la librería vegan ################### ########################################################################### #Vamos a utilizar unos datos clásicos de la diversidad de árboles de Barro Colorado. Los datos consisten de 50 parcelas (filas) de 1 hectárea donde se realiza el conteo de los árboles en cada parcela y la identificación de su especie (columnas). Se incluyeron todos los árboles con diámetros por encima de 10cm. Se encontraron 225 especies. data(BCI) # Barro Colorado Island tree data # 50 1-ha plots, counts of 225 spp. View(BCI) #para ver los datos fix(BCI) ### Riqueza específica: índices cts= c(20,12,5,3,2,1,1,0,0,0) # counts for 10 spp #1. Sin utilizar vegan sum(cts>0) #7: sums a logical vector: T=1,F=0 #2. Utilizando vegan specnumber(BCI) # función para encontrar el número de especies, trabaja con matrices specnumber(cts) # error specnumber(matrix(cts,nrow=1)) #Diversidad: otras métricas. Función diversity() diversity(BCI, index="shannon") # Shannon-Wiener, H diversity(BCI, index="simpson") # Simpson, 1-D diversity(BCI, index="invsimpson") # Inverse Simpson, 1/D fisher.alpha(BCI) # parameter alpha of the log-series Vegan no tiene índices de equitatividad o igualdad, peru podemos calcular el más común de ellos, el índice de Pielou: J = H′/ log(S) de una manera sencilla J <- H/log(specnumber(BCI))  # Rarefacción: # min N in BCI data is 340 bci.340 = rarefy(BCI, sample=340, se=FALSE) plot(bci.340) #Curvas de rarefacción. #? Based on rarefying to many different subsample sizes. #? rarefy() can only take a single value for sample. #? Can call it lots of times, and collect the data: tedious! #? In Session 4, we?ll make our own function to do this.  #### También existen las funciones: specaccum #Species accumulation renyiaccum; tsallisaccum #diversity accumulation fisherfit; prestonfit; radfit #Species abundance models specpool, estimateR # Extrapolated species richness beta diver; betadisper #Indices of betadiversity taxondive #Taxonomic diversity treedive #one option of functional diversity Null models and nestedness (oecosimu with support functions) Technical and advanced functions like dissimilarity indices (vegdist, designdist) and their extensions (step across) #### Bibliografía y material para profundizar: http://entomologia.rediris.es/sea/manytes/metodos.pdf http://duncanjg.files.wordpress.com/2008/02/clasediversidad1.pdf http://www.um.es/docencia/geobotanica/ficheros/practica2.pdf http://prof.usb.ve/eklein/labeco2/biblio/Estructura%20de%20las%20comunidades.pdf http://palaeo-electronica.org/2011_1/238/238.pdf http://cran.r-project.org/web/packages/vegan/vignettes/diversity-vegan.pdf http://www.bio.utk.edu/fesin/msa2009/snowbeam4.pdf Read more... ### Clase 5: Remuestreo en R-project ## Remuestreo Los métodos de remuestreo o submuestreo (i.e. resampling) evalúan los estadísticos en remuestras o submuestras obtenidas a partir de los datos originales, y mediante estos valores obtienen estimadores de lmedidas de exactitud o de la distribución muestral del estadístico. El remuestreo son métodos que permiten: • Estimar la precisión de muestras estadísticas (según su media, mediana, varinza, etc..). • Intercambiar marcadores de puntos de datos al realizar test de significación (test de permutación, tambien denominados tests exactos, tests de aleatoriedad, o pruebas de re-aleatoriedad) • Validar modelos para el uso de subconjuntos aleatorios (bootstrapping, validación cruzada) #### Tipos de remuestreo. Existen al menos 4 grupos principales de remuestreo: 1. aleatorización de pruebas exactas (randomization exact test) 2. validación cruzada (cross-validation): se usa para evaluar los resultados de un análisis estadístico y garantizar que son independientes de la partición entre datos de entrenamiento y prueba. Definida una muestra, se realizan particiones de las mismas y se evalúan los estadísticos para cada una de las particiones. El caso más sencillo consiste en tomar 2 particiones, una de entrenamiento y la otra de prueba. 3. jackknife: para la estimación del sesgo y del error estándar de un estadístico. El nombre alude a las navajas, y más precisamente a aquellas multiuso, que sirven pa tó. Definida una muestra de observaciones de tamaño n, el método jacknife consiste en suprimir cada vez un conjunto de observaciones g y calcular sobre el conjunto (n-g) restante de datos el estadístico de interés. Una vez obtenidas las n estimaciones del estadístico para cada una de las muestras Jackknife, se puede calcular una estimación del sesgo asociado a la muestra original, así como una estimación no paramétrica del error estándar del estadístico. 4. boostrap: su nombre refiere a su autosuficiencia, ya que significa “levantarse tirando hacia arriba de las propias correas de las botas”. El Bootstrap se basa en el remuestreo con reposición a partir de la muestra original para obtener nuevas muestras, calcular en cada una de ellas el estadístico deseado, y así obtener una distribución de valores para el estadístico en la que podemos calcular su dispersión (análogo del error estándar) y determinar unos límites de confianza utilizando esa distribución. #### Descripción y aplicación en R #### 1) Aleatorización de pruebas exactas o test de permutación (Randomization exact test or permutation test):    NOTA: Also known as the permutation test, this test was developed by R. A. Fisher (1935/1960), the founder of classical statistical testing. However, in his later years Fisher lost interest in the permutation method because there were no computers in his days to automate such a laborious method. El propósito es simular artificialmente la aleatoriedad. Los datos se re-asignan aleatoriamente (datos permutados) de tal manera de obtener p-valores basados en estos nuevos datos. EJ: evaluar la diferencia entre dos m?todos. Luego de evaluar cada conjunto de datos permutados, por ejemplo 100, el investigador puede poner juntos los t-valores y obtener un "p-valor exacto" (la probabilidad de que la diferencia solo ocurra por azar). Por ejemplo, si el t-valor de 1.55 ocurre solo 5 veces en las 100, el p-valor exacto ser?a .05. Problemas: d?bil ante tama?os de muestras peque?os. #### 2) Cross-validation:    NOTA: Simple cross-validation was proposed by Kurtz (1948) as a remedy for the Rorschach test, a form of personality test which was criticized by psychometricians for its lack of common psychometric attributes such as data normality. Based on Kurtz's simple cross-validation, Mosier (1951) developed double cross-validation, which was later extended to multicross-validation by Krus and Fuller (1982).  El propósito es averiguar si los resultados son replicables o solo una cuesti?n de fluctuaciones aleatorias, es decir, verificar la replicabilidad de los resultados. La muestra se divide aleatoriamente en 2 o más subconjuntos y los resultados de la prueba sn validados comparando las sub-muestras. Varias categor?as: simple, doble o m?ltiple: *Simple cross-validation. Take regression as an example. In the process of implementing a simple cross-validation, the first sub-sample is usually used for deriving the regression equation while another sub-sample is used for generating predicted scores from the first regression equation. Next, the cross-validity coefficient is computed by correlating the predicted scores and the observed scores on the outcome variable. *Double cross-validation. Double cross-validation is a step further than its simple counterpart. Take regression as an example again. In double cross-validation regression equations are generated in both sub-samples, and then both equations are used to generate predicted scores and cross-validity coefficients. *Multicross-validation. Multicross-validation is an extension of double cross-validation. In this form of cross-validation, double cross-validation procedures are repeated many times by randomly selecting sub-samples from the data set. In the context of regression analysis, beta weights computed in each sub-sample are used to predict the outcome variable in the corresponding sub-sample. Next, the observed and predicted scores of the outcome variable in each sub-sample are used to compute the cross validated coefficient. Problemas: d?bil ante tama?os de muestras peque?os. Ang (1998) argued that cross-validation is problematic because splitting an already small sample increases the risk that the beta weights are just artifacts of the sub-sample. Thus, Ang endorsed use of Jackknife, which will be discussed next. Es ?til cuando tenemos que elegir entre modelos (por ejemplo, anidados) ya que la validaci?n cruzada nos permite elegir el modelo que "mejor predice" los valores, a diferencia de los m?todos de m?nimos cuadrados o m?xima verosimilitud que eligen los par?metros dentro de un conjunto de modelos que "mejor ajustan" los valores. NOTA: a pesar de que la validaci?n cruzada es ?til para elegir entre modelos, no tenemos idea de si la diferencia entre ellos es estad?sticamente significativa. La manera m?s f?cil de hacer validaci?n cruzada es dividir los datos en dos partes: el conjunto de entrenamiento (training set) que ser? usado para ajustar varios modelos a los datos, y el conjunto de prueba (test set) que se usar? para elegir entre modelos. un inconveniente del m?todo anterior es que reduce el poder de los datos para ajustar un modelo ocultando algunas observaciones para validaci?n. Podemos minimizar esta p?rdida tomando un conjunto de entrenamiento tan grande como sea posible, pero en tal caso no tendr?amos suficientes puntos para la validaci?n. Evitamos esto utilizando una estratgia "leave one out", superficialmente similar al jackknife. (el jackknife trabaja con un modelo particular paara estimar un estimador de sesgo y varianza muestral, mientras que en la validaci?n cruzada estamos eligiendo entre modelos). Ejemplo: an?lisis de regersi?n ? Leaving out the ith pair of points (Xi, Yi), use least squares methods to find parameter estimates ?a-i and ?b-i. ? Hence, find the best predictor of the left out Yi: ?Y-i = ?a-i + ?b-i*Xi ? Repeat the above for each point and calculate the cross validation score, or estimate predictive error (1/n)*sum(Yi-?Y-i)^2 ? If we also calculate the cross validation score for the log(X) model we can again use predictive error as a basis to choose between the models. #### 3) Jackknife: Also known as the Quenouille-Tukey Jackknife, this tool was invented by Maurice Quenouille (1949) and later developed by John W. Tukey (1958). As the father of EDA, John Tukey attempted to use Jackknife to explore how a model is influenced by subsets of observations when outliers are present. The name "Jackknife" was coined by Tukey to imply that the method is an all-purpose statistical tool. El prop?sito es detectar valores extremos y c?mo cada sub-muestra afecta al modelo. Se repite el mismo test dejando un sujetopor vez ("leave one out"). Este procedimiento es especialmente ?til cuando la dispersi?n de la distribuci?n es amplia o existen valores extremos en el conjunto de datos. En estos casos se espera que el Jackknife retorne una estimaci?n de menor sesgo. Se lo puede ver como una versi?n m?s organizada o sistem?tica del bootstrap. Algoritmo: Given a data set of n observations we create n new data sets, each of size n - 1, by leaving out one point at a time. ? Suppose we have a sample of observations X = {X1, . . . Xn} that we will use to estimate a parameter Theta. ? Suppose that T = T(X) is an estimator of Theta. ? Draw n jackknife samples X-1, X-2, ...., X-n by leaving out the first, second . . . nth observation respectively. ? Calculate T-i=T(X-i) for each i. ? Calculate T* =(1/n)*sum(T-i) ? The jackknife estimate of variance is ((n-1)/n)*sum(T-i-T*)^2 ? And the jackknife estimate of bias is (n - 1)(T*- T) which leads to a bias corrected estimate of nT - (n - 1)T* when we subtract this from T. Ejemplo: dat=c(3.9140011, 0.1819694, 2.5604846, 0.3703330, 6.7950865, 4.0200424, 0.660441) # queremos estimar su media. jk = rep(0,length(dat)) for (i in 1:length(jk)) jk[i] = mean(dat[-i]) mean(jk) #[1] 2.643194 (length(jk)-1)/length(jk) * sum((jk-mean(jk))^2) # [1] 0.855013 mean(dat) #[1] 2.643194, observemos que la media de la muestra jackknife es igual a la media de los datos originales. var(dat)/length(dat) #[1] 0.855013, observemos que la estimaci?n usual de la varianza (insesgada) de la media muestra es la misma que la del estimados jackknife de la varianza. Sin embargo, la estimaci?n bootstrap de la varianza de ma media muestral conveerge al MLE (sesgado). # f?rmula general de jackknife jackknife=function(T,x) { n = length(x) jk = rep(0,n) for (i in 1:n) jk[i] = T(x[-i]) Tstar = mean(jk) jkvar = (n-1)/n * sum((jk-Tstar)^2) list(var = jkvar, bias = (n-1)*(Tstar - T(x)), corrected = n*T(x)-(n-1)*Tstar) } salida=jackknife(T=mean,x=dat) # la estimaci?n de la varianza por jackknife puede utilizarse para realizar intervalos de confianza y test de hip?tesis, tal como se realiza con la manera usual de estimar la varianza. #### 4) Bootstrap: This technique was invented by Bradley Efron (1979, 1981, 1982) and further developed by Efron and Tibshirani (1993). "Bootstrap" means that one available sample gives rise to many others by resampling (a concept reminiscent of pulling yourself up by your own bootstrap). While the original objective of cross-validation is to verify replicability of results and that of Jackknife is to detect outliers, Efron (1981, 1982) developed bootstrap with inferential purposes. El prop?sito es la inferencia. Genera resultados con menor sesgo y m?s consistentes, respecto al Jackknife (Fan & Wang, 1996). Comparado con la t?cnica de Jackknife, la estrategia de remuestreo del Bootstrao es m?s rigurosa en t?rminos de la magnitud de la replicaci?n. En Jackknife, el n?mero de remuestreos est? limitado por el n?mero de observaciones (n-1), pero en bootstrap, la muestra original puede ser duplicada tantas veces como que quiera y por tanto esta muestra expandida se trata como una poblaci?n virtual. Las muestras se obtienen de esta poblaci?n para veerificar los estimadores. A diferencia de la validaci?n cruzada o el Jackknife, el Bootstrap emplea un muestreo con reemplazamiento, que -en t?rminos de simulaci?n de aleatoriedad- es m?s adecuado que el muestreo sin reemplazo. A diferencia de la validaci?n cruzada o el Jackknife, donde el n de las submuestras es menor que el de la muestra original, en Bootstrap cada remuestreo tiene el mismo n?mero de observacinones que la muestra original. Por tanto, tiene la ventaja de modelar los impactos del tama?o real de la muestra (Fan & Wang, 1996). Algoritmo: ? Suppose we have a sample of observations X = {X1, . . . Xn} that we want to use to estimate a parameter Theta. ? Suppose also that we have derived T(X), an estimator of Theta, and want to evaluate its sampling variance. ? Draw a bootstrap sample X1 of size n by sampling with replacement from X. ? Calculate T1=T(X1) ? Repeat the resampling b times to get T1, T2, ..., Tb ? Calculate the mean of the bootstrapped statistics meanT=(1/b)*sum(Ti) ? Finally, calculate the bootstrap estimate of variance (1/(b-1))*sum(Ti-meanT)^2 Ejemplo: dat=c(3.9140011, 0.1819694, 2.5604846, 0.3703330, 6.7950865, 4.0200424, 0.660441) #sea el conjunto de datos dat, queremos usarlo para estimar la media de la poblaci?n de la que proviene mean(dat) #para evaluar la varianza de esta estimaci?n, podemos comenzar tomando una muestra bootstrap y calculando la media en ella. boot1=sample(dat, replace=TRUE) boot1 mean(boot1) # esto nos da T1=2.461595 # si lo ponemos en un loop podemos generar una muestra completa de bootstrap y evaluar la varianza. bs=rep(0,100000) for (i in 1:length(bs)) bs[i] = mean(sample(dat,replace=TRUE)) var(bs) # mi estimaci?n de la varianza muestra del estimador es 0.7376547 sum((dat-mean(dat))^2) / length(dat) / length(dat) # la comparamos con la forma usual de obtener la varianza muestral #cuantas muestras tomar depende en qu? tan r?pido queramos calcular nuestras estimaciones. # f?rmula general de bootsrap bootstrap=function(T,x,b=10000) { bs = rep(0,b) for (i in 1:b) bs[i] = T(sample(x,replace=TRUE)) bs } salida=bootstrap(T=mean,x=dat,b=10000) #mi estimaci?n de la varianza muestra del estimador es 0.7312076 # la estimaci?n de la varianza por bootstrap puede utilizarse para realizar intervalos de confianza y test de hip?tesis, tal como se realiza con la manera usual de estimar la varianza. ## Variantes del bootstrap # i) smoothed bootstrap. Se agrega a cada remuestreo una cantidad comparativamente peque?a de valores generados aleatoriamente. Esto tiene el mismo efecto que ajustar una estimaci?n del n?cle (kernel) de densidad a los datos, y remuestrar en ellos. # ii) parametric bootstrap. Se ajusta una densidad param?trica a los datos observados y se realiza el remuestreo en ella. # iii) wild bootstrap. Para datos distribuidos sim?tricamente respecto a 0. Cada valor remuestrado se multiplica por -1 con una probabilidad 0.5 antes de ser usado. # iv) balanced bootstrap. concatenamos b copias del conjunto original de datos en un vector y permutamso aleatoriamente en ?l. Los primeros n elementos se vuelven la primer muestra bootstrap, los degundos n elementos la segunda muestra y as? sucesivamente. #### 6) Markov chain Monte Carlo or MCMC (o Metropolis & Hastings samplers) Dada una distribuci?n phi en un n?mero finito de estados, podemos construir alg?n proceso equivalente a una cadena de Markov irreducible P de tal manera que la distribuci?n erg?dica de la cadena sea phi. mirar tambi?n: http://cran.r-project.org/web/packages/mcmc/vignettes/demo.pdf Mirar en: http://www.iiap.res.in/astrostat/School07/R/MCMC.html A toy example to calculate P(-1 < X < 0) when X is a Normal(0,1) random variable: xs = rnorm(10000) # simulate 10,000 draws from N(0,1) xcount = sum((xs>-1) & (xs<0)) # count number of draws between -1 and 0 xcount/10000 # Monte Carlo estimate of probability pnorm(0)-pnorm(-1) # Compare it to R's answer (cdf at 0) - (cdf at -1) #i) Metropolis sampling Implica una etapa de proposici?n y una etapa de aceptaci?n, de tal manera que para movernos desde el estado i al estado j debemos primero proponer el movimiento y luego aceptar el movimiento. Algoritmo: ? Proposals are generated according to a symmetric transition matrix Q, so that Qi,j = Qj,i ,or the probability of proposing a move to state j from state i is the same as proposing a move to state i from state j. ? A proposed move from state i to state j is accepted with probability min(1,phij/phii) ? If the proposal is accepted we move to state j, and if it is rejected, we stay in state i. para correr esta cadena solo neceistamos definir el mecanismo de propuesta, que no depende de phi, y luego evaluar la propuesta utilizando la proporci?n phij/phii. Note que no necesitariamente tenemos que conocer Qi,j, solo necesitamos conocer que Q sea sim?trica. Note que para muestras v?lidas de phi, la cadena de markov debe ser irreducible. Ejemplo: implementar Metropolis sampler dado un vector de probabilidades phi. f.metro=function(pi,n=10000) { s = 1:length(pi) x = rep(0,n) x[1] = sample(s,1) for (i in 2:n) { if (runif(1) < 0.5) z = x[i-1] - 1 else z = x[i-1] + 1 if (z < 1 | z > length(s)) a = 0 else a = min(1, pi[z]/pi[x[i-1]]) if (runif(1) < a) x[i] = z else x[i] = x[i-1] } x } f.runmc(P) pi= z=f.metro(pi,100000) table(z)/length(z) #z #1 2 3 4 5 6 #0.20470 0.18060 0.04159 0.04179 0.12459 0.40673 # esta salida concuerda bastante bien con el vector phi de entrada. # From Robert & Casells: A Metropolis{Hastings sample is then generated with the following R code: a=2.7; b=6.3; c=2.669 # initial values Nsim=5000 X=rep(runif(1),Nsim) # initialize the chain for (i in 2:Nsim) { Y=runif(1) rho=dbeta(Y,a,b)/dbeta(X[i-1],a,b) X[i]=X[i-1] + (Y-X[i-1])*(runif(1)<rho) } #i) Hasting?s sampling Este m?todo generaliza el m?todo de Metropolis al caso en el cual Q no es sim?trico. Tiene la misma estructura de propuesta y aceptaci?n, pero la probabilidad de aceptaci?n viene dada por max(1,(phij*Qj,i)/(phii*Qi,j)). Para realizar este muestreo necesitamos conocer m?s expl?citamente en detalle la matriz Q, especificamente necesitamos la proporci?n Qi,j/Qj,i. ### M?s ejemplos: http://www.mayin.org/ajayshah/KB/R/documents/boot.html http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-bootstrapping.pdf http://www.statmethods.net/advstats/bootstrapping.html http://cran.r-project.org/doc/Rnews/Rnews_2002-3.pdf ##### Types of resampling: # bibliograf?a: http://pareonline.net/getvn.asp?v=8&n=19 http://episun7.med.utah.edu/~alun/teach/compstats/notes.pdf http://cran.r-project.org/web/packages/mcmc/vignettes/demo.pdf  ### EJEMPLOS ################################################################################# #### Resample data ################################################################################# ## Basado en http://www.zoology.ubc.ca/~schluter/zoo502stats/Rtips.resample.html #Let's assume that the data are a sample of measurements for a single variable stored in a vector "x". The data may be numeric or categorical. #A single bootstrap replicate is obtained as follows. The replace option is used to indicate that sampling is carried out with replacement. xboot = sample(x, replace=TRUE) #Calculate the statistic of interest (for example, the mean) on the resampled data in "xboot" and store the result in a vector created for this purpose. z = vector() # initialize z (do only once) z[1] = mean(xboot) # first bootstrap replicate estimate #Repeat steps (1) and (2) many times. The result will be a large collection of bootstrap replicate estimates for subsequent analysis. #In other cases, two or more variables are measured on individuals (e.g., stem height, leaf area, petal diameter, etc). Assume that each row of a data frame "mydata" is a different individual, and each column a different variable. #To resample individuals (i.e., rows), iboot = sample(1:nrow(mydata), replace=TRUE)) bootdata = mydata[iboot,] #The data frame "bootdata" will contain a single bootstrap replicate including all the variables. #Calculate the statistic of interest on the resampled data and store the result in vector created for this purpose. For example, to calculate the correlation between two variables x and y in "bootdata", z = vector() # initialize z (do only once) z[1] = cor(bootdata$x, bootdata$y) #Repeat steps (1) and (2) many times. The result will be a large collection of bootstrap replicate estimates for subsequent analysis. # Basado en http://www.zoology.ubc.ca/~schluter/zoo502stats/Rtips.resample.html ### Randomization test #A randomization test uses resampling and the computer to generate a null distribution for a test statistic. The test statistic is a measure of association between two variables or difference between groups. Each randomization step involves randomly resampling without replacement the values of one of the two variables in the data. The two variables may be categorical (character or factor), numeric, or one of each. ## Both variables categorical #R has a built-in randomization procedure for a contingency test of association when both of two variables are categorical (call them A1 and A2). To apply it, execute the usual command for the ?2 contingency test, but set the "simulate.p.value" option to TRUE. The number of replicates in the randomization is set by the option "B" (default is 2000). Each randomization rearranges the values in the contingency table while keeping all the row and column totals fixed to their observed values. chisq.test(A1, A2, simulate.p.value = TRUE, B = 5000) ## Data are numeric #If one or both of the variables is numeric, then you will need to create a short loop to carry out the resampling necessary for the randomization test. Choose one of the two variables to resample (call it "x"). It doesn't matter which of the two variables you choose. Keep the other variable (call it "y") unchanged (there is no benefit to resampling both variables). #Resample "x" without replacement to create a new vector (call it x1). x1 = sample(x, replace = FALSE) #Calculate the test statistic to measure association between y and the randomized variable x1. Store the result in a vector created for this purpose. For example, to calculate the correlation between the two variables, z = vector() # initialize z (do only once) z[1] = cor(x1, y) # first randomization result #Repeat steps (1) and (2) many times. The result will be a large collection of replicates representing the null distribution of your test statistic. #Use this null distribution to calculate an approximate P-value for your test. This involves calculating the proportion of values in the null distribution that are as extreme or more extreme than the observed value of the test statistic (calculated from the real data). Whether this calculation must take account of just one or both tails of the null distribution depends on your test statistic. ################################################################################# ### Permutation test: randomization test for small size samples ################################################################################# ## Basado en http://www.burns-stat.com/pages/Tutor/bootstrap_resampling.html # The idea: Permutation tests are restricted to the case where the null hypothesis really is null -- that is, that there is no effect. If changing the order of the data destroys the effect (whatever it is), then a random permutation test can be done. The test checks if the statistic with the actual data is unusual relative to the distribution of the statistic for permuted data. # Our example permutation test is to test volatility clustering of the S&P returns. Below is an R function that computes the statistic for Engle's ARCH test. engle.arch.test = function (x, order=10) { xsq = x^2 nobs = length(x) inds = outer(0:(nobs - order - 1), order:1, "+") xmat = xsq[inds] dim(xmat) = dim(inds) xreg = lm(xsq[-1:-order] ~ xmat) summary(xreg)$r.squared * (nobs - order)
}
# All you need to know is that the function returns the test statistic and that a big value means there is volatility clustering, but here is an explanation of it if you are interested. The test does a regression with the squared returns as the response and some number of lags (most recent previous data) of the squared returns as explanatory variables. (An estimate of the mean of the returns is generally removed first, but this has little impact in practice.) If the last few squared returns have power to predict tomorrow's squared return, then there must be volatility clustering. The tricky part of the function is the line that creates inds. This object is a matrix of the desired indices of the squared returns for the matrix of explanatory variables. Once the explanatory matrix is created, the regression is performed and the desired statistic is returned. The default number of lags to use is 10.
#    A random permutation test compares the value of the test statistic for the original data to the distribution of test statistics when the data are permuted. We do this below for our example.
spx.arch.perm = numeric(1000)
for(i in 1:1000) {
spx.arch.perm[i] = engle.arch.test(sample(spxret))
}
hist(spx.arch.perm, col='yellow')
abline(v=engle.arch.test(spxret), col='blue')
# The simplest way to get a random permutation in R is to put your data vector as the only argument to the sample function. The call:
sample(spxret)
#is equivalent to:
sample(spxret, size=length(spxret), replace=FALSE)

#    A simple calculation for the p-value of the permutation test is to count the number of statistics from permuted data that exceed the statistic for the original data and then divide by the number of permutations performed. In R a succinct form of this calculation is:
mean(spx.arch.perm >= engle.arch.test(spxret))
#    In my case I get 0.016. That is, 16 of the 1000 permutations produced a test statistic larger than the statistic from the real data. A test using 100,000 permutations gave a value of 0.0111.

## Basado en ayuda del paquete BHH2
## Permutation test for means and variance comparisons.
library(BHH2)

# Permutation test for Tomato Data
data(tomato.data)
cat("Tomato Data (not paired):\n")
attach(tomato.data)
a = pounds[fertilizer=="A"]
b = pounds[fertilizer=="B"]
print(round(test = permtest(b,a),3))
detach()

# Permutation test for Boy's Shoes Example
data(shoes.data)
cat("Shoes Data (paired):\n")
attach(shoes.data)
x = matB-matA
print(round(test = permtest(x),3))
detach()

#################################################################################
###############################################################

library(bootstrap)

# to do a leave-on-out jackknife estimate for the mean of the
# data ?jackknife gives an example
# Having a look at the jackknife function we see that it demands
# two parameters: x and theta. x is supposed to contain the data
# and theta the function that is applied to the data

x = rnorm(20)
theta = function(x){mean(x)}
results = jackknife(x,theta)

###############################################################
#### Jackknife
###############################################################

#Bootstrap
#Now I want to program the estimation for one coefficient of a linear regression model.
DF = as.data.frame(matrix(rnorm(250), ncol=5))    # my data
model.lm = formula(V1 ~ V2 + V3 + V4)             # my model
# Now I need to specify the theta function. Here x is not the
# data itself but is used as the row index vector to select
# a subset from the data frame (xdata). Also the coefficient
# to be returned is specified.
theta = function(x, xdata, coefficient)
{
coef(lm(model.lm, data=xdata[x,]))[coefficient]
}
# So now at each leave-on-out run the model is calculated with
# a subset defined by the vector x (here one is left out) and one
# coefficient is returned:
results = jackknife(1:50, theta, xdata=DF, coefficient="(Intercept)")
## To expand this code onto the estimation of all the regression coefficients is only a small step now. As the theta function is supposed to return a scalar and not a list of estimates for each coefficient, the following workaround is used: The sapply function calls the jackknife function four times prompting a different parameter estimate at each run. The prompted coeffient is passed on to the jackknife function by the three point (?) argument .
# The following function calculates all then coefficients
jackknife.apply = function(x, xdata, coefs)
{
sapply(coefs,
function(coefficient) jackknife(x, theta, xdata=xdata,
coefficient=coefficient),
simplify=F)
}
# now jackknife.apply() can be called
results = jackknife.apply(1:50, DF, c("(Intercept)", "V2", "V3", "V3"))

###############################################################
### Cross Validation
#################################################################################

#The idea: Models should be tested with data that were not used to fit the model. If you have enough data, it is best to hold back a random portion of the data to use for testing. Cross validation is a trick to get out-of-sample tests but still use all the data. The sleight of hand is to do a number of fits, each time leaving out a different portion of the data.
#Cross validation is perhaps most often used in the context of prediction. Everyone wants to predict the stock market. So let's do it.
#Below we predict tomorrow's IBM return with a quadratic function of today's IBM and S&P returns.
predictors = cbind(spxret, ibmret, spxret^2, ibmret^2,
spxret * ibmret)[-251,]
predicted = ibmret[-1]
predicted.lm = lm(predicted ~ predictors)
#The p-value for this regression is 0.048, so it is good enough to publish in a journal. However, you might want to hold off putting money on it until we have tested it a bit more.
#In cross validation we divide the data into some number of groups. Then for each group we fit the model with all of the data that are not in the group, and test that fit with the data that are in the group. Below we divide the data into 5 groups.
group = rep(1:5, length=250)
group = sample(group)
mse.group = numeric(5)
for(i in 1:5) {
group.lm = lm(predicted[group != i] ~
predictors[group != i, ])
mse.group[i] = mean((predicted[group == i] -
cbind(1, predictors[group == i,]) %*% coef(group.lm))^2)
}
#The first command above repeats the numbers 1 through 5 to get a vector of length 250. We do not want to use this as our grouping because we may be capturing systematic effects. In our case a group would be on one particular day of the week until a holiday interrupted the pattern. Hence the next line permutes the vector so we have random assignment, but still an equal number of observations in each group. The for loop estimates each of the five models and computes the out-of-sample mean squared error.
#The mean squared error of the predicted vector taking only its mean into account is 0.794. The in-sample mean squared error from the regression is 0.759. This is a modest improvement, but in the context of market prediction might be of use if the improvement were real. However, the cross validation mean squared error is 0.800 -- even higher than from the constant model. This is evidence of overfitting (which this regression model surely is doing).

#################################################################################
#### Boostrap
#################################################################################

There are a number of R packages that are either confined to or touch upon bootstrapping or its relatives. These include:
boot: This package incorporates quite a wide variety of bootstrapping tricks.
bootstrap: A package of relatively simple functions for bootstrapping and related techniques.
coin: A package for permutation tests (which are discussed below).
Design: This package includes bootcov for bootstrapping the covariance of estimated regression parameters and validate for testing the quality of fitted models by cross validation or bootstrapping.
MChtest: This package is for Monte Carlo hypothesis tests, that is, tests using some form of resampling. This includes code for sampling rules where the number of samples taken depend on how certain the result is.
meboot: Provides a method of bootstrapping a time series.
permtest: A package containing a function for permutation tests of microarray data.
resper: A package for doing restricted permutations.
scaleboot: This package produces approximately unbiased hypothesis tests via bootstrapping.
simpleboot: A package of a few functions that perform (or present) bootstraps in simple situations, such as one and two samples, and linear regression.
There are a large number of R packages that include bootstrapping. Examples include multtest that has the boot.resample function, and Matching which has a function for a bootstrapped Kolmogorov-Smirnov test (the equality of two probability distributions).

################################################################################################
library(boot)
### Single variable
#To use the boot package you will need to write a function to calculate the statistic of interest. The format is illustrated below for the sample mean, but any univariate function would be handled similarly. We'll call our function "boot.mean". When you have finished writing the script for a function you will need to cut and paste it into the command window so that R can access it (you'll need to do this just once in an R session). Here, "x" refers to the vector of data. "i" refers to a vector of indices, which must be included as an argument in the function and employed as shown.
boot.mean = function(x,i){boot.mean = mean(x[i])}
#The command "boot" will automatically carry out all the resampling and computations required. For this example, "x" is the vector of original data and "boot.mean" is the name of the function we created above to calculate the statistic of interest. R specifies the number of bootstrap replicate estimates desired.
z = boot(x, boot.mean, R=2000)
#The resulting object (which here named "z") is a boot object containing all the results. Use the following additional commands to pull out the results.
print(z)            # Bootstrap calculation of bias and SE
sd(z$t) # Another way to get the standard error hist(z$t)           # Histogram of boostrap replicate estimates
qqnorm(z$t) # Normal quantiles of replicate estimates boot.ci(z,type="bca") # 95% confidence interval using BCa boot.ci(z,type="perc") # Same using percentile method ### Two or more variables #Here's how you would use the boot command If the statistic of interest must be calculated from two (or more) variables (for example, a correlation, regression slope, or odds ratio). Assume that there are two variables of interest, "x" and "y". If not done already, put the two variables into a data frame (here called "mydata"), mydata = cbind.data.frame(x, y, stringsAsFactors=FALSE) # Then create a function to calculate the statistic of interest on the variables in "mydata". For example, to create a function that calculates the correlation coefficient between the two variables "x" and "y", use boot.cor = function(mydata,i) { x = mydata$x[i]
y = mydata$y[i] boot.cor = cor(x,y) } #Here, "i" refers to a vector of indices, which must be included as an argument in the function and employed as shown. #Finally, pass your data frame and function to the boot command, z = boot(mydata, boot.cor, R=2000) #See the previous section for a list of commands to pull results from the boot object (here named "z"). ################################################################################################ Basado en el difmediasTestBoot.pdf #Uso de library( boot ) Tres ejemplos sencillos ################################################################################################ #1) Estimaci?n del par?metro exp( media ) mediante xbar=mean( x ): ################################################################################################ set.seed( 89723 ) x = c(-1.034,-1.044,-1.445,-0.954,0.223,0.096,0.487,-0.485,0.491,-1.122) #Los datos x xbar = mean(x) n = length(x) expm.fun = function(x) { exp(mean(x)) } expm = expm.fun(x) cbind( n, xbar, expm ) #n xbar expm #[1,] 10 -0.4787 0.6195883 library( boot ) R = 999 #Es preciso definir una funci?n de los datos, x, y de los ?ndices, que ser?n tomados al azar al hacer bootstrap, que calcule el estad?stico a "bootstrapear": expm.i = function( x, indices ) { expmb = exp( mean( x[indices] ) ) c( expmb ) } #Bootstrap del estad?stico exp(media): expm.boot = boot( data=x, statistic=expm.i, R=R ) expm.boot #ORDINARY NONPARAMETRIC BOOTSTRAP #Call: #boot(data = x, statistic = expm.i, R = R) #Bootstrap Statistics : #original bias std. error #t1* 0.6195883 0.01076195 0.1456947 #Boostrap del intervalo de confianza. La funci?n boot.ci aplicada a un objeto de tipo boot permite calcular diferentes IC bootstrap: boot.ci( expm.boot, conf=0.95, type=c("perc", "bca") ) #BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS #Based on 999 bootstrap replicates #CALL : #boot.ci(boot.out = expm.boot, conf = 0.95, type = c("perc", "bca")) #Intervals : #Level Percentile BCa #95% ( 0.3928, 0.9627 ) ( 0.4067, 0.9822 ) #Calculations and Intervals on Original Scale #Las r?plicas bootstrap se encuentran en el elemento expm.boot$t del objeto expm.boot:
expmb = expm.boot$t ### valores bootstrap de expm sd( expmb ) [1] 0.1456947 ### se bootstrap de expm quantile( expmb, c( 0.025, 0.975) ) #2.5% 97.5% #0.3946921 0.9600645 #La funci?n plot aplicada al objeto de tipo boot, expm.boot, produce por defecto un histograma de expmb y un plot de normalidad: windows( width=7, height=4 ) par( pty="s" ) plot( expm.boot ) ### funci?n plot.boot: produce dos plots title( main="Resultados de expm.boot", font.main=4, line=3, col.main="blue" ) x11() hist( expmb, breaks=50, probability=T ) abline( v=c(expm,quantile(expmb,c(0.025,0.975))), lty=3 ) rug( expm, lwd=2, col="blue" ) lines( density( expmb ) ) ################################################################################################ #### 2) Regresi?n Bootstrap sobre la muestra de pares. IC para beta1 ################################################################################################ ### y respuesta; x explicativa; IC para beta1 x = c(109,88,96,96,109,116,114,96,85,100,113,117,107,104,101,81) y = c(116,77,95,79,113,122,109,94,91,88,115,119,100,115,95,90) lmyx = lm( y ~ x ) summary( lmyx ) #Call: #lm(formula = y ~ x) #Coefficients: #Estimate Std. Error t value Pr(>|t|) #(Intercept) -9.3016 20.0530 -0.464 0.65 #x 1.0826 0.1955 5.537 7.32e-05 *** #Residual standard error: 8.414 on 14 degrees of freedom #Multiple R-squared: 0.6865, Adjusted R-squared: 0.6641 #F-statistic: 30.66 on 1 and 14 DF, p-value: 7.319e-05 library( boot ) R = 999 datos = cbind(x, y) beta1.i = function( datos, i ) { coef( lm( datos[,2][i] ~ datos[,1][i] ) )[2] } beta1.b = boot( data=datos, statistic=beta1.i, R=R ) beta1.b # ORDINARY NONPARAMETRIC BOOTSTRAP # Call: # boot(data = datos, statistic = beta1.i, R = R) # Bootstrap Statistics : # original bias std. error # t1* 1.082613 0.02671124 0.2162240 boot.ci( beta1.b, type=c("perc", "bca") ) # BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS # Based on 999 bootstrap replicates # CALL : # boot.ci(boot.out = beta1.b, type = c("perc", "bca")) # Intervals : # Level Percentile BCa # 95% ( 0.795, 1.653 ) ( 0.774, 1.628 ) # Calculations and Intervals on Original Scale windows( width=7, height=4 ) par( pty="s" ) plot( beta1.b, cex.main=0.8 ) title( main="Regresi?n lineal: beta1.b", font.main=4, line=3, col.main="blue") x11() hist( beta1.b$t, probability=T, breaks=30, main="Histograma de beta1.b", cex.main=0.9 )
abline(v=c( coef(lmyx)[2], quantile(beta1.b$t,c(0.025,0.975))), lty=3) Histograma de beta1.bbeta1.b$tDensity0.60.81.01.21.41.61.80.00.51.01.52.0
rug( coef(lmyx)[2], lwd=2, col="blue" )
lines( density(beta1.b$t) ) ################################################################################################ #### 1) Bootstrap para la comparaci?n de medias (IC y tests bootstrap) ################################################################################################ #Dadas dos muestras, x e y, de sendas poblaciones independientes, con distribuciones F y G respectivamente, se quiere calcular un IC para la diferencia de medias. #Por ejemplo: ?Cu?l es la diferencia de medias del consumo diario de energ?a de individuos normales y obesos? #Consumo diario (MJ/d?a) en un grupo de 13 individuos de peso normal: 6.13, 7.05, 7.48, 7.48, 7.53, 7.58, 7.90, 8.08, 8.09, 8.11, 8.40, 10.15, 10.88. #Consumo en un grupo de 9 individuos obesos: 7.67, 8.19, 8.21, 8.79, 9.11, 9.21, 10.34, 11.53, 12.10. ### Aplicaci?n: IC percentil bootstrap para la diferencia de medias x = c(6.13, 7.05, 7.48, 7.48, 7.53, 7.58, 7.90, 8.08, 8.09, 8.11, 8.40, 10.15, 10.88) ### x: consumo diario de energia en individuos de peso normal y = c(7.67, 8.19, 8.21, 8.79, 9.11, 9.21, 10.34, 11.53, 12.10) ### y: consumo diario de energia en individuos obesos n = length(x) m = length(y) mx = mean(x) my = mean(y) se.x = sd(x)/sqrt(n) se.y = sd(y)/sqrt(m) ### Estimaci?n de la diferencia de medias ### dif = my-mx ### diferencia de medias observadas se1.dif = sqrt((var(y)/m)+(var(x)/n)) ### se varianzas distintas se2.dif = sqrt((((n-1)*var(x))+((m-1)*var(y)))/(n+m-2))*sqrt((1/m)+(1/n)) ### se2.dif: se para varianzas iguales round(cbind( n, mx, se.x, m, my, se.y, dif, se1.dif, se2.dif ), 3) #n mx se.x m my se.y dif se1.dif se2.dif #[1,] 13 8.066 0.343 9 9.461 0.514 1.395 0.618 0.593 ###Se puede obtener f?cilmente un IC basado en la t de student para la diferencia de medias: alfa = 0.05; confianza = 1-alfa; infdif = dif - qt(1-alfa/2, n+m-2) * se2.dif supdif = dif + qt(1-alfa/2, n+m-2) * se2.dif round(cbind( dif, se2.dif, infdif, supdif, confianza ), 3) #dif se2.dif infdif supdif confianza #[1,] 1.395 0.593 0.157 2.632 0.95 #Mediante bootstrap podemos obtener un IC para la diferencia de medias utilizando simplemente su estimador, la diferencia de medias muestrales. Usando R=1000 r?plicas bootstrap, realizadas por separado, para x, y para y, se simula la distribuci?n bootstrap de la diferencia de medias muestrales, de la que obtenemos el IC bootstrap percentil para la diferencia de medias. ### bootstrap para la diferencia de medias: IC percentil bootstrap ### R = 1000 set.seed( 48761 ) tm = matrix(sample(y,R*m,replace=T), R, m) ### r?plicas bootstrap separado cm = matrix(sample(x,R*n,replace=T), R, n) tb = apply(tm, 1, mean) cb = apply(cm, 1, mean) difb = tb-cb ### replicas de diferencia de medias bootstrap c( mean(difb), sd(difb)) #[1] 1.3928984 0.5880266 ppp = c(0.025,0.05,0.5,0.95,0.975) round(quantile(difb, ppp), 4) ### percentiles bootstrap de diferencia de medias #2.5% 5% 50% 95% 97.5% #0.2375 0.3918 1.4011 2.3413 2.5311 x11() hist( difb, breaks=25, probability=T, main="histograma de difb", cex.main=0.9, xlab="difb: diferencia de medias bootstrap" ) rug( dif, lwd=2, col="blue" ) abline( v=c(quantile(difb,0.025), quantile(difb,0.975), dif), lty=3 ) curve( dnorm(x, mean(difb), sd(difb)), add=T, lty=3 ) ################################################################################################ #### 2) Bootstrap de igualdad de distribuciones de dos muestras independientes ################################################################################################ #Problema .Hay diferencias en la distribucion del consumo diario de energia de individuos normales y obesos? H0: F=G. #El problema puede ser resuelto mediante un test de aleatorizacion, pero la metodologia bootstrap nos proporciona otros procedimientos para resolver el mismo problema. #Supongamos que para hacer el contraste de igualdad de distribuciones utilizamos el estadistico diferencia de medias muestrales D=meanY-meanX, estimador de la diferencia de medias mu1-mu2. de ambas poblaciones, aunque D sea mas apropiado para valorar solo la diferencia de medias. #Fijado un numero R de replicas, de la muestra conjunta z = (x, y) extraemos al azar y con reemplazamiento n+m valores en cada replica. #Con cada replica r, tenemos un vector en el que llamaremos a los n primeros elementos e a los m ultimos, calculando con ellos el valor de la replica r bootstrap de la diferencia de medias *rz*rx*ry***rxryrd.=. #Disponemos asi de las R simulaciones de la distribucion bootstrap de la diferencia de medias con las que podemos aproximar el valor-p bootstrap de nuestro test: p-valor bootstrap~#(dr*>d)/R, donde d es el valor observado de la diferencia de medias muestrales D. #El proceso es similar al llevado a cabo en el test de aleatorizacion. La diferencia estriba en que ahora se extraen muestras con reemplazamiento de la muestra conjunta z = (x, y). #Aplicaci?n del procedimiento: z = c(x, y) zm = matrix( sample(z, R*(n+m), replace=T), R, n+m ) ### remuestreo bootstrap conjunto xzm = apply( zm[ , 1:n], 1, mean ) yzm = apply( zm[ , (n+1):(n+m)], 1, mean ) difzb = yzm-xzm ### replicas de diferencia de medias bootstrap pv.difzb = mean( difzb >= dif ) ### p-valor bootstrap de la difer de medias cbi #nd( dif, pv.difzb) dif pv.difzb #[1,] 1.394957 0.026 x11() #par( mfrow=c(2,2), pty="s" ) hist( difzb, breaks=20, probability=T, main="histograma de difzb", xlab="difzb: diferencia de medias bootstrap", cex.main=0.9 ) rug( dif, lwd=2, col="blue" ) abline( v=dif, lty=3 ) curve( dnorm(x, mean(difzb), sd(difzb)), add=T, lty=3 ) #El mismo procedimiento puede ser aplicado a cualquier otro test para contrastar la igualdad de distribuciones. Veamos el uso del test t (dise?ado para la comparaci?n de medias, cuando las varianzas de ambas poblaciones, supuestamente normales, son iguales): ### R?plicas del test t bootstrap basado en la muestra conjunta ### sumx2zm = apply( zm[ , 1:n], 1, var )*(n-1) sumy2zm = apply( zm[ , (n+1):(n+m)], 1, var )*(m-1) se2b = sqrt( (sumx2zm+sumy2zm)/(n+m-2) )*sqrt((1/m)+(1/n)) t2 = difzb/se2b ### replicas del estadistico t bootstrap tobs2 = dif/se2.dif ### estadistico t observado tobs2 #[1] 2.351313 pv.tboot = mean( t2>=tobs2 ) ### p-valor bootstrap del t bootstrap (un lado) pv.t = 1-pt( tobs2,n+m-2 ) ### p-valor del test t para igualdad de medias round(cbind( dif, pv.difzb, tobs2, pv.t, pv.tboot ), 4) #dif pv.difzb tobs2 pv.t pv.tboot #[1,] 1.395 0.026 2.3513 0.0145 0.024 x11() hist( t2, breaks=25, probability=T, main="histograma de t2", cex.main=0.9, xlab="t2: estad?stico t bootstrap para diferencia de medias" ) abline( v=tobs2, lty=3 ) curve( dnorm(x, mean(t2), sd(t2)), add=T, lty=3 ) #Obs?rvese que no hay diferencia de procedimiento en la obtenci?n del p-valor del test basado en la diferencia de medias o en el estad?stico t de student, y el resultado es similar (ver pv.difzb y pv.tboot). En cambio, la diferencia con el p-valor del test t de student exacto es notable (ver pv.t). Hay que insistir en que contrastar F=G no es equivalente a contrastar la igualdad de sus medias. ################################################################################################ #### 3) Bootstrap de igualdad de medias de muestras independientes ################################################################################################ #Si unicamente tuvieramos interes en contrastar la hipotesis 21:0?g?g=H, mas debil que la anterior H0 : F = G, entonces lo mas adecuado es utilizar el test t-Student si se cumplen las condiciones de independencia, normalidad e igualdad de varianzas. #Si las varianzas fueran diferentes podria plantearse el uso de otro estadistico t con una estimacion de las varianzas por separado, 221211smsnXYt+.=, donde en el denominador tenemos las varianzas muestrales de x e y. #Dicho estadistico no tiene bajo H0 una distribucion t-Student y solo se dispone de soluciones aproximadas para obtener el valor-p. Es conocido en la estadistica clasica como ??problema de Behrens-Fisher??. #El uso del bootstrap permite soslayar este problema, cuando tratamos de contrastar 21:0?g?g=H. Bajo esta hipotesis, es admisible que las distribuciones de ambas poblaciones, F y G, puedan ser diferentes. #Los datos x nos dan la distribucion empirica para estimar F y los datos y la G para estimar G. Para calcular el valor-p de un test mediante bootstrap necesitaremos que ambas distribuciones empiricas sean del tipo indicado por F..21:0?g?g=H, esto es tengan iguales medias. Una forma razonable de conseguirlo sin alterar otras caracterirticas de ambas distribuciones consiste en relocalizar los datos x e y observados, de forma que mediante dicha traslacion ambas muestras presenten la misma media, por ejemplo la media de la muestra conjunta. #Si z = (x, y) es la muestra conjunta, con media muestral )(1ymxnmnz++=, consideramos zxxxloc+.= e zyyyloc+.= y adoptamos ahora xloc e yloc como si fueran nuestros datos de partida, teniendo ambas muestras como media muestral el valor z. #Con xloc e yloc extraemos por separado muestras bootstrap xloc* e yloc*, calculando con cada una de ellas sus medias y varianzas bootstrap para obtener las replicas bootstrap del estadistico t, t1*, ..., tR*. #El valor-p bootstrap es Robstrt)*(#., donde es el valor observado del estadistico t. obst #A continuacion aplicamos el procedimiento, relocalizando las muestras de modo que ambas tengan la misma media comun global. #Aplicaci?n xloc = x-mean(x)+mean(z) ### controles relocalizados en la media global yloc = y-mean(y)+mean(z) ### tratados relocalizados ### bootstrap separado para test de igualdad de medias xlocm = matrix(sample(xloc,R*n,replace=T),R,n) ylocm = matrix(sample(yloc,R*m,replace=T),R,m) xlocb = apply(xlocm,1,mean) ylocb = apply(ylocm,1,mean) dyxlocb = ylocb-xlocb ### diferencias de medias bootstrap tras relocalizacion mean( dyxlocb >= dif ) ### p-valor bootstrap de test basado en diferencia de medias #[1] 0.01 x11() hist(dyxlocb, breaks=20, probability=T) rug( dif, lwd=2, col="blue" ) abline( v=dif, lty=3 ) curve( dnorm(x, mean(dyxlocb), sd(dyxlocb)), add=T, lty=3 ) ### usando el estad?stico t vxlocb = apply(xlocm, 1, var) vylocb = apply(ylocm, 1, var) ## replicas bootstrap del estadistico t con varianzas distintas tyxlocb = dyxlocb/sqrt((vxlocb/n)+(vylocb/m)) mean(tyxlocb>=tobs2) ### p-valor bootstrap de tyxlocb #[1] 0.01 x11() hist( tyxlocb, breaks=20, probability=T ) rug( tobs2, lwd=2, col="blue" ) abline( v=tobs2, lty=3 ) curve( dnorm(x, mean(tyxlocb), sd(tyxlocb)), add=T, lty=3 ) ################################################################################################ #### 4) CORRELACI?N BOOTSTRAP ################################################################################################ #Se quiere conocer la correlaci?n entre las puntuaciones de expresi?n verbal y capacidad matem?tica de los estudiantes de un curso. Para ello se ha pasado un test a 16 estudiantes y se han obtenido los resultados de x:expresi?n verbal e y:capacidad matem?tica. (Pawitan) ### Datos IQ: x expresi?n verbal, y capacidad matem?tica ### x = c(109,88,96,96,109,116,114,96,85,100,113,117,107,104,101,81) y = c(116,77,95,79,113,122,109,94,91,88,115,119,100,115,95,90) n = length(x) xm = mean(x); ym = mean(y); s2x = var(x); sdx = sqrt(s2x) s2y = var(y); sdy = sqrt(s2y) cxy = cor(x, y) plot(x, y, xlab="verbal", ylab="matematica", main="datos IQ") abline(lsfit(x,y)$coef, lty=3)

round( cbind(n, xm, ym, s2x, sdx, s2y, sdy, cxy), 4)
#n xm ym s2x sdx s2y sdy cxy
#[1,] 16 102 101.125 123.4667 11.1116 210.7833 14.5184 0.8286

#Bajo normalidad conjunta de (X, Y), un estimador clasico del error estandar de la correlacion muestral rho. viene dado por (1-rho^2)/sqrt(n-3).

#Si se considera la transformacion de Fisher tau=0.5log((1+rho)/(1-rho)) entonces el estimador ?n. es aproximadamente normal de media tau y desviacion estandar 1/sqrt(n-3).
#Esto permite obtener IC's para tau y aplicando la transformacion inversa de Fisher a los extremos del IC, obtendremos un I.C. para la correlacion.

### se de la correlacion bajo normalidad conjunta de (X,Y)
se.cxy = (1-cxy^2)/sqrt(n-3)

### tau: transformacion de Fisher de la correlacion
tau = 0.5*log( (1+cxy)/(1-cxy) )

alfa = 0.05; confianza = 1-alfa;
### IC para tau
tau.inf = tau-qnorm(1-alfa/2)*(1/sqrt(n-3))
tau.sup = tau+qnorm(1-alfa/2)*(1/sqrt(n-3))
### IC para la correlacion
cxy.inf = (exp(2*tau.inf)-1)/(exp(2*tau.inf)+1)
cxy.sup = (exp(2*tau.sup)-1)/(exp(2*tau.sup)+1)

round(cbind( n, cxy, se.cxy, cxy.inf, cxy.sup, confianza ), 4)
#n cxy se.cxy cxy.inf cxy.sup confianza
#[1,] 16 0.8286 0.0869 0.5649 0.9387 0.95

library( boot )
R = 999
datos = cbind(x,y)
correla = function( datos, indices )
{
x = datos[,1][indices]
y = datos[,2][indices]
cor( x[indices], y[indices] )
}

corxy.b = boot( data=datos, statistic=correla, R=R )
corxy.b
#ORDINARY NONPARAMETRIC BOOTSTRAP
#Call:
#boot(data = datos, statistic = correla, R = R)
#Bootstrap Statistics :
#original bias std. error
#t1* 0.8285718 0.01169248 0.08414495

boot.ci( corxy.b, conf=c(0.9, 0.95), type="all" )
#BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#Based on 999 bootstrap replicates
#CALL :
#boot.ci(boot.out = corxy.b, conf = c(0.9, 0.95), type = "all")
#Intervals :
#Level Normal Basic
#90% ( 0.6785, 0.9553 ) ( 0.7048, 0.9575 )
#95% ( 0.6520, 0.9818 ) ( 0.6929, 1.0093 )
#Level Percentile BCa
#90% ( 0.6997, 0.9523 ) ( 0.5930, 0.9213 )
#95% ( 0.6478, 0.9642 ) ( 0.5170, 0.9407 )
#Calculations and Intervals on Original Scale
#Some BCa intervals may be unstable
#Warning message:
#In boot.ci(corxy.b, conf = c(0.9, 0.95), type = "all") : bootstrap variances needed for studentized intervals

windows( width=7, height=4 )
par( pty="s" )
plot( corxy.b, cex.main=0.8 )
title( main="Correlaci?n: corxy.b", font.main=4, line=3, col.main="blue" )

##### Bootstarap No Parametrico
B = 1000
set.seed( 98725 )
corxy.boot = NULL
for(i in 1:B)
{
star = sample(1:n, n, replace=T)
xstar = x[star]
ystar = y[star]
corxy.boot[i] = cor(xstar, ystar)
}
qcorxy.boot = quantile(corxy.boot, c(0.025, 0.05, 0.5, 0.95, 0.975))
round(qcorxy.boot, 4)
#2.5% 5% 50% 95% 97.5%
#0.7136 0.7341 0.8415 0.9258 0.9365
#Los percentiles 2.5% y 97.5% determinan los extremos de un IC bootstrap percentil del 95%.
bias.corxy = cxy - mean(corxy.boot) ### sesgo bootstrap
se.corxy.boot = sd(corxy.boot) ### se bootstrap
round(cbind(cxy, se.corxy.boot, bias.corxy, mean(corxy.boot)), 4)
#cxy se.corxy.boot bias.corxy
#[1,] 0.8286 0.0574 -0.0091 0.8377

#El error estandar bootstrap de rho. se ha calculado mediante...
#El estimador de la correlacion es sesgado, pero en este caso el sesgo estimado es pequeno, de modo que la correccion del sesgo apenas cambiara el valor estimado de la correlacion. El estimador bootstrap del sesgo es (1/B)*sum(rhp^p*-rho), y por lo tanto, un estimador de la correlacion con correccion del sesgo viene dado por: 2rho-(1/B)*sum(rho^b).
#Si la distribucion de la correlacion muestral fuera aproximadamente normal, un IC del tipo Wald, seria muy sencillo de calcular, utilizando la estimacion bootstrap de su error estandar:
### IC tipo Wald del 95% ###
wic.inf = cxy - qnorm(0.975) * se.corxy.boot
wic.sup = cxy + qnorm(0.975) * se.corxy.boot
round(cbind( wic.inf, cxy, wic.sup), 4)
#wic.inf cxy wic.sup
#[1,] 0.7161 0.8286 0.941
x11()
hist( corxy.boot, probability=T, breaks=25, main="" )
rug( qcorxy.boot )
abline(v=c(qcorxy.boot[c(1, 5)], cxy), lty=3) ### IC percentil bootstrap 95%
rug( cxy, lwd=2, col="blue" )

#La transformaci?n tau de Fisher, bootstrap:
tau.boot = 0.5*log( (1+corxy.boot)/(1-corxy.boot) )
x11()
hist( tau.boot, probability=T, breaks=25, main="" )
curve( dnorm(x, tau, 1/sqrt(n-3)), lty=3, add=T )
qb.025 = quantile( corxy.boot, 0.025 )
qb.975 = quantile( corxy.boot, 0.975 )
abline( v=c(tau,0.5*log((1+qb.025)/(1-qb.025)),0.5*log((1+qb.975)/(1-qb.975))), lty=3 )

#calculating the standard error of the median
#creating the data set by taking 100 observations
#from a normal distribution with mean 5 and stdev 3
#we have rounded each observation to nearest integer
data = round(rnorm(100, 5, 3))
data[1:10]
[1] 6 3 3 4 3 8 2 2 3 2
#obtaining 20 bootstrap samples
#display the first of the bootstrap samples
resamples = lapply(1:20, function(i)
sample(data, replace = T))
resamples[1]
[[1]]:
[1]  5  1  7  6  5  2  2  6  9  5  4  6  6  3  5  4 10  7  8  1  8  0  5  2
[25]  8  3  0  9  3  2  3 10  5  8  5  4  0  4  7  3  5  6  3  6  3  2  9  7
[49]  2  4  9  6  6  0  7  5  9  3  0  6  8  5  2  3  3  3  4  3  2  9  3  3
[73]  2  3  8  2  8  3  9  6  5  2  4  3  3  7  1  3  5  9  4  3  4  2  9  0
[97]  3  6  9  7
#calculating the median for each bootstrap sample
r.median = sapply(resamples, median)
r.median
[1] 4.0 4.5 4.0 5.0 4.0 5.0 5.0 5.0 5.0 4.0 5.0 5.0 5.0 5.0 5.0 4.0 5.0 5.0
[19] 6.0 5.0
#calculating the standard deviation of the distribution of medians
sqrt(var(r.median))
[1] 0.5250313
#displaying the histogram of the distribution of the medians
hist(r.median)

#We can put all these steps into a single function where all we would need to specify is which data set to use and how many times we want to resample in order to obtain the adjusted standard error of the median. For more information on how to construct functions please consult the R library pages on introduction to functions and advanced functions.
#function which will bootstrap the standard error of the median
b.median = function(data, num) {
resamples = lapply(1:num, function(i) sample(data, replace=T))
r.median = sapply(resamples, median)
std.err = sqrt(var(r.median))
list(std.err=std.err, resamples=resamples, medians=r.median)
}
#generating the data to be used (same as in the above example)
data1 = round(rnorm(100, 5, 3))
#saving the results of the function b.median in the object b1
b1 = b.median(data1, 30)
#displaying the first of the 30 bootstrap samples
b1$resamples[1] [[1]]: [1] 6 6 6 0 9 3 5 6 2 7 5 3 3 2 5 3 0 4 3 7 3 1 6 3 [25] 6 3 5 9 6 4 8 4 3 1 3 2 4 5 2 3 1 7 6 3 7 10 7 2 [49] 3 3 5 10 9 6 2 3 4 5 1 3 5 0 8 1 4 2 7 8 2 2 10 6 [73] 4 8 6 3 5 2 10 5 0 7 6 5 4 9 6 0 6 6 3 4 8 10 7 6 [97] 3 3 2 3 #displaying the standard error b1$std.err
[1] 0.5155477
#displaying the histogram of the distribution of medians
hist(b1$medians) #we can input the data directly into the function and display #the standard error in one line of code b.median(rnorm(100, 5, 2), 50)$std.err
[1] 0.5104178
#displaying the histogram of the distribution of medians
hist(b.median(rnorm(100, 5, 2), 50)$medians) ### Basado en http://www.statoo.com/en/publications/bootstrap_scgn_v131.pdf # los datos son de un experimento en el cual se utilizan aleatoriamente dos tratamientos de l?ser en los ojos de los pacientes. La respuesta es la adecuaci?n visual, medida por el n?mero de letras identificadas correctamente en un test est?ndard de ojos. # algunos pacientes tienen solo un ojo sano, y reciben un tratamiento localizado aleatoriamente. Existen 20 apcientes con datos pareados y 20 pacientes con con solo 1 observaci?n disponible, por lo que tenemos una mezcla de comparaciones pareadas y dos muestras. blue=c(4, 69, 87, 35, 39, 79, 31, 65, 95, 68, 62, 70, 80, 84, 79, 66, 75, 59, 77, 36, 86, 39, 85, 74, 72, 69, 85, 85, 72) red=c(62, 80, 82, 83, 0, 81, 28, 69, 48, 90, 63, 77, 0, 55, 83, 85, 54, 72, 58, 68, 88, 83, 78, 30, 58, 45, 78, 64, 87, 65) acui=data.frame(str=c(rep(0,20), rep(1,10)), red, blue) # solo consideraremos los datos pareados y construiremos el cl?sico intervalo de confianza al 95% del estad?stico t-Student para la media de las diferencias. acu.pd=acui[acui$str==0,]
dif=acu.pd$blue-acu.pd$red
n=nrow(acu.pd)
tmp=qt(.025, n-1)*sd(dif)/sqrt(n)
c(mean(dif)+tmp, mean(dif)-tmp) # [1] -.9270336 15.770335
# pero un q-q plot de las diferencias parece ma?s una Cauchy que una normal, por lo cual el modelo usual parece ser inviable. Para realizar un boostrap no-param?trico en este caso, necesitamos definir la funci?n boostrap:
acu.pd.fun=function(data, i)
{
d=data[i,]
dif=d$blue-d$red
c(mean(dif), var(dif)/nrow(d))
}
# obtenemos R=999 r?plicas bootstrap
acui.pd.b=boot(acu.pd, acu.pd.fun, R=999)
# los IC no-param?tricos resultantes pueden calcularse como antes, o utilizando directamente:
boot.ci(acu.pd.b, type=c("norm", "basic", "stud"))
#Normal           Basic             Studentized
#(-8.20, 14.95)    (-8.10, 15.05)    (-8.66, 15.77)

# una alternatica es considerar solo los datos de dos muestras y comparar las medias de las dos poblaciones issuing los pacientes que solo tienen una observacion disponible
acu.ts=acui[acui$str==1,] # calculamos el intervalo de confianza al 05% para la diferencia de medias acu.ts=acui[acui$str==1,]
dif=mean(acu.ts$blue)-mean(acu.ts$red)
tmp=qnorm(0.025)*sqrt(var(acu.ts$blue)/nrow(acu.ts)+var(acu.ts$red)/nrow(acu.ts))
c(dif+tmp, dif-tmp) #[1] -13.76901   19.16901
# para construir IC boostrap generamos R=999 r?plicas del estimador y varianza estimada
y=c(acui$blue[21:30], acui$red[21:30])
acu=data.frame(col=rep(c(1,2), c(10,10)), y)
acu.ts.f=function(data, i)
{
d=data[i, ]
m=mean(d$y[1:10])-mean(d$y[11:20])
v=var(d$y[1:10])/10+var(d$y[11:20])/10
c(m,v)
}
acu.ts.boot=boot(acu, acu.ts.f, R=999, strata=acu$col) # strata=acu$col asegura la simulaci?n estratificada

# Here's a worked example for comparing group averages with bootstrap confidence intervals and allowing for different subsample sizes by calling the strata argument within the bootstrap function.
# The data is set up analogous to an before-after impact experiment conducted on plots across four different successional stages. Similarity was calculated for plot's composition before and after an impact. One hypothesis was that certain stages would show higher / smaller average similarities, that is, a higher / lower impact on composition. As plots within stages were situated within different subsites and the nr. of replicates was unbalanced, this had to be allowed for by use of the "strata" argument in the boot.ci call:
library(boot)
library(Hmisc)

### a factor 'stage' and a variable 'MH-Index'
stage=sort(as.factor(rep(LETTERS[1:4], 30)))
site=gl(12, 10)
MH.Index=runif(4*30, 0, 1)

### simulate effect:
n = 30
eff = c(rep(0.25, n), rep(0.4, n), rep(0.5, n), rep( 0.65,30))

MH.Index = MH.Index*eff

### dataframe:
print(sim=data.frame(stage, site, MH.Index))

### function needed for the bootstrap:
mean.fun = function(x, index){mean(x[index])}

### set up dataframe to collect results for plotting:
pldat=data.frame(mean=rep(NA,4),low=rep(NA,4),upp=rep(NA,4))
rownames(pldat)=c("A","B","C","D")
pldat

for(i in c("A","B","C","D")){
### the boot sample and the ci's:
boot=boot(sim$MH.Index[sim$stage==i], mean.fun, R=1000, strata=sim$site[sim$stage==i])
ci=boot.ci(boot,type = "norm")
### get ci's (method: normal)
ci[2]->pldat[i,"mean"] # mean
ci$normal[1,2]->pldat[i,"low"] # lower ci ci$normal[1,3]->pldat[i,"upp"]} # upper ci
pldat

### plot:
errbar(c(1:4),pldat$mean,pldat$upp,pldat\$low,ylab="MH-Similarity Index",xlab="Stage",
pch=15,xaxt="n",xlim=c(0.5,4.5))
axis(1,c(1:4),labels=row.names(pldat))
legend("topleft","errorbars =\n95%-normal\nbootstrap\nconf. int",bty="n")


### EJEMPLOS: en el siguiente link podrán ver aplicaciones del boostrap en índices de biodiversidad.

http://www.wcsmalaysia.org/analysis/Biod_indices_bstrap.htm