Prueba de hipótesis para la varianza

Sea $X_1, X_2, \ldots, X_n$ una muestra aleatoria proveniente de una población normal con varianza $\sigma^2$ y supóngase que se desea probar la hipótesis de que $\sigma^2$ es igual a un valor de referencia $\sigma^2_0$, es decir, la hipótesis nula a estudiar es $H_0: \sigma^2 = \sigma_0^2$. En este problema las sospechas sobre la varianza $\sigma^2$ se resumen por medio de la hipótesis alterna ($H_a$) en una de tres situaciones como se muestra a continuación:

$H_a: \sigma^2 < \sigma_0^2$,
$H_a: \sigma^2 \neq \sigma_0^2$,
$H_a: \sigma^2 > \sigma_0^2$,

El estadístico para realizar la prueba es:

\[ X_0^2=\frac{(n-1) S^2}{\sigma_0^2}, \]

donde $S$ corresponde a la desviación estándar muestral. Bajo la suposición de que $H_0$ es verdadera, el estadístico $X_0^2$ tiene distribución $\chi^2$ con $n-1$ grados de libertad.

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

Si $H_a: \sigma^2 < \sigma_0^2$ entonces valor-$P$=$P(\chi^2_{n-1} < \chi_0^2)$.
Si $H_a: \sigma^2 \neq \sigma_0^2$ entonces valor-$P$=$2 \times \min \left\{ P(\chi^2_{n-1} < \chi_0^2),~ P(\chi^2_{n-1} > \chi_0^2) \right\}$.
Si $H_a: \sigma^2 > \sigma_0^2$ entonces valor-$P$=$P(\chi^2_{n-1} > \chi_0^2)$.

