Prueba de hipótesis para el cociente de varianzas

Suponga que se tienen dos muestras aleatorias que provienen de poblaciones normales así:

$n_1$ observaciones $X_{11}, X_{12}, \ldots, X_{1,n1}$ de una población I con varianza $\sigma^2_1$,
$n_2$ observaciones $X_{21}, X_{22}, \ldots, X_{2,n2}$ de una población II con varianza $\sigma^2_2$,
ambas muestras son independientes entre sí.

En este problema se desea estudiar la hipótesis nula $H_0: \sigma_1^2 / \sigma_2^2 = 1$ de igualdad de varianzas y se sospecha que el cociente de varianzas podría estar en alguna de las siguientes situaciones:

$H_a: \sigma_1^2 / \sigma_2^2 < 1$
$H_a: \sigma_1^2 / \sigma_2^2 \neq 1$
$H_a: \sigma_1^2 / \sigma_2^2 > 1$

El estadístico para realizar la prueba es: \[ F_0=\frac{S_1^2}{S_1^2}, \] donde $S_1^2$ y $S_2^2$ son las varianzas de las muestras I y II respectivamente. El estadístico $F_0$, bajo la suposición de que $H_0$ es verdadera, tiene distribución $f$ con $n_1-1$ grados de libertad en el numerador y $n_2-1$ grados de libertad en el denominador.

Si el valor calculado para el estadístico dado en la ecuación anterior se denota por $f_0$, entonces el valor-$P$ de la prueba se calcula de acuerdo a la hipótesis alterna $H_a$ así:

Si $H_a: \sigma_1^2 / \sigma_2^2 < \sigma^2_0$ entonces valor-$P$=$P(f_{n_1-1,n_2-1} < f_0)$.
Si $H_a: \sigma_1^2 / \sigma_2^2 \neq \sigma^2_0$, valor-$P$=$2 \times \min \left\{ P(f_{n_1-1,n_2-1} < f_0), P(f_{n_1-1,n_2-1} > f_0) \right\}$.
Si $H_a: \sigma_1^2 / \sigma_2^2 > \sigma^2_0$ entonces valor-$P$=$P(f_{n_1-1,n_2-1} > f_0)$.

En esta prueba, al no rechazar la hipótesis nula $H_0$, se concluye que $\sigma_1^2 / \sigma_2^2 = 1$ lo que implica en términos prácticos que $\sigma_1^2 = \sigma_2^2$, es decir que las varianzas poblacionales se pueden considerar iguales.

