Cómo interpretar un gráfico QQ

Question

Estoy trabajando con un pequeño conjunto de datos (21 observaciones) y tengo el siguiente gráfico QQ normal en R:

Viendo que el gráfico no admite la normalidad, ¿qué podría inferir sobre la distribución subyacente? Me parece que una distribución más sesgada hacia la derecha se ajustaría mejor, ¿es así? Además, ¿qué otras conclusiones podemos sacar de los datos?

Answer 1

Si los valores se encuentran a lo largo de una línea, la distribución tiene la misma forma (hasta la ubicación y la escala) que la distribución teórica que hemos supuesto.

Comportamiento local : Al observar los valores muestrales ordenados en el eje Y y los cuantiles (aproximados) esperados en el eje X, podemos identificar a partir de cómo los valores de alguna sección del gráfico difieren localmente de una tendencia lineal general, viendo si los valores están más o menos concentrados de lo que supondría la distribución teórica en esa sección de un gráfico:

Como vemos, los puntos menos concentrados aumentan más y los puntos más concentrados de lo que se supone aumentan menos rápidamente de lo que sugeriría una relación lineal global, y en los casos extremos corresponden a una brecha en la densidad de la muestra (se muestra como un salto casi vertical) o un pico de valores constantes (valores alineados horizontalmente). Esto nos permite detectar una cola pesada o una cola ligera y, por tanto, una asimetría mayor o menor que la de la distribución teórica, etc.

Apariencia general:

Este es el aspecto de los gráficos QQ (para determinadas opciones de distribución) de media :

Pero la aleatoriedad tiende a oscurecer las cosas, especialmente con muestras pequeñas:

Tenga en cuenta que en $n=21$ los resultados pueden ser mucho más variables de lo que se muestra allí - generé varios conjuntos de seis gráficos y elegí un conjunto "bonito" en el que se podía ver la forma en los seis gráficos al mismo tiempo. A veces las relaciones rectas parecen curvas, las relaciones curvas parecen rectas, las colas pesadas parecen torcidas, etc. - con muestras tan pequeñas, a menudo la situación puede ser mucho menos clara:

Es posible discernir más características que esas (como la discreción, por ejemplo), pero con $n=21$ Incluso estos rasgos básicos pueden ser difíciles de detectar; no debemos tratar de "sobreinterpretar" cada pequeño meneo. A medida que aumenta el tamaño de las muestras, en general los gráficos se "estabilizan" y las características se interpretan más claramente en lugar de representar ruido. En el caso de algunas distribuciones con colas muy pesadas, el raro valor atípico puede impedir que la imagen se estabilice bien incluso con tamaños de muestra bastante grandes].

También puede encontrar la sugerencia aquí útil cuando se trata de decidir cuánto debe preocuparse por una determinada cantidad de curvatura o ondulación.

Una guía más adecuada para la interpretación en general incluiría también visualizaciones con tamaños de muestra más pequeños y más grandes.

Answer 2

Hice una aplicación brillante para ayudar a interpretar la trama normal de QQ. Prueba este enlace.

En esta aplicación, puede ajustar la asimetría, la cola (curtosis) y la modalidad de los datos y puede ver cómo cambian el histograma y el gráfico QQ. A la inversa, se puede utilizar de manera que dado el patrón de la parcela QQ, a continuación, comprobar cómo la asimetría, etc debe ser.

Para más detalles, consulte la documentación correspondiente.

---

Me di cuenta de que no tengo suficiente espacio libre para proporcionar esta aplicación en línea. Como petición, voy a proporcionar los tres trozos de código: sample.R , server.R y ui.R aquí. Aquellos que estén interesados en ejecutar esta aplicación pueden simplemente cargar estos archivos en Rstudio y luego ejecutarlos en su propio PC.

El sample.R archivo:

# Compute the positive part of a real number x, which is $\max(x, 0)$.
positive_part <- function(x) {ifelse(x > 0, x, 0)}

# This function generates n data points from some unimodal population.
# Input: ----------------------------------------------------
# n: sample size;
# mu: the mode of the population, default value is 0.
# skewness: the parameter that reflects the skewness of the distribution, note it is not
#           the exact skewness defined in statistics textbook, the default value is 0.
# tailedness: the parameter that reflects the tailedness of the distribution, note it is
#             not the exact kurtosis defined in textbook, the default value is 0.