library(shiny) source('var.test.R') shinyServer(function(input,output,session){ observe({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) updateSelectInput(session, "variable", choices=names(dt)) }) output$statistic <- renderTable({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) y <- na.omit(dt[, input$variable]) # Para sacar los NA de la variable res <- data.frame(Min=min(y), Var=var(y), Max=max(y), n=length(y)) colnames(res) <- c('Mínimo', 'Varianza', 'Máximo', 'Número obs') res }, align='c', bordered=TRUE) output$appPlot <- renderPlot({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) par(mfrow=c(1, 2), bg='gray98') y <- na.omit(dt[, input$variable]) # Para sacar los NA de la variable hist(y, col='deepskyblue3', freq=F, las=1, xlab=as.character(input$variable), main='Histograma y densidad', ylab='Densidad') lines(density(y), lwd=4, col='firebrick3') qqnorm(y, las=1, main='QQplot', pch=19, col='deepskyblue3', ylab=as.character(input$variable)) qqline(y) ks <- ks.test(x=y, y=pnorm) legend('topleft', bty='n', col='red', text.col='deepskyblue3', legend=paste('Valor P=', round(ks$p.value,2))) }) output$inputData <- renderTable({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) dt }) output$resul1 <- renderText({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) y <- na.omit(dt[, input$variable]) # Para sacar los NA de la variable ph <- var.test(x=y, alternative=input$h0, null.value=input$sigma20, conf.level=input$alfa) paste0('El estadístico de prueba es ', round(ph$statistic, 2), ' con un valor-P de ', round(ph$p.value, 4), '.') }) output$resul2 <- renderText({ inFile <- input$file1 if(is.null(inFile)) dt <- read.table('var_data.txt', header=T, sep='\t') else dt <- read.csv(inFile$datapath, header=input$header, sep=input$sep) y <- na.omit(dt[, input$variable]) # Para sacar los NA de la variable ph <- var.test(x=y, alternative=input$h0, null.value=input$sigma20, 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 la varianza poblacional es ', intervalo) }) output$miteoria <- renderUI({ HTML(markdown::markdownToHTML(knit(input='include.md', quiet=TRUE))) }) })

library(shiny) library(markdown) shinyUI(pageWithSidebar( headerPanel(title=HTML("Prueba de hipótesis para la varianza σ2"), windowTitle="PH varianza"), sidebarPanel( h5('Esta aplicación realiza la prueba de hipótesis para la varianza de una variable cuantitativa 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 siguiente botón 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="variable", label="Elija la variable cuantitativa para realizar la prueba de hipótesis.", choices=""), numericInput(inputId='sigma20', label=HTML("Ingrese el valor de referencia σ20 para probar H0: σ2 = σ 20"), value=20, min=0, step=0.1), selectInput(inputId="h0", label=HTML("Elija la hipótesis alternativa < , ≠ 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 la varianza σ2"), 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(), tags$a(href="https://srunal.github.io", "https://srunal.github.io") ), mainPanel( tabsetPanel(type = "pills", tabPanel(title="Resultados", h5('A continuación el histograma, la densidad, el QQplot y valor-P de la prueba de normalidad Kolmogorov-Smirnov para la muestra.'), plotOutput("appPlot", width='500px', height='300px'), h4("- Tabla de resumen con los estadísticos muestrales:"), tableOutput('statistic'), h4("- Resultados de la prueba de hipótesis:"), textOutput("resul1"), h4(HTML("- Intervalo de confianza para la varianza σ2:")), textOutput("resul2") ), tabPanel("Datos", "A continuación los datos que está usando la aplicación.", uiOutput('inputData')), tabPanel("Teoría", includeHTML("include.html")) ) ) ))

var.test <- function(x, alternative='two.sided', null.value=1, conf.level=0.95) { alpha <- 1 - conf.level n <- length(x) if (alternative == 'two.sided') { alt <- 'not equal to' quantiles <- c(qchisq(p=alpha/2, df=n-1, lower.tail=F), qchisq(p=1-alpha/2, df=n-1, lower.tail=F)) conf.int <- (n-1) * var(x) / quantiles statistic <- (n-1) * var(x) / null.value p.value <- 2 * min(c(pchisq(statistic, n-1, lower.tail=F), pchisq(statistic, n-1, lower.tail=T))) } if (alternative == 'less') { alt <- 'less than' quantiles <- c(qchisq(p=conf.level, df=n-1, lower.tail=F), 0) conf.int <- (n-1) * var(x) / quantiles statistic <- (n-1) * var(x) / null.value p.value <- pchisq(statistic, n-1) } if (alternative == 'greater') { alt <- 'greater than' quantiles <- c(0, qchisq(p=conf.level, df=n-1, lower.tail=T)) conf.int <- (n-1) * var(x) / quantiles statistic <- (n-1) * var(x) / null.value p.value <- pchisq(statistic, n-1, lower.tail=F) } res <- list(conf.int=conf.int, statistic=statistic, df=n-1, alternative=alternative, alt=alt, null.value=null.value, conf.level=conf.level, data.name=deparse(substitute(x)), p.value=p.value, sample.var=var(x)) class(res) <- "vartest" res } var_test_one <- function(varx, nx, alternative, conf.level, null.value) { alpha <- 1 - conf.level # Alternative two.sided if (alternative == 'two.sided') { quantiles <- c(qchisq(p=alpha/2, df=nx-1, lower.tail=F), qchisq(p=1-alpha/2, df=nx-1, lower.tail=F)) conf.int <- (nx-1) * varx / quantiles statistic <- (nx-1) * varx / null.value p.value <- 2 * min(c(pchisq(statistic, nx-1, lower.tail=F), pchisq(statistic, nx-1, lower.tail=T))) } # Alternative less if (alternative == 'less') { quantiles <- c(qchisq(p=conf.level, df=nx-1, lower.tail=T), 0) conf.int <- (nx-1) * varx / quantiles statistic <- (nx-1) * varx / null.value p.value <- pchisq(statistic, nx-1) } # Alternative greater if (alternative == 'greater') { quantiles <- c(Inf, qchisq(p=conf.level, df=nx-1, lower.tail=F)) conf.int <- (nx-1) * varx / quantiles statistic <- (nx-1) * varx / null.value p.value <- pchisq(statistic, nx-1, lower.tail=F) } # To ensure that the output values are in the correct form names(statistic) <- 'X-squared' parameter <- nx - 1 names(parameter) <- 'df' attr(conf.int, 'conf.level') <- conf.level estimate <- varx names(estimate) <- 'variance of x' method <- 'X-squared test for variance' data.name <- paste('varx =', varx, 'and nx =', nx) res <- list(statistic=statistic, parameter=parameter, p.value=p.value, conf.int=conf.int, estimate=estimate, null.value=null.value, alternative=alternative, method=method, data.name=data.name) return(res) } var_test_two <- function(varx, nx, vary, ny, alternative, conf.level, null.value) { alpha <- 1 - conf.level # Alternative two.sided if (alternative == 'two.sided') { quantiles <- c(qf(p=alpha/2, df1=nx-1, df2=ny-1, lower.tail=F), qf(p=1-alpha/2, df1=nx-1, df2=ny-1, lower.tail=F)) conf.int <- (varx / vary) / quantiles statistic <- (varx / vary) / null.value p.value <- 2 * min(c(pf(statistic, nx-1, ny-1, lower.tail=F), pf(statistic, nx-1, ny-1, lower.tail=T))) } # Alternative less if (alternative == 'less') { quantiles <- c(Inf, qf(p=conf.level, df1=nx-1, df2=ny-1, lower.tail=F)) conf.int <- (varx / vary) / quantiles statistic <- (varx / vary) / null.value p.value <- pf(q=statistic, df1=nx-1, df2=ny-1, lower.tail=T) } # Alternative greater if (alternative == 'greater') { quantiles <- c(qf(p=conf.level, df1=nx-1, df2=ny-1, lower.tail=T), 0) conf.int <- (varx / vary) / quantiles statistic <- (varx / vary) / null.value p.value <- pf(q=statistic, df1=nx-1, df2=ny-1, lower.tail=F) } # To ensure that the output values are in the correct form names(statistic) <- 'F' parameter <- c(nx-1, ny-1) names(parameter) <- c('num df', 'denom df') attr(conf.int, 'conf.level') <- conf.level estimate <- varx / vary names(estimate) <- 'ratio of variances' method <- 'F test to compare two variances' data.name <- paste('varx =', varx, ', nx =', nx, ', vary =', vary, 'and ny =', ny) res <- list(statistic=statistic, parameter=parameter, p.value=p.value, conf.int=conf.int, estimate=estimate, null.value=null.value, alternative=alternative, method=method, data.name=data.name) return(res) }