Para introducir en un modelo lineal una variable dicotómica (dos valores o niveles) es suficiente con considerar un vector de ceros y unos que representan cada uno de los niveles. Ese vector representará a la variable en la matriz de diseño del modelo lineal.
Sin embargo, cuando tenemos un factor con más de dos niveles, la cosa se complica.
Si tenemos un factor con \(k\) niveles, sea \(\mathbf{B} = (b_{ij})\) una matriz \(n\times k\) de ceros y unos donde \(b_{ij} = 1\) si la muestra \(i\) pertence a la clase \(j\) y cero en otro caso. Esta matriz \(\mathbf{B}\) se debería incorporar a la matriz de diseño del modelo, pero hay un inconveniente. La suma de las columnas de \(\mathbf{B}\) es uno en todas las muestras y, como todo modelo incorpora un vector de unos que representa la media general o la intercepción, la matriz de diseño no será de rango máximo y no podremos estimar todos los parámetros.
Una solución es eliminar de la matriz de diseño la columna de unos que representa a la media general, pero eso no funciona si tenemos más de un factor. La solución pasa por considerar una matriz de rango \(k-1\) que represente al factor. Para ello podemos simplemente eliminar una columna la matriz \(\mathbf{B}\) o utilizar una matriz de \(k-1\) columnas que evite la colinealidad. La elección se debería basar en la facilidad de cálculo y en la viabilidad de la interpretación de los resultados.
El código se fija con una matriz \(\mathbf{C}\) de dimensión \(k\times (k-1)\) llamada matriz de contraste. La contribución del factor a la matriz de diseño es exactamente \(\mathbf{B}\mathbf{C}\). La matriz \(\mathbf{B}\mathbf{C}\) será una submatriz de \(k-1\) columnas de la matriz de diseño. Lógicamente la matriz de diseño puede incorporar la columna de unos para la media general y seguramente otras columnas para otras variables o factores del modelo.
Ejemplo Supongamos que el vector \((\mu,\alpha_1,\alpha_2,\alpha_3,\alpha_4)'\) representa los parámetros de un modelo de un factor con cuatro niveles. Como sabemos, el problema es que tenemos cinco parámetros y sólo cuatro poblaciones.
La matriz de diseño reducida (con una única fila para cada situación experimental) es: \[ (\mathbf{1},\mathbf{B}_r)= \begin{array}{ccccc} \mu & \alpha_1 & \alpha_2 & \alpha_3 & \alpha_4 \\ 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 1 \end{array} \] Esta matriz de diseño es de rango 4 lo que impide estimar los parámetros de forma única. Así pues, es preciso (1) añadir alguna restricción sobre los parámetros o (2) reparametrizar con una matriz de contraste adecuada.
La matriz de contraste \(\mathbf{C}\) será una matriz \(4\times (4-1)\) de rango \(4-1\) donde cada columna suma 0 \[ \begin{pmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{pmatrix} =\mathbf{C} \begin{pmatrix} \gamma_1 \\ \gamma_2 \\ \gamma_3 \end{pmatrix} \] y la matriz de diseño reducida con los nuevos parámetros \((\mu,\gamma_1,\gamma_2,\gamma_3)\) es \((\mathbf{1},\mathbf{B}_r\mathbf{C})\).
Por otra parte, si \({\boldsymbol\mu} = (\mu_1,\mu_2,\mu_3,\mu_4)\) representa las medias poblacionales de cada nivel, entonces \[ {\boldsymbol\mu} = (\mathbf{1},\mathbf{C}) \begin{pmatrix} \mu \\ \gamma_1 \\ \gamma_2 \\ \gamma_3 \end{pmatrix} \] y como \((\mathbf{1},\mathbf{C})\) es una matriz no singular, podemos hallar los nuevos parámetros como combinación lineal de las medias: \[ \begin{pmatrix} \mu \\ \gamma_1 \\ \gamma_2 \\ \gamma_3 \end{pmatrix} = (\mathbf{1},\mathbf{C})^{-1}{\boldsymbol\mu} \]
El truco consiste en plantear una matriz \(\mathbf{C}\) de forma que el nuevo vector de parámetros, sin la media general, represente un conjunto de contrastes interesantes sobre las medias de cada nivel.
Una forma sencilla de proceder, aunque no la única, es hacer que las columnas de \(\mathbf{C}\) sean ortogonales entre sí. De este modo las filas de \((\mathbf{1},\mathbf{C})^{-1}\) serán proporcionales a las correspondientes columnas de \((\mathbf{1},\mathbf{C})\) y se podrá codificar directamente los contrastes de interés sobre las medias en las columnas de \(\mathbf{C}\).
Nada de lo propuesto requiere que el factor tenga igual número de observaciones en cada uno de los niveles, pero si el número de réplicas es el mismo, entonces no sólo las columnas de \(\mathbf{C}\) son ortogonales, sino que también los serán las columnas de la matriz de diseño \(\mathbf{X}\) construida con una tal \(\mathbf{C}\) ortogonal. De esta forma podremos partir la suma de cuadrados de la regresión para ese modelo en sus componentes con un grado de libertad para cada contraste.
Ejemplo Vamos a considerar un conjunto de datos a partir de un estudio (los datos provienen de Box, Hunter and Hunter, 1978) sobre el tiempo de coagulación de la sangre. En este estudio 24 animales se asignaron al azar a cuatro dietas y las muestras se tomaron en un orden aleatorio.
require(faraway)
data(coagulation)
str(coagulation)## 'data.frame': 24 obs. of 2 variables:
## $ coag: num 62 60 63 59 63 67 71 64 65 66 ...
## $ diet: Factor w/ 4 levels "A","B","C","D": 1 1 1 1 2 2 2 2 2 2 ...Nos vamos a centrar en el ANOVA de la variable categórica diet que tiene cuatro niveles y utilizaremos la variable coag como respuesta. Podemos ver dos representaciones de estos datos en la figura.
Las medias de las puntuaciones de la variable coag para cada dieta por separado son:
tapply(coagulation$coag, coagulation$diet, mean)## A B C D
## 61 66 68 61La pregunta es: ¿Son significativamente distintas?
Codificación dummy
La matriz de contraste más sencilla para un factor con 4 niveles es:
contr.treatment(4)## 2 3 4
## 1 0 0 0
## 2 1 0 0
## 3 0 1 0
## 4 0 0 1Se trata de considerar el primer nivel como el nivel estándar, basal o control y comparar los otros niveles con ese. Esta es la matriz de contraste por defecto en R.
Ejemplo Las estimaciones de los parámetros del modelo lineal son:
summary(lm(coag ~ diet, coagulation))##
## Call:
## lm(formula = coag ~ diet, data = coagulation)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.00 -1.25 0.00 1.25 5.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.100e+01 1.183e+00 51.554 < 2e-16 ***
## dietB 5.000e+00 1.528e+00 3.273 0.003803 **
## dietC 7.000e+00 1.528e+00 4.583 0.000181 ***
## dietD 2.991e-15 1.449e+00 0.000 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.366 on 20 degrees of freedom
## Multiple R-squared: 0.6706, Adjusted R-squared: 0.6212
## F-statistic: 13.57 on 3 and 20 DF, p-value: 4.658e-05Con esta codificación la matriz de diseño reducida es:
(xr <- cbind(rep(1,4),contr.treatment(4)))## 2 3 4
## 1 1 0 0 0
## 2 1 1 0 0
## 3 1 0 1 0
## 4 1 0 0 1y su inversa
solve(xr)## 1 2 3 4
## 1 0 0 0
## 2 -1 1 0 0
## 3 -1 0 1 0
## 4 -1 0 0 1con lo que queda claro lo que representa cada uno de los contrastes individuales en el summary anterior: \[
\left(
\begin{array}{rrrr}
1 & 0 & 0 & 0 \\
-1 & 1 & 0 & 0 \\
-1 & 0 & 1 & 0 \\
-1 & 0 & 0 & 1
\end{array}
\right)
\begin{pmatrix}
\mu_1 \\
\mu_2 \\
\mu_3 \\
\mu_4
\end{pmatrix} =
\begin{pmatrix}
\mu_1 \\
\mu_2-\mu_1 \\
\mu_3-\mu_1 \\
\mu_4-\mu_1
\end{pmatrix}
\]
En SAS el nivel base es el último de forma que la matriz de contraste es:
contr.SAS(4)## 1 2 3
## 1 1 0 0
## 2 0 1 0
## 3 0 0 1
## 4 0 0 0que es equivalente a la que se obtiene con el código:
contr.treatment(4, base=4)ya que el parámetro base sirve para fijar el nivel basal o control sin necesidad de reordenar los niveles del factor.
Ejemplo Si queremos cambiar el tipo de codificación para un factor específico, lo podemos hacer directamente sobre ese factor así:
contrasts(coagulation$diet) <- contr.treatment(4, base=4)Ahora el modelo lineal para el contraste proporciona el siguiente resultado:
summary(lm(coag ~ diet, coagulation))##
## Call:
## lm(formula = coag ~ diet, data = coagulation)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.00 -1.25 0.00 1.25 5.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.100e+01 8.367e-01 72.909 < 2e-16 ***
## diet1 4.378e-15 1.449e+00 0.000 1.000000
## diet2 5.000e+00 1.278e+00 3.912 0.000864 ***
## diet3 7.000e+00 1.278e+00 5.477 2.32e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.366 on 20 degrees of freedom
## Multiple R-squared: 0.6706, Adjusted R-squared: 0.6212
## F-statistic: 13.57 on 3 and 20 DF, p-value: 4.658e-05Código de desviación
Este código compara la media de cada nivel con la media de las medias de la variable respuesta.
Este esquema se consigue en R así:
contr.sum(4)## [,1] [,2] [,3]
## 1 1 0 0
## 2 0 1 0
## 3 0 0 1
## 4 -1 -1 -1Como vemos en este caso, los códigos asignan 1 al nivel 1 para la primera comparación (que compararemos con todos los otros niveles), un 1 también al nivel 2 para la segunda comparación (ya que el nivel 2 se comparará con todos los otros) y un 1 en el nivel 3 para la tercera comparación (ya que también se comparará con todos los otros). En cambio, asignamos un -1 al nivel 4 para las tres comparaciones (ya que es el nivel que no compararemos con los demás).
Creo que se puede entender mejor si calculamos la matriz \((\mathbf{1},\mathbf{C})^{-1}\)
solve(cbind(rep(1,4),contr.sum(4)))## 1 2 3 4
## [1,] 0.25 0.25 0.25 0.25
## [2,] 0.75 -0.25 -0.25 -0.25
## [3,] -0.25 0.75 -0.25 -0.25
## [4,] -0.25 -0.25 0.75 -0.25que significa \[ \frac 14 \left( \begin{array}{rrrr} 1 & 1 & 1 & 1 \\ 3 & -1 & -1 & -1 \\ -1 & 3 & -1 & -1 \\ -1 & -1 & 3 & -1 \end{array} \right) \begin{pmatrix} \mu_1 \\ \mu_2 \\ \mu_3 \\ \mu_4 \end{pmatrix} = \begin{pmatrix} \bar\mu \\ \mu_1-\bar\mu \\ \mu_2-\bar\mu \\ \mu_3-\bar\mu \end{pmatrix} \] donde \(\bar\mu=(\mu_1+\mu_2+\mu_3+\mu_4)/4\).
Ejemplo Con los datos del ejemplo del tiempo de coagulación tenemos:
contrasts(coagulation$diet) <- contr.sum(4)
summary(lm(coag ~ diet, coagulation))##
## Call:
## lm(formula = coag ~ diet, data = coagulation)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.00 -1.25 0.00 1.25 5.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.0000 0.4979 128.537 < 2e-16 ***
## diet1 -3.0000 0.9736 -3.081 0.005889 **
## diet2 2.0000 0.8453 2.366 0.028195 *
## diet3 4.0000 0.8453 4.732 0.000128 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.366 on 20 degrees of freedom
## Multiple R-squared: 0.6706, Adjusted R-squared: 0.6212
## F-statistic: 13.57 on 3 and 20 DF, p-value: 4.658e-05Observemos que efectivamente la estimación de \(\bar\mu\) es \((61+66+68+61)/4=64\) y que este valor no es la media general \(\hat\mu\) de todos los datos. El estimador para el primer contraste (segunda fila del summary) es \(61 - 64\), el segundo contraste es \(66-64\) y el tercero \(68-64\). Cada una de las filas del summary contrasta si ese valor es cero.
Código simple
Este código muy parecido al dummy en el que se compara cada nivel con un nivel de referencia. La diferencia está en la constante o intercept. En la codificación dummy la constante representa la media del nivel de referencia, mientras que con la codificación simple la constante es la media de todas las celdas \(\bar\mu\).
Este código hay que definirlo explícitamente así:
contr.simple <- function(k) contr.treatment(k) - matrix(rep(1/k, k*(k-1)), ncol=k-1)Si calculamos la matriz \((\mathbf{1},\mathbf{C})^{-1}\), veremos claramente qué es cada contraste.
solve(cbind(rep(1,4),contr.simple(4)))## 1 2 3 4
## 0.25 0.25 0.25 0.25
## 2 -1.00 1.00 0.00 0.00
## 3 -1.00 0.00 1.00 0.00
## 4 -1.00 0.00 0.00 1.00Ejemplo Ahora lo podemos aplicar a nuestro ejemplo:
contrasts(coagulation$diet) <- contr.simple(4)
summary(lm(coag ~ diet, coagulation))##
## Call:
## lm(formula = coag ~ diet, data = coagulation)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.00 -1.25 0.00 1.25 5.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.400e+01 4.979e-01 128.537 < 2e-16 ***
## diet2 5.000e+00 1.528e+00 3.273 0.003803 **
## diet3 7.000e+00 1.528e+00 4.583 0.000181 ***
## diet4 2.719e-15 1.449e+00 0.000 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.366 on 20 degrees of freedom
## Multiple R-squared: 0.6706, Adjusted R-squared: 0.6212
## F-statistic: 13.57 on 3 and 20 DF, p-value: 4.658e-05Código de Helmert
Otra posibilidad es la matriz de contraste de Helmert que SAS llama Reverse Helmert coding y es la elección por defecto en S-Plus.
contr.helmert(4)## [,1] [,2] [,3]
## 1 -1 -1 -1
## 2 1 -1 -1
## 3 0 2 -1
## 4 0 0 3Se trata de comparar cada nivel con la media de los anteriores.
Ejemplo En la primera comparación se observa la diferencia entre el primer y el segundo nivel. En la segunda comparación es el tercer nivel que se compara con la media de los dos anteriores. Finalmente el cuarto nivel se compara con la media de los tres anteriores.
contrasts(coagulation$diet) <- contr.helmert(4)
summary(lm(coag ~ diet, coagulation))##
## Call:
## lm(formula = coag ~ diet, data = coagulation)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.00 -1.25 0.00 1.25 5.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.0000 0.4979 128.537 < 2e-16 ***
## diet1 2.5000 0.7638 3.273 0.003803 **
## diet2 1.5000 0.4105 3.654 0.001577 **
## diet3 -1.0000 0.2578 -3.880 0.000932 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.366 on 20 degrees of freedom
## Multiple R-squared: 0.6706, Adjusted R-squared: 0.6212
## F-statistic: 13.57 on 3 and 20 DF, p-value: 4.658e-05En el ejemplo anterior realmente la matriz de contraste debería ser
(contr.rev.helmert <- matrix(c(-1/2, 1/2, 0, 0, -1/3, -1/3, 2/3, 0,
-1/4, -1/4,-1/4, 3/4), ncol = 3))## [,1] [,2] [,3]
## [1,] -0.5 -0.3333333 -0.25
## [2,] 0.5 -0.3333333 -0.25
## [3,] 0.0 0.6666667 -0.25
## [4,] 0.0 0.0000000 0.75que es equivalente a la matriz de contraste de R salvo por una constante para cada columna.
xr <- cbind(rep(1,4),contr.rev.helmert)
solve(xr)## [,1] [,2] [,3] [,4]
## [1,] 0.2500000 0.2500000 0.2500000 0.25
## [2,] -1.0000000 1.0000000 0.0000000 0.00
## [3,] -0.5000000 -0.5000000 1.0000000 0.00
## [4,] -0.3333333 -0.3333333 -0.3333333 1.00Este código no es fácil de interpretar excepto en algunos casos especiales. En realidad, no tiene mucho sentido en el ejemplo anterior.
Si hay el mismo número de observaciones en cada nivel (un diseño balanceado), entonces las variables codificadas serán ortogonales entre sí y con la columna que representa la media general.
Otras codificaciones
Es evidente que hay muchas otras codificaciones interesantes que nos pueden ayudar en ciertos experimentos. R dispone de un código con polinomios ortogonales que se utiliza para factores ordinales:
contr.poly(4)## .L .Q .C
## [1,] -0.6708204 0.5 -0.2236068
## [2,] -0.2236068 -0.5 0.6708204
## [3,] 0.2236068 -0.5 -0.6708204
## [4,] 0.6708204 0.5 0.2236068Entre los recursos del Institute for Digital Research and Education (IDRE) de la Universidad de California (UCLA) hay una página dedicada a los sistemas de codificación de variables categóricas con R:
http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm
También en la ayuda de SAS podemos ver algunas de las matrices de contraste más útiles:
http://www.ats.ucla.edu/stat/sas/webbooks/reg/chapter5/sasreg5.htm
La elección de un código no afecta al cálculo de \(R^2\), \(\hat\sigma^2\) o al contraste \(F\) global de una regresión. Sólo modifica la estimación de los parámetros coeficientes de las variables predictoras, luego debemos conocer cual es el código si queremos interpretar dichos coeficientes.
En algunas ocasiones incluso nos puede interesar la utilización de un esquema especial. Veamos un ejemplo.
Ejemplo Baumann and Jones (ver la sección 4.4 del libro de Fox, 2002) realizaron un experimento en el que cada uno de 66 niños se asignaron al azar a uno de tres grupos experimentales, 22 a cada grupo. Los tres grupos representan diferentes métodos de enseñar a leer: un método estándar (llamado “Basal”) y dos nuevos métodos (llamados “DRTA” y “Strat”). Los investigadores realizaron dos pre-tests y tres post-tests de compresión lectora. Aquí nos centraremos en el tercer post-test. Los datos del estudio están en la base de datos del paquete .
require(car)
data(Baumann)
str(Baumann)## 'data.frame': 66 obs. of 6 variables:
## $ group : Factor w/ 3 levels "Basal","DRTA",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ pretest.1 : int 4 6 9 12 16 15 14 12 12 8 ...
## $ pretest.2 : int 3 5 4 6 5 13 8 7 3 8 ...
## $ post.test.1: int 5 9 5 8 10 9 12 5 8 7 ...
## $ post.test.2: int 4 5 3 5 9 8 5 5 7 7 ...
## $ post.test.3: int 41 41 43 46 46 45 45 32 33 39 ...attach(Baumann)Como es lógico, los investigadores estaban interesados en ver si los nuevos métodos obtenían mejores resultados que el método estándar y si los nuevos métodos eran distintos entre sí.
Un poco de estadística descriptiva:
tapply(post.test.3, group, mean)## Basal DRTA Strat
## 41.04545 46.72727 44.27273tapply(post.test.3, group, sd)## Basal DRTA Strat
## 5.635578 7.388420 5.766750En la figura se observa que puede haber diferencias entre los grupos, especialmente entre los nuevos métodos y el estándar.
En este caso es natural definir dos contrastes: (1) Basal versus la media de DRTA y Strat y (2) DRTA versus Strat:
(contrasts(group) <- matrix(c(1,-0.5,-0.5,0,1,-1), 3, 2))## [,1] [,2]
## [1,] 1.0 0
## [2,] -0.5 1
## [3,] -0.5 -1summary(lm(post.test.3 ~ group))##
## Call:
## lm(formula = post.test.3 ~ group)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.727 -3.614 1.114 3.954 12.954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.0152 0.7772 56.632 < 2e-16 ***
## group1 -2.9697 1.0991 -2.702 0.00885 **
## group2 1.2273 0.9519 1.289 0.20201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.314 on 63 degrees of freedom
## Multiple R-squared: 0.1245, Adjusted R-squared: 0.09675
## F-statistic: 4.481 on 2 and 63 DF, p-value: 0.01515Los estadísticos \(t\) para group1 y group2 contrastan justamente las dos hipótesis planteadas. Queda claro que los nuevos métodos son mejores que el anterior, pero parece que no hay diferencias significativas entre los dos nuevos métodos.
Aunque hay otras formas de hacer estos dos contrastes, este ejemplo muestra que la codificación puede ser de gran ayuda.
Bibliografia
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.
D.C Montgomery, Design and Analysis of Experiments, John Wiley & Sons, 2008.
