viernes, 30 de diciembre de 2016

Diferentes sumas de cuadrados en ANOVA

Introducción

En los diseños factoriales no balanceados, es decir con un distinto número de réplicas por cada casilla o situación experimental, resulta que las llamadas “sumas de cuadrados” de distintos “tipos” pueden dar diferentes resultados en los contrastes de un ANOVA excepto para el efecto de la interacción. Para calcular la suma de cuadrados o SS de un efecto hay que calcular la diferencia entre la suma de cuadrados de los errores entre los valores observados y sus predicciones mínimo-cuadráticas a partir de un modelo restringido y la suma de cuadrados de los errores entre los valores observados y sus predicciones por un modelo más general. Lo que difiere entre los ``tipos de SS’’ es la elección de los modelos restringido y más general para contrastar un efecto.

Los diseños balanceados son ortogonales y eso disimula completamente las diferencias ya que entonces las distintas hipótesis conducen a los mismos cálculos.

El término “tipos de sumas de cuadrados” fué introducido por SAS y actualmente son cuatro: I, II, III y IV. En SAS y SPSS, entre otros, se utilizan por defecto las SS tipo III, mientras que R utiliza por defecto el tipo I. Esto significa que se pueden obtener resultados distintos cuando se analizan unos mismos datos en función de paquete estadístico utilizado.

Diseños ortogonales

Para que podamos extraer conclusiones definitivas de un experimento, nuestro diseño debería tener dos características fundamentales como son la aleatorización y la ortogonalidad.