library(shiny) shinyServer(function(input,output,session){ observe({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('unequal_var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) updateSelectInput(session, "variable1", choices=names(dt[!sapply(dt, is.factor)])) # Asegurar cuanti updateSelectInput(session, "variable2", choices=names(dt[sapply(dt, is.factor)])) # Asegurar cuali }) output$summary <- renderTable({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('unequal_var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) dt <- na.omit(dt) # Para eliminar obs con NA dt }) output$statistic <- renderTable({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('unequal_var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) dt <- na.omit(dt) # Para eliminar obs con NA x <- dt[, input$variable1] # Variable de interes group <- dt[, input$variable2] # Variable de clasificacion if (nlevels(group) != 2) group <- dt[, sapply(dt, is.factor)] xx <- split(x, group) # Lista con variable interes resumen <- function(x) c(mean(x), var(x), length(x)) res <- sapply(xx, resumen) rownames(res) <- c('Media', 'Varianza', 'Número obs') t(res) }, rownames = TRUE, align='c', bordered = TRUE) # Para obtener tabla con rownames output$appPlot <- renderPlot({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('unequal_var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) # Aqui inicia la figura par(mfrow=c(1, 2), bg='gray98') dt <- na.omit(dt) # Para eliminar obs con NA x <- dt[, input$variable1] group <- dt[, input$variable2] if (nlevels(group) != 2) group <- dt[, sapply(dt, is.factor)] # Para dibujar las densidades xx <- split(x, group) den <- lapply(xx, density) plot(den[[1]], lwd=4, col='deepskyblue3', main='Densidad', las=1, xlab=as.character(input$variable1), ylab='Densidad', xlim=range(range(den[[1]]$x), range(den[[2]]$x)), ylim=c(0, max(c(den[[1]]$y, den[[2]]$y)))) lines(den[[2]], lwd=4, col='firebrick3') # Leyenda para distiguir las densidades legend('topright', bty='n', lwd=4, col=c('deepskyblue3', 'firebrick3'), legend=unique(group)) # Para dibujar los qqplot qq1 <- qqnorm(xx[[1]], plot.it=FALSE) qq2 <- qqnorm(xx[[2]], plot.it=FALSE) plot(qq1, las=1, main='QQplot', pch=19, col='deepskyblue3', xlim=range(c(qq1$x, qq2$x)), ylim=range(c(qq1$y, qq2$y)), xlab='Cuantiles teóricos N(0, 1)', ylab=as.character(input$variable1)) points(qq2, pch=19, col='firebrick3') # Para construir los qqplot qqline(xx[[1]], col='deepskyblue3') qqline(xx[[2]], col='firebrick3') # Para incluir el valor P de Shapiro shapi <- lapply(xx, shapiro.test) leyenda <- c(paste('Valor P=', round(shapi[[1]]$p.value, 2)), paste('Valor P=', round(shapi[[2]]$p.value, 2))) legend('topleft', bty='n', text.col=c('deepskyblue3', 'firebrick3'), legend=leyenda) }) output$resul1 <- renderText({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('unequal_var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) dt <- na.omit(dt) # Para eliminar obs con NA x <- dt[, input$variable1] group <- dt[, input$variable2] if (nlevels(group) != 2) group <- dt[, sapply(dt, is.factor)] xx <- split(x, group) ph <- var.test(x=xx[[1]], y=xx[[2]], alternative=input$h0, ratio=1, conf.level=input$alfa) paste0('El estadístico de prueba es f0=', round(ph$statistic, 4), ' con un valor-P de ', round(ph$p.value, 2), '.') }) output$resul2 <- renderText({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('unequal_var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) dt <- na.omit(dt) # Para eliminar obs con NA x <- dt[, input$variable1] group <- dt[, input$variable2] if (nlevels(group) != 2) group <- dt[, sapply(dt, is.factor)] xx <- split(x, group) ph <- var.test(x=xx[[1]], y=xx[[2]], alternative=input$h0, ratio=1, conf.level=input$alfa) intervalo <- paste("(", round(ph$conf.int[1], digits=4), ", ", round(ph$conf.int[2], digits=4), ").", sep='') paste0('El intervalo de confianza del ', 100*input$alfa, '% para el cociente de varianzas poblacionales es ', intervalo) }) })

library(shiny) library(markdown) shinyUI(pageWithSidebar( headerPanel(title=HTML("Prueba de hipótesis para el cociente de varianzas σ21 / σ22"), windowTitle="PH coc varianzas"), sidebarPanel( h5('Esta aplicación sirve para realizar prueba de hipotesis para el cociente de varianzas de variables cuantitativas con distribución normal.'), h6('La aplicación usa una base de datos de ejemplo pero el usuario puede cargar su propia base de datos.'), fileInput(inputId='file1', label='Use el botón siguiente para cargar su base de datos.', accept = c( 'text/csv', 'text/comma-separated-values', 'text/tab-separated-values', 'text/plain', '.csv', '.tsv' )), checkboxInput(inputId='header', label='¿Tiene encabezado la base de datos?', value=TRUE), selectInput(inputId="sep", label = "¿Cuál es la separación de los datos?", choices = list(Tab='\t', Comma=',', Semicolon=';', 'Space'=' '), selected = ';'), selectInput(inputId="variable1", label=p("Elija la variable", span("cuantitativa", style = "color:red"), "para realizar la prueba de hipótesis."), choices=""), selectInput(inputId="variable2", label=p("Elija la variable", span("cualitativa", style = "color:blue"), "de agrupacion, DEBE tener 2 niveles y ser un factor."), choices=""), selectInput(inputId="h0", label=HTML("La hipótesis nula de la prueba es Ho: σ21 / σ22 = 1, elija el tipo de hipotesis alterna < , ≠ o >"), choices=list("Menor" = "less", "Diferente" = "two.sided", "Mayor" = "greater"), selected = "two.sided"), sliderInput(inputId='alfa', label=HTML("Opcional: elija un nivel de confianza para construir el intervalo de confianza para el cociente σ21 / σ22"), min=0.90, max=0.99, value=0.95, step=0.01), img(src="logo.png", height = 60, width = 120), img(src="udea.png", height = 25, width = 70), img(src="cua.png", height = 40, width = 110), br(), br(), tags$a(href="https://srunal.github.io", "https://srunal.github.io") ), mainPanel( tabsetPanel(type = "pills", tabPanel(title="Resultados", h5('A continuación el histograma, densidad, QQplot y valor-P para la prueba de normalidad Shapiro-Wilk de cada una de las dos muestras.'), plotOutput("appPlot", width='500px', height='300px'), h4("- Tabla de resumen con estadísticos muestrales:"), tableOutput('statistic'), h4("- Resultados de la prueba de hipótesis:"), textOutput("resul1"), h4(HTML("- Intervalo de confianza para el cociente σ21 / σ22:")), textOutput("resul2")), tabPanel("Datos", "A continuación los datos que está usando la aplicación.", uiOutput('summary')), tabPanel("Teoría", includeHTML("include.html")) ) ) ))