]]>
1Departamento de Matemáticas
Resumen
Objetivo. Presentar una metodología basada en el concepto de productos Kronecker que facilita la construcción de la matriz de varianzas y covarianzas para diseños con estructura balanceada de datos a 2 y 3-vías y una aplicación realizada en R para facilitar su cálculo y aplicación en diferentes áreas. Materiales y métodos. Se proporciona un punto de partida para personas interesadas en utilizar R en el análisis de varianza. Resultados. Se utiliza una aplicación realizada en R donde se desarrolla la metodología basada en productos Kronecker mediante la cual se construye la matriz de varianzas y covarianzas cuando se trabaja en diseños con estructura balanceada de datos desarrollada por Moya 2003. De igual forma se presenta una aplicación del método con datos reales. Conclusiones. La metodología expuesta permite agilizar el desarrollo y solución de algunos problemas prácticos. El método propuesto puede ser aplicado a modelos mixtos con efectos fijos o aleatorios con cualquier número de factores.
Palabras clave: productos Kronecker, matriz de varianzas y covarianzas, diseños balanceados, modelos lineales, R Gui.
Objective. To present a methodology based on the concept of Kronecker products that facilitates the construction of the variance and covariance matrix for designs with balanced data structure for 2 and 3 ways, and an application in R to facilitate its calculation and application in different areas. Materials and methods. We provide a starting point for people interested in using R in the analysis of variance. Results. We use an application made in R for a methodology based on Kronecker products through which we build the covariance matrix for working with designs with balanced data structure developed by Moya (2003). We also present an application of the method with real data. Conclusions. With this methodology we can accelerate the development and solution of some practical problems. The proposed methodology can be applied to mixed models with fixed or random effects with any number of factors.
Key words: Kronecker products, variance and covariance matrix, balanced designs, linear models, R Gui.
Resumo
Objetivo. Apresentar uma metodologia baseada no conceito de produtos de Kronecker para facilitar a construção da matriz de variância-covariância para desenhos com estrutura balanceada de dados em 2 e 3 vias e uma aplicação realizada em R para o cálculo fácil e aplicação em diferentes áreas. Materiais e métodos. É fornecido um ponto de partida para as pessoas interessadas em utilizar R na análise de variância. Resultados. Usa-se uma aplicação realizada em R, onde é desenvolvida a metodologia baseada nos produtos de Kronecker com a qual se constrói a matriz de variância-covariância quando se trabalha em desenhos com estrutura balanceada de dados segundo Moya (2003). Da mesma forma, se apresenta uma aplicação do método com dados reais. Conclusões. A metodologia descrita pode acelerar o desenvolvimento e solução de alguns problemas práticos. O método proposto pode ser aplicado a modelos mistos com efeitos fixos ou aleatórios com qualquer número de fatores.
Palavras-chave: produtos de Kronecker, matriz de variância-covariância, desenhos balanceados, modelos lineares, R Gui.
Introducción
Hoy en día el avance de los procesadores y el desarrollo de software estadístico han permitido que los investigadores en diversas áreas del conocimiento tengan a su disposición una gran cantidad de herramientas computacionales que les permiten un análisis óptimo de la información. Dentro de estas herramientas, se encuentran aplicaciones de distribución libre, como es el caso del R, el software estadístico más empleado a nivel mundial en investigación. Lamentablemente, en diferentes áreas aplicadas, su uso se ha visto limitado por el desconocimiento de sus procedimientos y su lenta divulgación.
El objetivo de este artículo es proporcionar un punto de partida para personas interesadas en utilizar R en el análisis de varianza, y presentar una aplicación realizada en el mismo software, donde se desarrolla la metodología basada en productos Kronecker mediante la cual se construye la matriz de varianzas y covarianzas cuando se trabaja en diseños con estructura balanceada de datos a 2 y 3-vías desarrollada por Moya (1).
Los algoritmos presentados pueden ser empleados en modelos a 1 y 2-vías con cualquier número de factores. La matriz de varianzas y covarianzas se construye bajo los supuestos usuales para modelos finitos utilizando el método propuesto por Moya (1). Se ha escogido utilizar R, debido a la universalidad de su uso y libre distribución. R se ofrece gratuitamente bajo los términos de la GNU General Public Licence; su desarrollo y distribución son llevados a cabo por varios estadísticos conocidos como el Grupo Nuclear de Desarrollo de R. Se ha tratado de simplificar las explicaciones al máximo con el fin de hacerlas lo más comprensivas posibles para que usuarios de cualquier nivel puedan utilizar los desarrollos aquí planteados.
]]> Para profundizar en la teoría subyacente al análisis de la varianza y sus aplicaciones se puede consultar, entre otros, Carmona (2), McCulloch (3), Sudhir y Saha (4) o los clásicos Hartley y Rao (5) y Searle (6). Si se quiere profundizar en el aprendizaje de R o desarrollar otro tipo de aplicaciones, en el link http://www.sciviews.org/_rgui/ se encuentra toda la información necesaria. También se pueden consultar Crawley (7) y Dalgaard (8), donde se presentan algunas aplicaciones de R.Materiales y métodos
Con el fin de realizar una adecuada descripción de los métodos de estimación de componentes de varianza, se presenta la forma más general de un modelo lineal denominado modelo mixto general. La formulación del modelo mixto general originalmente propuesto por Hartley y Rao (5) y aceptada como estándar en el campo de la estimación de las componentes de varianza esta dada por:
Donde:
yNx1: Vector de observaciones
XNxK: Matriz asociada a los efectos fijos
βKx1: Vector de parámetros asociados a los efectos fijos
Uqx1: Vector asociado a los efectos aleatorios.
ZNxq: Usualmente es una matriz de incidencia asociada a efectos aleatorios
]]> eNx1: Vector de términos del error.Como:
Donde: Zi son los vectores columna de Z.
por lo tanto:
De esta forma, la varianza de y puede ser definida como:
A e se le atribuye la estructura usual de varianza-covarianza de un componente de error, dado por:
E(e) = 0 Var(e) = ae IN [5]
]]> En primer lugar se presentan los algoritmos realizados en R para realizar análisis de varianza a 1 y 2-vías.Diseño Aleatorio a 1-vía
Con el fin de ilustrar la utilización de R en el análisis de varianza, se puede definir el modelo a 1-vía como:
Donde: τ1 denota el factor 1 con τ1 ~ N 
y el termino e corresponde al error definido anteriormente
con 
En primer lugar se deben ingresar los datos (note que existen k niveles), para tal fin se usan las siguientes instrucciones:
Sin embargo debido a la estructura de R, existe una forma adecuada para organizar la información que se pretende analizar. Primero que todo se organizan los datos en un único vector y se agrega una variable cualitativa o factor que nos indique la población correspondiente a cada caso, así:
Es necesario notar que en R, es imprescindible definir el vector tratamiento como un factor, ya que en caso contrario se podría confundir con un vector numérico y esto afectaría el análisis, además en este caso n1= n2=...= nk. Existe una instrucción que implementa los dos pasos anteriores:
Una vez se tienen los datos organizados en R, se pueden obtener algunas estadísticas descriptivas o un gráfico de medias. De igual forma, resulta útil para análisis de esta naturaleza, realizar un boxplot múltiple, mediante las instrucciones:
Si se asume que la variable Data sigue una distribución normal con igual varianza en las k poblaciones, el ANOVA o la tabla de análisis de varianza, se obtiene con:
Si se quiere ser más especifico o más formal, se puede utilizar la definición del modelo lineal y obtener los mismos resultados mediante
Además, las estimaciones de los respectivos parámetros, pueden ser calculadas con la instrucción:
Un resumen de todo el modelo con la información pertinente se obtiene con
Diseño Aleatorio a 2-vías
En primer lugar, es necesario definir el modelo. Así, El modelo a 2-vías puede ser expresado como:
Donde: τ1 denota el factor 1 con τ1 ~ N , τ2
denota el factor 2 con τ2 ~N 
, y el termino ei1, i2 corresponde al error definido anteriormente con .
En este caso los datos son introducidos mediante:
Usando comandos similares a los presentados en el caso a 1-vía, se puede generar un resumen de los datos y los respectivos gráficos que acompañan este análisis utilizando las instrucciones:
Y el modelo lineal y la tabla del análisis de la varianza se construye con:
Como se puede observar en los anteriores procedimientos, se utilizan instrucciones similares a las presentadas en el caso de 1-vía, es decir al lector interesado en realizar análisis a mayor nivel, le resultará sencillo establecer un patrón de generalización de las instrucciones de acuerdo al número de componentes que requiera.
A continuación se presenta una aplicación en R que permite la construcción de la matriz de varianzas y covarianzas cuando se trabaja en diseños con estructura balanceada de datos a 1 y 2-vías basada en la metodología propuesta por Moya (2003). En este tipo de modelo es importante determinar los mejores estimadores lineales insesgados (MELI), y los mejores predictores lineales (BLUP). Estos estimadores juegan un importante papel a la hora de realizar inferencia puntual y por intervalos sobre los parámetros. Tanto los MELIS como los BLUPS dependen del conocimiento de la estructura apropiada de la matriz de varianzas y covarianzas, así como de su inversa y es allí donde radica la importancia de la aplicación aquí presentada.
Resultados
Construcción de la matriz de varianzas y covarianzas
Un diseño factorial es aquel en el que se investigan todas las posibles combinaciones de los niveles de los factores en cada ensayo completo o réplica del experimento. Por ejemplo, si existen n1 niveles del factor τ1 y n2 niveles del factor τ2, entonces cada replica del experimento contiene todas las n1n2 combinaciones de los tratamientos. A menudo, se dice que los factores están cruzados cuando estos se arreglan en un diseño factorial. En este artículo se determina la matriz de varianza y covarianzas (V) para diseños a dos y tres vías sin tener en cuenta el efecto de réplica. Para el cálculo de la matriz V en el caso de dos vías con efecto de réplica ver Moya (1).
Modelo Aleatorio a 2-vías
A continuación de presenta la metodología mediante la cual, se construye la matriz V de manera sencilla y directa. De igual forma, se propone una metodología para hallar la inversa de esta matriz.
En la determinación de la matriz de estructura
] de matrices Is (matriz identidad de tamaño SxS) y Js (matriz de unos con tamaño SxS) como se presenta a continuación, y σi2 es la varianza de yi1,i2. Sea P el número de subíndices asociados a la variable respuesta (en este caso, P = 2), los cuales dependen de los efectos e interacciones fijas y/o aleatorias del modelo, entonces se tienen 2p particiones de productos Kronecker con matrices Is y Js , las cuales son resumidas por la expresión:
y definiendo la función indicadora:
Se puede reescribir 
, para el caso P = 2 como:
Caracterizando el vector 
, se construye el vector asociado con las componentes de varianza 
,
con esto
]]>
Donde: es la matriz final de permutaciones y θ0 representa una permutación de θp dada por
Donde: 
son los respectivos componentes de varianza.
Por otro lado, definiendo
se tiene que
Donde:
A continuación se construye una función en R que realiza este cálculo de manera sencilla con el fin de facilitar su aplicación práctica.
El primer paso en la construcción de la matriz V , es la definición del vector asociado a las componentes de varianza θ0, como:
Donde: son los respectivos componentes de varianza.
Este vector será el insumo del algoritmo realizado en R, al igual que n1 y n2 que corresponden a los niveles del primer y segundo factor aleatorio (τ1 y τ2 ), respectivamente. Para calcular la matriz V, primero que todo, se debe compilar la función (V_2vias) que aparece a continuación:
Luego para obtener la matriz V, basta con utilizar la instrucción:
Donde: Se corresponde a 
, n1 y n2 corresponden al total de niveles de τ1 y τ2, respectivamente.
Finalmente para calcular la inversa de la matriz V, primero que todo, se debe compilar la función (Vinv_2vias) que aparece a continuación:
> Vinv_2vias <- function(Se,S1,S2,n1,n2) + {
+ Theta <- matrix(c( Se, S1, S2, 0) ,nrow = 1, ncol = 4)
+ Tp <- matrix(c(1,1,1,1,0,n2,0,n2,0,0,n1,n1,0,0,0,(n1 *n2)),nrow = 4, ncol = 4)
+ Lambda <- Theta %*% t(Tp)
]]> + Lambdai <- 1/ Lambda+ Tpi <- solve(Tp)
+ M<- Tpi %*% t(Lambdai)
+ D1 <- diag(n1*n2)*M[1,1]
+ D2 <- (kronecker ((diag(n1)),matrix( 1, ncol=n2, nrow=n2)))* M[2,1]
+ D3 <- (kronecker(matrix( 1, ncol=n1, nrow=n1),(diag(n2))))* M[3,1]
+ D4 <- (matrix( 1, ncol=(n1*n2), nrow=(n1*n2)))* M[4,1]
+ D1+D2+D3+D4
+ }
Donde: Se corresponde a s2e, S1 a s2τ1, S2 a s2τ2, n1 y n2 corresponden al total de niveles de τ1 y τ2 respectivamente.
Modelo 3-vías
Utilizando una metodología similar a la aplicada en el modelo a 2-vías, V puede ser calculada como
Donde:
Donde: s2e, s2τ1, s2τ2 y s2τ3 son los componentes de varianza, y
Donde: s2e, s2τ1, s2τ2 y s2τ3 son los componentes de
]]> Además
Donde: λpes definido como en [16],
y
A continuación se construye la respectiva función en R mediante la cual se realizan los cálculos anteriormente expuestos.
Como en el caso a 2-vías, el primer paso en la construcción de la matriz de V es la definición del vector asociado a las componentes de varianza θ0, así:
Este vector será utilizado por el algoritmo realizado en R, al igual que n1, n2 y n3, que corresponden a los niveles del primer, segundo y tercer factor aleatorio (τ1, τ2 y τ3), respectivamente. Para calcular la matriz se varianzas y covarianzas, primero que todo, se debe compilar la función (V_3vias) que aparece a continuación:
> V_3vias <- function(Se,Sl,S2,S3,nl,n2,n3){ ]]>
+ (Se*diag(nl*n2*n3)) +S3* kronecker(matrix( 1, ncol=nl*n2, nrow=nl*n2),(diag(n3)))+Para obtener la matriz V, se utiliza la instrucción:
Donde: Se corresponde a σ2e, S1 a σ2τ1, S2 a σ2τ2, S3 a σ2τ3, n1, n2 y n3 corresponden a los niveles de τ1, τ2 y τ3, respectivamente.
Para calcular la inversa de la matriz V , primero que todo, se debe compilar la función (Vinv_3vias) que aparece a continuación:
> Vinv_3vias <- function(Se,S1,S2,S3,n1,n2,n3)
]]> + {+ Theta <- matrix(c( Se, 0,0,0,S3,S2, S1, 0) ,nrow = 1, ncol = 8)
+ Tp1 <- kronecker(matrix(c(1,1,0,n2),nrow = 2, ncol = 2), matrix(c(1,1,0,n3),nrow = 2, ncol = 2))
+ Tp<- kronecker(matrix(c(1,1,0,n1),nrow = 2, ncol = 2),Tp1)
+ Lambda <- Theta %*% t(Tp)
+ Lambdai <- 1/ Lambda
+ Tpi <- solve(Tp)
+ M<- Tpi %*% t(Lambdai)
+ D1 <- diag(n1*n2*n3)*M[1,1]
+ D2 <- (kronecker ((diag(n1*n2)),matrix( 1, ncol=n3, nrow=n3)))* M[2,1]
]]> + D3 <-kronecker(diag(n1), (kronecker(matrix( 1, ncol=n2, nrow=n2),(diag(n3)))))* M[3,1]+ D4 <- kronecker(matrix( 1, ncol=n1, nrow=n1) , (matrix( 1, ncol=(n2*n3), nrow=(n2*n3))))* M[4,1]
+ D5 <- kronecker(matrix( 1, ncol=(n1*n2), nrow=(n1*n2)) , (matrix( 1, ncol=n3, nrow=n3)))* M[5,1]
+ D6 <-kronecker(matrix( 1, ncol=n1, nrow=n1), (kronecker(diag(n2),matrix( 1, ncol=n3, nrow=n3))))* M[6,1]
+ D7 <- (kronecker ((diag(n2)),matrix( 1, ncol=(n2*n3), nrow=(n2*n3))))* M[7,1]
+ D8 <- matrix ( 1, ncol=(n1*n2*n3), nrow=(n1*n2*n3))* M[8,1]
+ D1+D2+D3+D4+D5+D6+D7+D8
+ }
Así, para determinar la matriz V-1, se utiliza la instrucción:
Donde: Se corresponde a σ2e, S1 a σ2τ1, S2 a σ2τ2, S3 a σ2τ3, n1, n2 y n3 corresponden a los niveles de τ1, τ2 y τ3, respectivamente.
Discusión
La estimación de componentes de varianza es útil en la identificación de fuentes de variación. Existen varios métodos de estimación de componentes de varianza. El problema de la estimación de los componentes de varianza se ha abordado históricamente desde dos perspectivas diferentes: estructuras balanceadas y desbalanceadas de datos. Para la primera, la estimación recae exclusivamente en el método propuesto por Henderson, mientras que para la segunda, existen varios procedimientos, entre ellos, tres extensiones distintas del método utilizado para datos balanceados. Se debe distinguir entre estas dos estructuras puesto que ambas se basan en la posibilidad que existe de encontrar formas cuadráticas no definidas positivas en una estructura desbalanceada de datos.
Después de tener un algoritmo que permite deducir la matriz de varianzas y covarianzas y la correspondiente inversa, resultan de gran utilidad las funciones realizadas en R, mediante las cuales se pueden llevar a cabo múltiples aplicaciones. La discusión entorno a este tema queda abierta y este trabajo permite explorar diferentes tipos de modelos con estructuras diferentes a los aquí expuestos.
Aplicación del modelo
González (9), evalúa el efecto del agua residual urbana con diferentes niveles de tratamiento y laminas de riego sobre las características físicas, químicas y carga bacterióloga de un suelo agrícola. Los tratamientos se establecieron a través del diseño completamente al azar con arreglo combinatorio factorial 4x3, correspondiente a 12 tratamientos. Los factores fueron calidad del agua (suelo regado con agua cruda, suelo regado con agua tratada a nivel primario, suelo regado con agua tratada a nivel secundario y suelo regado con agua blanca), y el segundo factor fue las laminas de riego 8,8 10 mm/día. El autor realizó un análisis de varianza a cada una de las variables evaluadas para determinar el efecto de los tratamientos. Para esta aplicación se tomó el análisis de varianza para el contenido de arcilla en el suelo regado con agua residual. De esta forma los componentes de varianza están dados por:
Entonces aplicando la función V_2vias (Para este caso V_2vias (18.20, 28.4 , 8.53, 4.3)), se obtiene la matriz de varianzas y covarianzas (V ) dada por:
]]>
De igual forma, aplicando la función Vinv_2vias se obtiene la inversa (V -1 ) dada por:
Conclusiones
La metodología expuesta permite agilizar el desarrollo y solución de algunos problemas prácticos mediante la estimación de la matriz de varianzas y covarianzas. El método propuesto puede ser aplicado a modelos mixtos con efectos fijos o aleatorios con cualquier número de factores y mediante la aplicación realizada en R se logran implementar los cálculos de manera óptima y eficiente. En este tipo de modelo es muy importante determinar los mejores estimadores lineales insesgados (MELI), y los mejores predictores lineales (BLUP), con el objeto de realizar inferencia puntual y por intervalos para tomar decisiones apropiadas en muchos problemas de aplicaciones. Como los MELIS y los BLUPS dependen del conocimiento de la estructura apropiada de la matriz de varianzas y covarianzas y de su inversa, el algoritmo expuesto permite optimizar el cálculo de estos estimadores.
Financiación
Este trabajo no tuvo ningún tipo de financiación externa y fue realizado dentro de las labores académicas de los autores.
Conflicto de intereses
]]> No existe conflicto de intereses.Referencias
1. Moya LM. Obtención de la matriz de varianzas y covarianzas a través de los productos Kronecker para modelos balanceados. Universitas Scientiarum 2003; 8 (2): 39-44. [ Links ]
2. Carmona F. Modelos lineales. Volumen 3 de Electronic-university mathematical books. Editor Publicacions i Edicions, Universitat de Barcelona, Barcelona, España. 2005, 266 p. [ Links ]
3. MacCulloch, C.E. Generalized mixed models. Volume 3. NSF-CMBS Regional Conference Series in Probability and Statistics. 2003, 84 p. [ Links ]
4. Sudhir P, Saha K. The generalized linear model and extensions: A review and some biological and environmental applications. Environmetrics 2007; 18:421-443. [ Links ]
5. Hartley HO, Rao JNK. Maximum likelihood estimation for the mixed analysis of variance models. Biometrika 1967; 54: 93-108. [ Links ]
6. Searle SR. Built in restrictions on best linear unbiased predictors (BLUP) of random effects in mixed models. The American Statistician 1997; 51: 19-21. [ Links ]
7. Crawley M. Biodiversidad statistics: an introduction using R. Primera edición. John Wiley & Sons. London. England. 2005, 327 p. [ Links ]
8. Dalgaard P. Introductory statistics with R. Primera edición. Springer. New York, USA. 2008, 363 p. [ Links ]
9. González J. Efecto de la adición de agua residual urbana sobre las características de un suelo Agrícola. Tesis Doctoral, área: Biotecnología. Universidad de Colima, Tecomán, México, 2003, 133 p. [ Links ] ]]>