# When all arguments take their default values, the data will be generated from standard 
# normal distribution.

random_sample <- function(n, mu = 0, skewness = 0, tailedness = 0){
  sigma = 1

  # The sampling scheme resembles the rejection sampling. For each step, an initial data point
  # was proposed, and it will be rejected or accepted based on the weights determined by the
  # skewness and tailedness of input. 
  reject_skewness <- function(x){
      scale = 1
      # if `skewness` > 0 (means data are right-skewed), then small values of x will be rejected
      # with higher probability.
      l <- exp(-scale * skewness * x)
      l/(1 + l)
  }

  reject_tailedness <- function(x){
      scale = 1
      # if `tailedness` < 0 (means data are lightly-tailed), then big values of x will be rejected with
      # higher probability.
      l <- exp(-scale * tailedness * abs(x))
      l/(1 + l)
  }

  # w is another layer option to control the tailedness, the higher the w is, the data will be
  # more heavily-tailed. 
  w = positive_part((1 - exp(-0.5 * tailedness)))/(1 + exp(-0.5 * tailedness))

  filter <- function(x){
    # The proposed data points will be accepted only if it satified the following condition, 
    # in which way we controlled the skewness and tailedness of data. (For example, the 
    # proposed data point will be rejected more frequently if it has higher skewness or
    # tailedness.)
    accept <- runif(length(x)) > reject_tailedness(x) * reject_skewness(x)
    x[accept]
  }

  result <- filter(mu + sigma * ((1 - w) * rnorm(n) + w * rt(n, 5)))
  # Keep generating data points until the length of data vector reaches n.
  while (length(result) < n) {
    result <- c(result, filter(mu + sigma * ((1 - w) * rnorm(n) + w * rt(n, 5))))
  }
  result[1:n]
}

multimodal <- function(n, Mu, skewness = 0, tailedness = 0) {
  # Deal with the bimodal case.
  mumu <- as.numeric(Mu %*% rmultinom(n, 1, rep(1, length(Mu))))
  mumu + random_sample(n, skewness = skewness, tailedness = tailedness)
}

El server.R archivo:

library(shiny)
# Need 'ggplot2' package to get a better aesthetic effect.
library(ggplot2)

# The 'sample.R' source code is used to generate data to be plotted, based on the input skewness, 
# tailedness and modality. For more information, see the source code in 'sample.R' code.
source("sample.R")

shinyServer(function(input, output) {
  # We generate 10000 data points from the distribution which reflects the specification of skewness,
  # tailedness and modality. 
  n = 10000

  # 'scale' is a parameter that controls the skewness and tailedness.
  scale = 1000

  # The `reactive` function is a trick to accelerate the app, which enables us only generate the data
  # once to plot two plots. The generated sample was stored in the `data` object to be called later.
  data <- reactive({
    # For `Unimodal` choice, we fix the mode at 0.
    if (input$modality == "Unimodal") {mu = 0}

    # For `Bimodal` choice, we fix the two modes at -2 and 2.
    if (input$modality == "Bimodal") {mu = c(-2, 2)}

    # Details will be explained in `sample.R` file.
    sample1 <- multimodal(n, mu, skewness = scale * input$skewness, tailedness = scale * input$kurtosis)
    data.frame(x = sample1)})

  output$histogram <- renderPlot({
    # Plot the histogram.
    ggplot(data(), aes(x = x)) + 
      geom_histogram(aes(y = ..density..), binwidth = .5, colour = "black", fill = "white") + 
      xlim(-6, 6) +
      # Overlay the density curve.
      geom_density(alpha = .5, fill = "blue") + ggtitle("Histogram of Data") + 
      theme(plot.title = element_text(lineheight = .8, face = "bold"))
  })

  output$qqplot <- renderPlot({
    # Plot the QQ plot.
    ggplot(data(), aes(sample = x)) + stat_qq() + ggtitle("QQplot of Data") + 
      theme(plot.title = element_text(lineheight=.8, face = "bold"))
    })
})

Por último, el ui.R archivo:

library(shiny)

# Define UI for application that helps students interpret the pattern of (normal) QQ plots. 
# By using this app, we can show students the different patterns of QQ plots (and the histograms,
# for completeness) for different type of data distributions. For example, left skewed heavy tailed
# data, etc. 