La ortogonalidad es una propiedad muy útil ya que nos permite interpretar más fácilmente el efecto de una variable predictora o factor sin que otras variables interfieran. Supongamos que podemos partir la matriz de diseño \(\mathbf{X}\) en dos \(\mathbf{X} = (\mathbf{X}_1 \mathbf{X}_2)\) tales que \(\mathbf{X}_1'\mathbf{X}_2=\mathbf{0}\), entonces \[ \mathbf{Y} = \mathbf{X}\mathbf{\beta} + \mathbf{\epsilon} = \mathbf{X}_1\mathbf{\beta}_1 + \mathbf{X}_2\mathbf{\beta}_2 + \mathbf{\epsilon} \] y \[ \mathbf{X}'\mathbf{X} = \left( \begin{array}{cc} \mathbf{X}_1'\mathbf{X}_1 & \mathbf{X}_1'\mathbf{X}_2 \\ \mathbf{X}_2'\mathbf{X}_1 & \mathbf{X}_2'\mathbf{X}_2 \\ \end{array} \right)= \left( \begin{array}{cc} \mathbf{X}_1'\mathbf{X}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_2'\mathbf{X}_2 \\ \end{array} \right)= \] de manera que \[ \hat{\mathbf{\beta}}_1 = (\mathbf{X}_1'\mathbf{X}_1)^{-1}\mathbf{X}_1'\mathbf{Y} \qquad \hat{\mathbf{\beta}}_2 = (\mathbf{X}_2'\mathbf{X}_2)^{-1}\mathbf{X}_2'\mathbf{Y} \] Así pues la estimación de \(\mathbf{\beta}_1\) será la misma tanto si \(\mathbf{X}_2\) está en el modelo como si no está. De modo que podemos interpretar el efecto de \(\mathbf{X}_1\) sin preocuparnos de \(\mathbf{X}_2\). Aunque en realidad la separación no es perfecta, ya que si queremos contrastar la hipótesis \(H_0: \mathbf{\beta}_1=\mathbf{0}\) necesitaremos la suma de cuadrados residual del modelo (que depende también de \(\mathbf{X}_2\)).

En resumen, la ortogonalidad es una buena propiedad muy aconsejable en cualquier diseño experimental que podamos controlar.

Diferentes sumas de cuadrados

En el artículo Esquemas de códigos para factores quedó claro que en un modelo ANOVA es necesario elegir la codificación de los factores. El objetivo es escribir el ANOVA como un modelo lineal ordinario y utilizar la matriz de diseño resultante como en regresión.

Así, si queremos contrastar el efecto de un factor o de una interacción, podemos eliminar algunas columnas de la matriz de diseño \(\mathbf{X}\) y calcular el incremento del error como una medida del efecto. La razón de esa medida relativa sobre un error global proporciona un valor \(F\) o significación del efecto.

Sin embargo, existen diferentes formas de mantener o eliminar columnas en representación de un efecto y esto ha sido y es motivo de controversia entre estadísticos.

Sumas de tipo I: secuencial

Para ilustrar el análisis de un diseño experimental con dos factores con la función aov de R vamos a utilizar los datos de un experimento de Hunter (1989).

El experimento consistió en quemar una determinada cantidad de combustible y calcular la cantidad de emisiones de CO liberadas. El factor \(A\) es la cantidad de etanol añadida al combustible estándar y el factor \(B\) representa la razón aire/combustible. La variable respuesta es la emisión de monóxido de carbono (CO) en gr/m\(^3\). Se hicieron dos réplicas para cada combinación de los niveles de los factores. Los datos pueden verse en la tabla y se encuentran a nuestra disposición en el data frame COdata del paquete daewr del libro de Lawson (2015).

Las instrucciones de R para analizar los datos son:

library(daewr)
mod1 <- aov(CO ~ Eth * Ratio, data=COdata)
summary(mod1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Eth          2  324.0   162.0   31.36 8.79e-05 ***
## Ratio        2  652.0   326.0   63.10 5.07e-06 ***
## Eth:Ratio    4  678.0   169.5   32.81 2.24e-05 ***
## Residuals    9   46.5     5.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De la tabla ANOVA resultante podemos concluir que los dos efectos y su interacción son significativos como muestran los \(p\)-valores a la derecha del test \(F\).

Veamos que la suma de cuadrados Sum Sq para cada factor es la diferencia de la suma de cuadrados del error cuando el factor se añade al modelo de regresión. Para obtener la suma de cuadrados del error utilizaremos la función deviance que calcula la suma de cuadrados de los residuos \(RSS\) para un modelo lineal dado.

deviance(lm(CO ~ 1, data=COdata)) - deviance(lm(CO ~ Eth, data=COdata))
## [1] 324
deviance(lm(CO ~ Eth, data=COdata)) - deviance(lm(CO ~ Eth + Ratio, data=COdata))
## [1] 652
deviance(lm(CO ~ Eth + Ratio, data=COdata)) - deviance(lm(CO ~ Eth * Ratio, data=COdata))
## [1] 678

Así pues, se trata de ver el efecto como si el factor se añadiera uno a uno, en el orden en el que se ha entrado al modelo. La suma de cuadrados del efecto se puede ver como la reducción en la suma de cuadrados residual \(RSS\) que se obtiene al añadir ese término a un modelo que ya incluye a los términos anteriores.

Una de las propiedades de esta tabla es que la suma de todos los efectos coincide con la suma de cuadrados total \(TSS\):

sum(summary(mod1)[[1]]["Sum Sq"])
## [1] 1700.5
(18-1)*var(COdata$CO)
## [1] 1700.5

Esta propiedad se verifica aunque el modelo sea no balanceado:

mod2 <- aov(CO ~ Eth * Ratio, data=COdata[-18, ])
sum(summary(mod2)[[1]]["Sum Sq"])
## [1] 1467.529
(17-1)*var(COdata$CO[-18])
## [1] 1467.529

Como veremos esta propiedad no se verifica con otros tipos de sumas de cuadrados.


El hecho de que el efecto de un factor dependa del orden en que entra en el modelo puede tener sus ventajas y sus inconvenientes.

Será una ventaja cuando algunos factores (como en el caso jerárquico) deberían eliminarse antes que otros. Por ejemplo, con un número desigual de hombres y de mujeres, el factor “sexo” debería preceder al “sujeto” en un diseño no balanceado.

Como inconveniente: en el caso no balanceado las hipótesis dependen del orden en el que se han especificado los factores. En un ANOVA con dos factores podemos tener dos modelos, uno con \(A\) y \(B\), el otro con \(B\) y \(A\). No solo la suma de cuadrados del efecto \(A\) será distinta en los dos modelos, sino que no podemos predecir si la suma de cuadrados crece o no cuando \(A\) está en segundo lugar y no es el factor principal. En el modelo mod2 tenemos:

##   A*B           Sum Sq B*A           Sum Sq
## 1 Eth         472.6627 Ratio       469.7294
## 2 Ratio       395.3282 Eth         398.2615
## 3 Eth:Ratio   555.0385 Ratio:Eth   555.0385
## 4 Residuals    44.5000 Residuals    44.5000

La falta de invariancia por el orden de entrada de los factores en el modelo limita la utilidad de las sumas de cuadrados de tipo I para contrastar hipótesis en muchos modelos. Por ejemplo, no resulta conveniente en modelos factoriales.

La función model.tables proporciona las estimaciones de la media general \(\hat\mu\), de las medias marginales \(\hat\mu_{i\cdot}\) y \(\hat\mu_{\cdot j}\) y de las medias de cada celda \(\hat\mu_{ij}\).

model.tables(mod1, type="means", se=T)
## Tables of means
## Grand mean
##          
## 72.83333 
## 
##  Eth 
## Eth
##   0.1   0.2   0.3 
## 66.83 75.83 75.83 
## 
##  Ratio 
## Ratio
##   14   15   16 
## 78.5 75.5 64.5 
## 
##  Eth:Ratio 
##      Ratio
## Eth   14   15   16  
##   0.1 64.0 69.5 67.0
##   0.2 79.5 80.5 67.5
##   0.3 92.0 76.5 59.0
## 
## Standard errors for differences of means
##           Eth Ratio Eth:Ratio
##         1.312 1.312     2.273
## replic.     6     6         2

Sin embargo, en el caso no balanceado estas estimaciones no son insesgadas como sí ocurre en el caso anterior.

Consideremos los datos COdata sin la última observación y supongamos que hemos impuesto la codificación contr.sum sobre los dos factores del modelo.

aire/comb
Etanol 14 15 16

0.1
66
62
\(\mu+\alpha_1+\beta_1+\alpha\beta_{11}\)
72
67
\(\mu+\alpha_1+\beta_2+\alpha\beta_{12}\)
68
66
\(\mu+\alpha_1+\beta_3+\alpha\beta_{13}\)

0.2
78
81
\(\mu+\alpha_2+\beta_1+\alpha\beta_{21}\)
80
81
\(\mu+\alpha_2+\beta_2+\alpha\beta_{22}\)
66
69
\(\mu+\alpha_2+\beta_3+\alpha\beta_{23}\)

0.3
90
94
\(\mu+\alpha_3+\beta_1+\alpha\beta_{31}\)
75
78
\(\mu+\alpha_3+\beta_2+\alpha\beta_{32}\)
60

\(\mu+\alpha_3+\beta_3+\alpha\beta_{33}\)


En la tabla se ven las observaciones de las dos réplicas por casilla, excepto en la última. También para cada casilla se muestra la función paramétrica que se estima con cada réplica.

La media marginal para las dos primeras columnas (también para las dos primeras filas) es un estimador insesgado del efecto. Por ejemplo, para la primera columna resulta que \[ \frac{1}{6}((2\mu+2\alpha_1+2\beta_1+2\alpha\beta_{11})+ (2\mu+2\alpha_2+2\beta_1+2\alpha\beta_{21}) + (2\mu+2\alpha_3+2\beta_1+2\alpha\beta_{31})) = \mu + \beta_1 \] ya que \(\alpha_1+\alpha_2+\alpha_3=0\) y \(\alpha\beta_{11}+ \alpha\beta_{21}+\alpha\beta_{31}=0\).

Lo mismo ocurre con la segunda columna \(\mu+\beta_2\) y con las dos primeras filas \(\mu+\alpha_1\) y \(\mu+\alpha_2\).

Sin embargo, en la tercera columna (también en la tercera fila) tenemos: \[ \frac{1}{5}((2\mu+2\alpha_1+2\beta_3+2\alpha\beta_{13})+ (2\mu+2\alpha_2+2\beta_3+2\alpha\beta_{23}) + (\mu+\alpha_3+\beta_3+2\alpha\beta_{33})) = \mu + \beta_3 + \gamma \] donde \(\gamma=(2\alpha_1+2\alpha_2+\alpha_3)/5 + (2\alpha\beta_{13}+2\alpha\beta_{23}+\alpha\beta_{33})/5\). Así, la media de la columna no es un estimador insesgado de \(\mu+\beta_3\) y la función model.tables sobre el modelo mod2 proporciona medias sesgadas de los parámetros \(\mu_{\cdot 3}\) y \(\mu_{3\cdot}\).

Sumas de tipo II: jerárquico o parcialmente secuencial

La suma de cuadrados del tipo II para el factor \(A\) en un diseño de dos factores con interacción es la diferencia entre la suma de cuadrados de los residuos del modelo reducido y ~ B y el modelo y ~ A+B. Del mismo modo, la suma de cuadrados del tipo II para el factor \(B\) es la diferencia entre la suma de cuadrados de los residuos del modelo reducido y ~ A y el modelo y ~ A+B. Para la interacción es la diferencia entre el modelo y ~ A + B y el modelo y ~ A + B + A:B.

En general, la suma de cuadrados de un efecto es la reducción en la suma de cuadrados del error al añadir ese término al modelo que ya contiene todos los otros efectos. La interacción solo aparece si están presentes todos los efectos de esa interacción. Por ejemplo, la suma de cuadrados del efecto principal del factor \(A\) no se debe ajustar por ninguna interacción que contenga \(A\).


En los siguientes cálculos vamos a utilizar la función lm con la opción contr.sum.

modB <- lm(CO ~ Ratio, data=COdata[-18,], contrasts=list(Ratio=contr.sum))
mods <- lm(CO ~ Eth + Ratio, data=COdata[-18,], 
                             contrasts=list(Eth=contr.sum, Ratio=contr.sum))
deviance(modB) - deviance(mods)
## [1] 398.2615
modA <- lm(CO ~ Eth, data=COdata[-18,], contrasts=list(Eth=contr.sum))
deviance(modA) - deviance(mods)
## [1] 395.3282
modc <- lm(CO ~ Eth * Ratio, data=COdata[-18,], 
                             contrasts=list(Eth=contr.sum, Ratio=contr.sum))
deviance(mods) - deviance(modc)
## [1] 555.0385

Para calcular las sumas de cuadrados de tipo II para las hipótesis sobre los efectos principales vamos a utilizar la función Anova del paquete car.

library(car)
Anova(modc, type="II")
## Anova Table (Type II tests)
## 
## Response: CO
##           Sum Sq Df F value    Pr(>F)    
## Eth       398.26  2  35.799 0.0001020 ***
## Ratio     395.33  2  35.535 0.0001048 ***
## Eth:Ratio 555.04  4  24.945 0.0001427 ***
## Residuals  44.50  8                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cuando el modelo es balanceado, las sumas de cuadrados de tipo I coinciden con las de tipo II.


Entre las ventajas de las sumas de cuadrados de tipo II tenemos:

  • Son invariantes respecto al orden de entrada de los efectos en el modelo.
  • Apropiadas para construir modelos y la elección natural en regresión.
  • Más potencia cuando no hay interacción.


Para obtener las medias insesgadas de los efectos columna \(\mu+\beta_1\), \(\mu+\beta_2\) y \(\mu+\beta_3\), deberemos calcular las medias ajustadas. Dichas medias ajustadas, que también se llaman medias mínimo-cuadráticas (lsmeans), se calculan con las medias marginales de las predicciones medias para cada casilla.

El siguiente código de R construye un data frame con todas las combinaciones para cada casilla y calcula las predicciones con el modelo completo en esas casillas. Luego calcula las medias marginales con la función tapply.

p <- data.frame( expand.grid(Eth=c(.1,.2,.3), Ratio=c(14,15,16)))
p[] <- lapply(p, factor)
yhat <- predict(modc, p)
p <- cbind(yhat , p)
with(p, tapply(yhat, Ratio, mean))
##       14       15       16 
## 78.50000 75.50000 64.83333

Podemos ver que las medias para las dos primeras columnas son las mismas que en el apartado anterior y coinciden con las medias aritméticas de las respuestas. Sin embargo, en la tercera columna hemos obtenido un estimador distinto y más acertado para \(\mu+\beta_3\).

El paquete lsmeans nos ayudará a calcular automáticamente las medias ajustadas o lsmeans.

library(lsmeans)
lsmeans(modc, ~ Eth)
##  Eth   lsmean        SE df lower.CL upper.CL
##  0.1 66.83333 0.9628517  8 64.61299 69.05367
##  0.2 75.83333 0.9628517  8 73.61299 78.05367
##  0.3 76.16667 1.1118053  8 73.60284 78.73049
## 
## Results are averaged over the levels of: Ratio 
## Confidence level used: 0.95
lsmeans(modc, ~ Ratio)
##  Ratio   lsmean        SE df lower.CL upper.CL
##  14    78.50000 0.9628517  8 76.27966 80.72034
##  15    75.50000 0.9628517  8 73.27966 77.72034
##  16    64.83333 1.1118053  8 62.26951 67.39716
## 
## Results are averaged over the levels of: Eth 
## Confidence level used: 0.95

Sumas de tipo III: marginal u ortogonal

Las sumas de cuadrados tipo III son las que se obtienen cuando cada variable se añade la última en el modelo. Es decir, el efecto de cada variable es el que se calcula después de que todos los otros factores ya se han tenido en cuenta. En consecuencia el resultado para cada término es equivalente al que se obtiene con el tipo I cuando el término entra el último en el modelo.

En el caso de dos factores con interacción, la suma de cuadrados del efecto del factor \(A\) es la diferencia entre las sumas de cuadrados del error del modelo y ~ B + A:B y del modelo completo y ~ A + B + A:B.

Anova(modc, type="III")
## Anova Table (Type III tests)
## 
## Response: CO
##             Sum Sq Df   F value    Pr(>F)    
## (Intercept)  86198  1 15496.351 1.939e-14 ***
## Eth            319  2    28.715 0.0002235 ***
## Ratio          511  2    45.973 4.105e-05 ***
## Eth:Ratio      555  4    24.945 0.0001427 ***
## Residuals       44  8                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Entre las ventajas de este tipo de sumas de cuadrados hay que destacar que no dependen del tamaño muestral de cada grupo y se pueden utilizar en modelos no balanceados, con distinto número de observaciones en cada grupo. Cuando no hay casillas vacías, las medias en cada subpoblación son las medias mínimo-cuadráticas o lsmeans, que son los mejores estimadores lineales insesgados (BLUE) para las medias marginales del diseño.

Por contra, con este tipo no se puede contrastar los efectos principales en presencia de interacciones.

Tampoco resultan apropiadas cuando hay casillas vacías y por ello se han definido las sumas de cuadrados de tipo IV, que son una variación de las de tipo III especialmente desarrolladas para esta situación.

Resumen

Consideremos un modelo de dos factores con interacción en el que los términos aparecen en el orden A, B, AB. Sea \(R(\cdot)\) la suma de cuadrados de los residuos o del error para un modelo. Por ejemplo, el modelo completo tiene una \(RSS=R(A,B,AB)\). El modelo con únicamente la media general tiene una suma de cuadrados del error que escribiremos \(R(1)\) y para el modelo con únicamente el factor A es \(R(A)\). Pues bien, con esta notación, en la tabla tenemos un resumen de los tres tipos de sumas de cuadrados que hemos explicado en los apartados anteriores.

Término SS tipo I SS tipo II SS tipo III
\(A\) \(R(1) -R(A)\) \(R(B) - R(A,B)\) \(R(B,AB) - R(A,B,AB)\)
\(B\) \(R(A) -R(A,B)\) \(R(A) - R(A,B)\) \(R(A,AB) - R(A,B,AB)\)
\(AB\) \(R(A,B) -R(A,B,AB)\) \(R(A,B) - R(A,B,AB)\) \(R(A,B) - R(A,B,AB)\)

El tipo de sumas de cuadrados solo influye en los cálculos cuando el diseño es no balanceado, ya que para diseños ortogonales no importa qué suma de cuadrados utilicemos pues todas son esencialmente iguales. Así, si es posible, los diseños balanceados son los preferidos en el análisis de grupos.

Por ejemplo, las sumas de cuadrados de tipo II para el diseño con los datos completos de COdata coincide con los resultados de la función aov de la página .

Anova(mod1, type="II")
## Anova Table (Type II tests)
## 
## Response: CO
##           Sum Sq Df F value    Pr(>F)    
## Eth        324.0  2  31.355 8.790e-05 ***
## Ratio      652.0  2  63.097 5.067e-06 ***
## Eth:Ratio  678.0  4  32.806 2.240e-05 ***
## Residuals   46.5  9                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El análisis del efecto de la interacción se hace con las mismas sumas de cuadrados en los tres tipos.

Cuando hay casillas vacías se recomienda el tipo IV.

No hay consenso sobre qué tipo de sumas de cuadrados debemos utilizar en diseños no balanceados, pero la mayoría de estadísticos recomiendan el tipo III, que es la elección por defecto en muchos programas como SAS, SPSS, JMP, Minitab, Stata, Statista, Systat y Unistat, mientras que otros como R, S-Plus, Genstat y Mathematica usan el tipo I. Sin embargo, Langsrud (2003) argumenta que el tipo II es preferible si consideramos la potencia de los tipos II y III.

Referencias

F. Carmona, Modelos lineales, Publicacions UB, 2005.

J.J. Faraway, Linear Models with R, Chapman & Hall/CRC, 2014.

J. Fox, An R and S-Plus companion to applied regression, Sage, 2002.

J.S. Hunter, Let’s All Beware the Latin Square, Quality Engineering, Vol. 1, 1989, pp. 453-465.

O. Langsrud, ANOVA for Unbalanced Data: Use Type II Instead of Type III Sums of Squares, Statistics and Computing, 2003, 13, 163-167.

J. Lawson, Design and Analysis of Experiments with R, Chapman & Hall/CRC, 2015.

D.C Montgomery, Design and Analysis of Experiments, John Wiley & Sons, 2008.

1 comentario:

  1. Muy interesante el blog. Quiero aprovechar para hacer una pregunta sobre R y analítica.
    Estoy analizando las ventas de un retailer. Estas ventas están relacionadas con una serie de variables: var1, var2, var3, var4. La mayoría de ellas son continuas.

    Para analizar la dependencia de las ventas respecto las variables creo una regresión lineal:

    rg<-lm(sales ~ var1 + var2 + var3 + var4, data=sales_2017)
    summary(rg)


    Posteriormente quiero analizar el peso de cada variable y para ello uso (paquete caret):

    varImp(rg, scale = FALSE)
    rsimp <- varImp(rg, scale = FALSE)
    plot(rsimp)

    Es este un buen método para analizar el peso e importancia de las variables?, alguna recomendación?

    Juan

    ResponderEliminar