# This app can be (and is encouraged to be) used in a reversed way, namely, show the QQ plot to the 
# students first, then tell them based on the pattern of the QQ plot, the data is right skewed, bimodal,
# heavy-tailed, etc.

shinyUI(fluidPage(
  # Application title
  titlePanel("Interpreting Normal QQ Plots"),

  sidebarLayout(
    sidebarPanel(
      # The first slider can control the skewness of input data. "-1" indicates the most left-skewed 
      # case while "1" indicates the most right-skewed case.
      sliderInput("skewness", "Skewness", min = -1, max = 1, value = 0, step = 0.1, ticks = FALSE),

      # The second slider can control the skewness of input data. "-1" indicates the most light tail
      # case while "1" indicates the most heavy tail case.
      sliderInput("kurtosis", "Tailedness", min = -1, max = 1, value = 0, step = 0.1, ticks = FALSE),

      # This selectbox allows user to choose the number of modes of data, two options are provided:
      # "Unimodal" and "Bimodal".
      selectInput("modality", label = "Modality", 
                  choices = c("Unimodal" = "Unimodal", "Bimodal" = "Bimodal"),
                  selected = "Unimodal"),
      br(),
      # The following helper information will be shown on the user interface to give necessary
      # information to help users understand sliders.
      helpText(p("The skewness of data is controlled by moving the", strong("Skewness"), "slider,", 
               "the left side means left skewed while the right side means right skewed."), 
               p("The tailedness of data is controlled by moving the", strong("Tailedness"), "slider,", 
                 "the left side means light tailed while the right side means heavy tailedd."),
               p("The modality of data is controlledy by selecting the modality from", strong("Modality"),
                 "select box.")
               )
  ),

  # The main panel outputs two plots. One plot is the histogram of data (with the nonparamteric density
  # curve overlaid), to get a better visualization, we restricted the range of x-axis to -6 to 6 so 
  # that part of the data will not be shown when heavy-tailed input is chosen. The other plot is the 
  # QQ plot of data, as convention, the x-axis is the theoretical quantiles for standard normal distri-
  # bution and the y-axis is the sample quantiles of data. 
  mainPanel(
    plotOutput("histogram"),
    plotOutput("qqplot")
  )
)
)
)

Answer 3

Una explicación muy útil (e intuitiva) es la que da el prof. Philippe Rigollet en el curso MOOC del MIT: 18.650 Statistics for Applications, otoño de 2016 - ver vídeo a los 45 min.

https://www.youtube.com/watch?v=vMaKx9fmJHE

He copiado burdamente su diagrama, que guardo en mis apuntes porque me resulta muy útil.

En el ejemplo 1, en el diagrama superior izquierdo, vemos que en la cola derecha el cuantil empírico (o muestral) es menor que el cuantil teórico

Qe < Qt

Esto puede interpretarse mediante las funciones de densidad de probabilidad. Para la misma $\alpha$ el cuantil empírico está a la izquierda del cuantil teórico, lo que significa que la cola derecha de la distribución empírica es más "ligera" que la cola derecha de la distribución teórica, es decir, cae más rápido hacia valores cercanos a cero.

Answer 4

Ya que este hilo ha sido considerado como un post definitivo de StackExchange sobre "cómo interpretar el gráfico normal q-q", me gustaría señalar a los lectores una bonita y precisa relación matemática entre el gráfico normal q-q y el estadístico de exceso de curtosis.

Aquí está:

https://stats.stackexchange.com/a/354076/102879

A continuación se ofrece un breve (y demasiado simplificado) resumen (véase el enlace para ver los enunciados matemáticos más precisos): En realidad se puede ver el exceso de curtosis en el gráfico normal q-q como la distancia media entre los cuantiles de los datos y los correspondientes cuantiles normales teóricos, ponderada por la distancia de los datos a la media. Así, cuando los valores absolutos en las colas del gráfico q-q se desvían generalmente de los valores normales esperados en las direcciones extremas, se tiene un exceso de curtosis positivo.

Como la curtosis es la media de estas desviaciones ponderadas por las distancias a la media, los valores cercanos al centro del gráfico q-q tienen poco impacto en la curtosis. Por lo tanto, el exceso de curtosis no está relacionado con el centro de la distribución, donde está el "pico". Más bien, el exceso de curtosis viene determinado casi por completo por la comparación de las colas de la distribución de datos con la distribución normal.