i UNIVERSIDAD DE COSTA RICA SISTEMA DE ESTUDIOS DE POSTGRADO EL DILEMA DE “SHRINKAGE” O FENÓMENO DE ENCOGIMIENTO EN MODELOS MIXTOS, UNA COMPARACIÓN ENTRE LOS ENFOQUES BAYESIANO Y FRECUENTISTA Tesis sometida a la consideración de la Comisión del Programa de Estudios de Posgrado en Estadística para optar al grado y título de Maestría Académica en Estadística ANDREA DE LOS ÁNGELES CÉSPEDES SOLÍS Cuidad Universitaria Rodrigo Facio, Costa Rica 2021 ii Dedicatoria y agradecimientos Este trabajo se lo dedico a mi familia por ser mi apoyo incondicional en esta travesía y aventura universitaria, especialmente a mi mamá por ser mi pilar, confiar en mí, apoyarme, amarme en todo momento y ser parte de mis decisiones, y a quien es el amor de mi corazón, a mi tía Carmen, quien fue una segunda mamá y empezó acompañarme desde el cielo hace unos meses, espero estén orgullosas; a mi esposo que ha estado siempre ahí y es el amor de mi vida, a mi pequeña bendición, Horacio esto toma más fuerza por ti. A todos los colegas que estuvieron anuentes a colaborarme ante las dudas que me surgieron en el camino, que aportaron a mi formación profesional en este proceso de aprendizaje, a mi tutora por darme un reto académico y ser parte de esto, que además de profesionalmente, me permitió crecer personalmente. Finalmente, a los todos los que de una u otra manera formaron parte de esto y me alentaron hasta su conclusión; son parte importante de mi corazón. iii “Esta tesis fue aceptada por la Comisión del Programa de Estudios de Posgrado de la Universidad de Costa Rica, como requisito parcial para optar al grado y título de Maestría Académica en Estadística” _______________________________ Dr. Ricardo Alvarado Barrantes Representante del Decano Sistema de Estudios de Postgrado _______________________________ Dra. Eiliana Montero Rojas Directora de Tesis _______________________________ Dr. Guaner Rojas Rojas Lector _______________________________ Dr. Juan Carlos Brenes Sáenz Lector _______________________________ M.Sc. Johnny Madrigal Pana Director del Programa de Posgrado en Estadística _______________________________ Andrea de los Ángeles Céspedes Solís Sustentante iv Tabla de contenidos Dedicatoria................................................................................................................ii Agradecimientos........................................................................................................ii Hoja de aprobación...................................................................................................iii Tabla de contenidos..................................................................................................iv Resumen..................................................................................................................vi Abstract...................................................................................................................vii Lista de cuadros…..................................................................................................viii Lista de gráficos........................................................................................................ix Lista de imágenes……………………………………………………….....…………......xi 1. Introducción....................................................................................................1 2. Justificación y objetivos..................................................................................7 2.1 Objetivo general........................................................................................9 2.2 Objetivos específicos................................................................................9 3. Marco teórico..................................................................................................9 4. Marco metodológico.....................................................................................10 4.1 Población bajo estudio............................................................................10 4.2 Análisis de datos originales.....................................................................10 4.3 Análisis de datos en esta investigación...................................................27 4.4 Técnicas a emplear.................................................................................29 4.5 Modelos analizados….............................................................................34 4.6 Alternativas de estimación......................................................................40 4.6.1 Regresión LASSO..................................................................40 4.6.2 Laplace...................................................................................43 4.6.3 Matriz de dispersión................................................................45 4.7 Estudio de simulaciones..........................................................................48 4.7.1 Diseño del estudio de simulaciones........................................48 4.7.2 Generación y análisis de los datos..........................................50 5. Análisis de datos y resultados......................................................................53 v 5.1 Análisis descriptivo.............................................................................53 5.2 Modelos generados................................................................................57 5.3 Gráficas de coeficientes estimados........................................................72 5.4 Del estudio de simulaciones...................................................................79 6.Conclusiones………………….……….............................................................91 6.1 Conclusiones relevantes del caso de estudio..........................................91 6.2 Hallazgos adicionales.............................................................................96 Bibliografía..............................................................................................................99 Anexos..................................................................................................................106 Anexo 1: Estadísticas descriptivas.................................................................106 Anexo 2: Selección del modelo.......................................................................109 Anexo 3: Paquetes de R utilizados…………...................................................111 3.1 glmer….................................................................................................111 3.2 MCMCGLMM........................................................................................112 3.3 glmmLasso...........................................................................................115 3.4 glmmsr..................................................................................................116 3.5 hglm2....................................................................................................117 3.6 glmmTBM.............................................................................................119 Anexo 4: Estimación adicional de modelos a partir del enfoque frecuentista y bayesiano.......................................................................................................121 4.1 Estimación frecuentista con R...............................................................121 4.2 Estimación frecuentista y bayesiana con STATA 15............................122 Anexo 5: Selección del valor para Lambda en la regresión LASSO................130 5.1 Uso de Validación cruzada....................................................................130 5.2 Uso de BIC y AIC..................................................................................131 5.3 Resumen de estimaciones....................................................................131 Anexo 6: Estimación de coeficientes a partir del estudio de simulaciones.......134 6.1 Estimación del SESGO.........................................................................134 6.2 Estimación del CME..............................................................................137 vi Resumen En los modelos mixtos (multinivel) la especificación matemática del modelo puede impactar de manera relevante en los resultados sustantivos obtenidos. El fenómeno de shrinkage o fenómeno de encogimiento provoca que los efectos fijos comiencen a reducirse al punto de no presentar significancia estadística, diversos autores señalan como los principales causantes de este dilema analítico al aumento en la cantidad de parámetros, la cantidad de elementos en los niveles (n), la varianza entre y dentro de los conglomerados, y el enfoque de estimación. En este estudio se presentaron dos propuestas de estimadores shrinkage, una correspondiente a la ponderación de los parámetros del modelo y otra ajustando las matrices de varianza y covarianza. A partir de estas dos se utilizaron los siguientes métodos de estimación LASSO, Laplace y ajuste por matriz de dispersión, además de las estimaciones frecuentista y bayesiana consideradas preliminarmente. Los datos utilizados correspondieron a experimentos realizados en el Centro de Investigación en Neurociencias de la Universidad de Costa Rica. Por la estructura de los datos, se utilizó un modelo lineal generalizado multinivel con distribución Poisson de enlace logaritmo, y de tres niveles. Así mismo, se realizó un estudio de simulaciones considerando un modelo base y distintos escenarios a partir de los posibles causantes del efecto shrinkage, con el fin de comparar las técnicas de estimación y determinar cuál es la mejor. El análisis indica que la estimación bayesiana es la que tiene mayor estabilidad en las estimaciones de los coeficientes, también se tiene que el tamaño de muestra no está relacionado con el efecto shrinkage, pero, la varianza, y la cantidad de coeficientes aleatorios sí. En el caso de las estimaciones bayesianas se utilizaron 30 000 iteraciones con un burn-in de 10 000. Para el análisis, generación de bases de datos y estudio de simulación, se utiliza la plataforma de programación R 3.5.2 y el paquete STATA 15. vii Abstract In mixed (multilevel) models, the mathematical specification of the model can have an important impact on the substantive results obtained. The phenomenon of shrinkage causes the fixed effects to begin to reduce to the point of not presenting statistical significance, various authors point out as the main causes of this analytical dilemma to the increase in the number of parameters, the elements in the levels (n), the variance between and within the clusters, and the estimation approach. In this study, they present two proposals for shrinkage estimators, one corresponding to the weighting of the model parameters and the other adjusting the variance and covariance matrices. From these two, the following LASSO, Laplace, and dispersion matrix adjustment methods are used, in addition to the frequenter and Bayesian estimates previously considered. The data used correspond to experiments carried out at the Neuroscience Research Center of the University of Costa Rica. Due to the structure of the data, a multilevel generalized linear model with a Poisson distribution of logarithmic link, and three levels, was used. Likewise, a simulation study is carried out considering a base model and different scenarios from the possible causes of the shrinkage effect, to compare the estimation techniques and determine which the best is. The analysis indicates that the Bayesian estimate is the one with the greatest stability in the estimations of the coefficients; it is also found that the sample size is not related to the shrinkage effect, but the variance, and the number of random coefficients, is. In the case of Bayesian estimates, 30,000 iterations with a burn-in of 10,000 are used. For analysis, database generation and simulation study, the R 3.5.2 programming platform, and the STATA 15 package are used. viii Lista de cuadros Cuadro 1: Descripción de los escenarios aplicados en las simulaciones……......…50 Cuadro 2: Estadísticos generales de las vocalizaciones emitidas según día del experimento................................................................…………………………….…57 Cuadro 3: Estimación frecuentista para los coeficientes de los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos………………...……...…..…61 Cuadro 4: Estimación bayesiana para los coeficientes de los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos………………….………...…65 Cuadro 5: Estadísticos de resumen para los modelos estimados con R según enfoque frecuentista y bayesiano……………………………………….…….……..…66 Cuadro 6: Estimación con regresión LASSO para los coeficientes en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos…...…………..…….68 Cuadro 7: Estimación con cadenas de Markov con aproximación Laplace para los coeficientes de los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos………..…………………….…..……………..……………………….……...70 Cuadro 8: Estimación con matriz dispersa para los coeficientes de los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos…………………….71 Cuadro A1: Estadísticas generales de las vocalizaciones emitidas según sujetos del grupo control…………………………………………..……………………...……..….107 Cuadro A2: Estadísticos generales de las vocalizaciones emitidas según sujetos del grupo tratamiento…………………………………………………..………..........…...108 Cuadro A3: Ranking de modelos seleccionados preliminarmente para el análisis……………………………………………………………………………..…....110 ix Cuadro A4: Estimación frecuentista de coeficientes para los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, utilizando el paquete glmmTMB……………………………………..…………………………………..……..121 Cuadro A5: Estimación de coeficientes fijos para los modelos de efectos fijos, según valores de Lambda…….....................................................................................…131 Cuadro A6: Estimación de los errores estándar para los coeficientes fijos en los modelos de efectos fijos, según valores de Lambda………………………….…….132 Cuadro A7: Estimación de coeficientes fijos para los modelos de efectos aleatorios completos, según valores de Lambda…...…………...………………...…..….….…132 Cuadro A8: Estimación de los errores estándar para los coeficientes fijos en los modelos de efectos aleatorios completos, según valores de Lambda….……..…133 Lista de gráficos Gráfico 1: Vocalizaciones emitidas por ratas según grupo………..…..…….…….....53 Gráfico 2: Vocalizaciones emitidas por día del experimento…………..……….…....54 Gráfico 3: Vocalizaciones emitidas durante el experimento por ratas del grupo control…..……………………………………………………………………..…….…....55 Gráfico 4: Vocalizaciones emitidas durante el experimento por ratas del grupo tratamiento…..…………………………………………………………………..….…....56 Gráfico 5: Coeficientes del intercepto en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación….………….…73 Gráfico 6: Coeficientes de DIA2 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación……………..….……..…..73 x Gráfico 7: Coeficientes de DIA3 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..……….…….………...74 Gráfico 8: Coeficientes de M1 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..………...…………….....75 Gráfico 9: Coeficientes de M4 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..………………….……..75 Gráfico 10: Coeficientes de M5 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..………….…...……...….76 Gráfico 11: Coeficientes de M6 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..……………………..….76 Gráfico 12: Coeficientes de M9 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..……….….......…..….…77 Gráfico 13: Coeficientes de M10 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..………....……..…….78 Gráfico 14: Coeficientes de M11 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación…..………....…….……..78 Gráfico A1: Cantidad de vocalizaciones emitidas por los sujetos del grupo control según día………………………………………………………………………….…….106 Gráfico A2: Conteo de vocalizaciones emitidas por las ratas del grupo control durante el experimento…………………………………………………….………..…106 Gráfico A3: Cantidad de vocalizaciones emitidas por los sujetos del grupo tratamiento según día…………………………………………………………….....…107 Gráfico A4: Conteo de vocalizaciones emitidas por las ratas del grupo tratamiento durante el experimento………………………………………..……………...….……108 xi Lista de imágenes Imagen 1: Diagrama de la estructura de los datos…………………………………....28 Imagen 2: Diagrama de la estructura de los datos analizada…..…..………………..37 Imagen 3: SESGO del coeficiente DIA3 en los modelos de coeficientes fijos y efectos aleatorios, según tipo de estimación, tamaño de muestra y varianza…...…81 Imagen 4: SESGO del coeficiente M1 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza……………………………………………………………………….82 Imagen 5: SESGO del coeficiente M5 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…………………………………………………………………….…83 Imagen 6: SESGO del coeficiente M10 en los modelos de coeficientes fijos y efectos aleatorios, según tipo de estimación, tamaño de muestra y varianza…………….…84 Imagen A1: Estimación frecuentista de coeficientes para el modelo de efectos fijos utilizando un modelo lineal generalizado mixto…………………….…………….….122 Imagen A2: Estimación frecuentista de coeficientes para el modelo de efectos aleatorios utilizando un modelo lineal generalizado mixto…………...………....….123 Imagen A3: Estimación frecuentista de coeficientes para el modelo de efectos fijos utilizando un modelo multinivel de regresión Poisson...……..……………...….…..124 Imagen A4: Estimación frecuentista de coeficientes para el modelo de efectos aleatorios utilizando un modelo multinivel de regresión Poisson.………..…....…..125 Imagen A5: Estimación bayesiana de coeficientes para el modelo de efectos fijos utilizando un modelo multinivel con regresión Poisson.………..…...................….126 xii Imagen A6: Estimación bayesiana de coeficientes para el modelo de efectos aleatorios utilizando un modelo multinivel con regresión Poisson……………......128 Imagen A7: SESGO del coeficiente intercepto en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza……………………………………………….………………....134 Imagen A8: SESGO del coeficiente DIA2 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza……….………………………………………………………….….134 Imagen A9: SESGO del coeficiente M4 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza……………………………………………………………………...135 Imagen A10: SESGO del coeficiente M6 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…….…………………………………………………………..……135 Imagen A11: SESGO del coeficiente M9 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…..…………………………………………………………....……136 Imagen A12: SESGO del coeficiente M11 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza……………….……………………………………………...…..…136 Imagen A13: CME del coeficiente intercepto en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…………………………………..……………………………...….137 xiii Imagen A14: CME del coeficiente M1 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza….……………………………………..…………………..…….…137 Imagen A15: CME del coeficiente M4 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…………………………..…………………………………....……138 Imagen A16: CME del coeficiente M5 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…………………………………………..…………………………138 Imagen A17: CME: coeficiente M6 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza………………………………………..…………………..……………...…….139 Imagen A18: CME del coeficiente M9 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza……………………………………..……………………..……..…139 Imagen A19: CME del coeficiente M10 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…….………………………………………………………….….…140 Imagen A20: CME del coeficiente M11 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…………………………………………………………………......140 Imagen A21: CME del coeficiente DIA2 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza………………………………………………………………..........141 xiv Imagen A22: CME del coeficiente DIA3 en los modelos de efectos fijos, efectos aleatorios y efectos aleatorios completos, según tipo de estimación, tamaño de muestra y varianza…………………………………………………………………..….141 xv 1 1. Introducción Los modelos de regresión son ampliamente utilizados en áreas del conocimiento como administración de empresas, economía, ingeniería, ciencias sociales, salud y biología, entre otras; este método de análisis genera para cada una de estas áreas un mayor conocimiento del fenómeno de interés y les brinda insumos para la toma de decisiones; sin embargo, la correcta aplicación de esos modelos implica una buena comprensión de la teoría que los sustenta, así como del contexto al cual se aplican (Kutner et al., 2005). Los supuestos básicos de los modelos de regresión son: linealidad, independencia de los residuos, homocedasticidad, normalidad en la distribución de los residuos y la no colinealidad. A pesar de que el investigador cuenta en la mayoría de los casos con técnicas que le permiten realizar el análisis de los datos cuando no se cumple alguno de estos supuestos, debe considerar el uso de otro método de análisis ante el no cumplimiento de independencia entre las observaciones, ya que los modelos de regresión lineal no consideran estas semejanzas o disimilitudes que se gestan en la estructura natural de los datos (Gaviria & Castro, 2005). La no independencia presente en los datos puede ser el reflejo de una compleja estructura de anidamiento; la cual es la principal particularidad de los modelos multinivel, y que los diferencia de los modelos de regresión lineal clásica (Gaviria & Castro, 2005). Los modelos multinivel de dos niveles respetan la estructura jerárquica1 de los datos ya que incorporan mecanismos de estimación que permiten realizar cálculos más precisos para los parámetros del modelo cuando se tienen distintos niveles. También permiten ventajas analíticas al incluir la interacción entre individuos de diferentes niveles, así como considerar más de un contexto y estimar posibles 1 Dentro de la jerarquía se identifican los niveles, un modelo mixto debe presentar al menos dos niveles, en donde el primer nivel corresponde a la unidad mínima estudiada y los niveles siguientes contienen la información del nivel anterior. (Aparicio & Morera, 2007) 2 diferencias entre ellos (Gaviria & Castro, 2005). Por tanto, la no comprensión de la estructura de los datos e incorrecta especificación matemática del modelo puede impactar de manera relevante en los resultados sustantivos obtenidos, por lo cual es importante que el investigador conozca estos modelos y los utilice en sus análisis con el fin de evitar estimaciones erróneas, así como conclusiones erradas. Dentro de la familia de los modelos multinivel se incluyen los modelos con datos longitudinales, en los cuales los conglomerados corresponden al sujeto de análisis y las unidades que conforman el conglomerado son las mediciones realizadas al sujeto a través de tiempo (Congdon, 2003, pp. 136-137). Por otra parte, los modelos de medidas repetidas hacen referencia a datos en los cuales el sujeto representa el conglomerado y las mediciones realizadas a este corresponden a las unidades del conglomerado, pero incluyendo la particularidad de que las condiciones de observación o del experimento puede someterse a un cambio (West et al., 2007). Algunas ramas del conocimiento en las cuales las estructuras jerárquicas predominan en los datos son, por ejemplo: la educación, en la cual los centros educativos son los conglomerados y los estudiantes son los sujetos dentro del conglomerado. En el área de la salud, por ejemplo, existen dos escenarios: el primero considera a los centros de salud como posibles conglomerados conformados por los individuos y el segundo hace referencia al mismo sujeto como conglomerado a quien se le realizan varias mediciones. En el área psicológica los modelos mixtos son utilizados en los análisis conductuales de sujetos. En agronomía se utilizan en el análisis de tratamientos en cultivos, así como la efectividad de plaguicidas o herbicidas, entre otros; sin embargo, debido a la complejidad de los problemas a resolver, y a la necesidad de contar con programas estadísticos que incorporen los algoritmos necesarios, el uso de los modelos multinivel se da siempre y cuando los programas computacionales respondan a las necesidades del analista (Catalán et al, 2003, citado por Aparicio & Morera, 2007). Además, gracias a las mejoras en los software en las últimas décadas la utilización de estos modelos ha ido en aumento, y es a finales de la década de los ochentas 3 cuando se produce una gran cantidad de publicaciones relacionadas con los modelos multinivel (Gaviria & Castro, 2005, pág. 8). A pesar de las ventajas analíticas que ofrecen modelos multinivel, también presentan algunas problemáticas, entre las que resalta la disyuntiva entre priorizar la estimación de los efectos fijos2 o la estimación de los efectos aleatorios3. Esta prioridad debe responder a los objetivos de la investigación, así como a los intereses del investigador; sin embargo, resulta importante saber que en la estimación de los efectos fijos y aleatorios del modelo se tiene un fenómeno llamado “shrinkage” o factor de encogimiento, el cual tiene un papel importante en las estimaciones de los efectos fijos. Si bien al realizar la comparación estadística entre el modelo con más parámetros versus el modelo con menos parámetros, la razón de verosimilitud arrojará resultados a favor del modelo más completo, pero debido a que la significancia estadística no implica importancia práctica, es necesario buscar otras alternativas para decidir cuál es el modelo más adecuado; también resulta importante conocer otros posibles métodos en la estimación del factor de encogimiento o de los coeficientes del modelo, con el fin de no limitar al investigador a un número determinado de coeficientes aleatorios o fijos. La problemática inicia, en términos matemáticos, con la presencia de correlación entre las variables predictoras. Algunas consecuencias de la colinealidad son predicciones poco fiables ya que la influencia de cada predictor se encuentra solapada por la relación con otros predictores, también se tendrían asociados “grandes” errores estándar y valores 𝑡 no significativos; además, los coeficientes serán sensibles a cambios en los datos (Gujarati & Porter, 2010, pág. 327.). De esta 2 Gaviria & Castro (2005) “Los parámetros fijos corresponden a los efectos medios en la población.” (pág. 59) 3 De la Cruz (2008) menciona que los efectos aleatorios se traducen en un modelo de coeficientes aleatorios que van a tomar en cuenta la variabilidad entre agrupamientos, que a su vez pueden tomar formas simples como a través de variabilidad a nivel del intercepto, o de formas más complejas, a través de variabilidades entre niveles o contextos. (pág. 3) 4 manera, al especificar los efectos aleatorios en el modelo, la estimación de los efectos fijos podría reducirse al punto de no presentar significancia estadística Gaviria & Castro (2005). Por ello, López-González (1998), Soto-Liria et al. (2000) proponen algunas alternativas para solventar este problema, entre ellas: 1- eliminar las variables problemáticas (con ello se puede cometer sesgo de especificación del modelo), 2- aumentar el tamaño de n (no siempre es posible), 3- realizar transformaciones en las variables (la interpretación de los resultados obtenidos puede ser compleja) y finalmente, 4- incorporar información externa a los datos originales, en este punto se hace referencia a un ajuste en el cálculo del coeficiente para disminuir la colinealidad de las variables mediante los llamados factores shrinkage o estimadores por contracción, que buscan reducir el error cuadrático medio para corregir la colinealidad. En el presente documento, se analizará el último punto señalado. Se parte de que los principales causantes del fenómeno shrinkage en las estimaciones son: un aumento en la cantidad de parámetros del modelo (Royle & Link, 2002), la cantidad de elementos en los niveles, así como la varianza entre y dentro de los conglomerados (Gaviria &Castro, 2005), y que esto ocurre de manera natural tanto en estimaciones de enfoque frecuentista como bayesiano. Existen métodos de estimación por contracción, así como de selección de variables y estimación por contracción que realizan de manera simultánea. Estos tienen un gran potencial para muchos problemas, especialmente aquellos que involucran pequeños tamaños de muestra y / o modelos altamente parametrizados (Royle-Link, 2002). El factor shrinkage puede ser incorporado en la estimación de los coeficientes del modelo, de las siguientes maneras: 5 López (1998) indica que mediante una matriz definida positiva como en el caso de la regresión ridge en los modelos de regresión lineal4, el estimador contraído es incorporado en el cálculo de una matriz definida positiva mediante una constante 𝑘, dicha modificación consiste en añadir pequeñas cantidades positivas a la diagonal de la matriz X'X5, los estimadores resultantes son sesgados, pero ya que su media cuadrática es menor, son más estables por lo que no se verán afectados por las variaciones de los datos. En este caso, una de las principales ventajas del uso de estimadores ridge es que con un 𝑘 lo suficientemente pequeño se logra reducir los efectos de la colinealidad y el estimador estaría más próximo al verdadero estimador mínimo cuadrático. Por otra parte, se puede incluir un ponderador 𝛌 como lo muestra Gaviria & Castro (2005, pág. 85). Este ponderador se denomina parámetro shrinkage o factor de encogimiento, su cálculo corresponde a la razón que involucra la varianza de las medias de los conglomerados, y la varianza entre los conglomerados, tal como se presenta a continuación: 𝜆𝑖 = 𝜎𝜇0 2 𝜎𝜇0 2 + 𝜎𝜀 2 𝑛𝑗 Este ponderador puede tomar valores de cero a uno siendo estos el mínimo valor posible y máximo valor posible respectivamente6. Otro ponderador utilizado en los modelos multinivel es el llamado James-Steins, el cual incorpora un componente 4 La ecuación correspondiente a la estimación del parámetro en una regresión lineal es �̂� = (𝑋´𝑋) −1 𝑋´𝑌. 5 La ecuación de regresión con la respectiva modificación, corresponde a �̂� = (𝑋´𝑋 + 𝐾𝑄) −1 𝑋´𝑌 6 En el caso de 𝜆𝑖 = 1, se tienen dos escenarios: 1. 𝜎𝜀 2 = 0 : en este caso los sujetos son iguales en el conglomerado. 2. 𝑛𝑗 = ∞ : la cantidad de sujetos dentro del conglomerado es infinita. Para 𝜆𝑖 = 0, se tienen dos escenarios: 1. 𝜎𝜇0 2 = 0 : todos los conglomerados tiene la misma media. 2a. 𝜎𝜀 2 > 𝜎𝜇0 2 : cuando las varianzas dentro de los conglomerados con mayores que las varianzas entre los conglomerados, los residuos medios de cada conglomerado no reflejan las diferencias del conglomerado respecto a la media general. 2b. 𝑛𝑗 → 0: si la cantidad de sujetos en el conglomerado es muy pequeña que tiende a cero, los residuos medios para estimar 𝜎𝜇0 2 introducen mucho error. (Gaviria & Castro, 2005) 6 denominado 𝜋, que corresponde a un factor de contracción que toma valores entre cero y uno (Ejaz Ahmed et al., 2007). En la revisión de literatura se describen algunos estudios relacionados con los estimadores insesgados para modelos multinivel, propuestas para el análisis de las matrices de covarianza, así como modelos bayesianos jerárquicos y análisis de factor de shrinkage o de encogimiento de casos aplicados; la literatura consultada y presentada en este trabajo aporta conocimiento al tema de interés desde distintas aristas. Finalmente, para las estimaciones de los modelos multinivel se utiliza la plataforma de programación R 3.5.2 y el programa estadístico STATA 15. 7 2. Justificación y objetivos El objetivo del presente trabajo es comparar los enfoques bayesiano y frecuentista al estimar modelos multinivel, cuando el investigador se enfrenta al shrinkage o fenómeno de encogimiento. Asimismo, se desea identificar posibles alternativas a la estimación de los coeficientes en los modelos multinivel cuando se tiene este dilema analítico. El fenómeno del shrinkage o encogimiento ocurre cuando en un modelo multinivel se incorpora una mayor cantidad de coeficientes aleatorios en el análisis, de modo que a mayor cantidad de coeficientes aleatorios mayor encogimiento o disminución en las estimaciones de los coeficientes fijos del modelo, de esta manera, el investigador se enfrenta ante la disyuntiva de decidir la cantidad determinada de coeficientes en el análisis o escoger variables importantes a incluir en el modelo, y de esta manera no afectar las estimaciones. El tema toma relevancia ya que en este caso específico se pretende aplicar en el área de la neurociencia, sentando las bases para estudios más sofisticados, pues en el país se ha trabajo poco este tema, y los estudios previos se han limitado al análisis de varianza (ANOVA). Por otra parte, el uso de los modelos multinivel adquiere protagonismo por su particularidad de considerar la estructura jerárquica de los datos para las estimaciones; por ejemplo, estos son muy utilizados en áreas como la educación, medicina y agronomía, entre otras; y, en el caso puntual de los datos recabados en neurociencias, la naturaleza de los mismos involucra una jerarquía que debe ser respetada. Aunado a lo anterior se incorpora en el análisis a los modelos lineales generalizados multinivel, debido a que la variable de respuesta no es normal, por lo que se requieren métodos más flexibles en términos de supuestos y que generen resultados certeros. 8 Además, se realiza un estudio de simulaciones considerando un modelo base que se especifica en el apartado de aspectos metodológicos, asociado a distintos escenarios – en los cuales la teórica indica la posible presencia de shrinkage en las estimaciones – de análisis que buscarán comparar las estimaciones frecuentistas y bayesianas obtenidas de los coeficientes, además de las alternativas propuestas y con ello determinar el mejor método de ser posible, para solventar la problemática del shrinkage en el modelo multinivel de interés. El estudio de simulaciones resulta útil en la generación de distintos escenarios, lo que permite al investigador comprender más el fenómeno de encogimiento; también, se considera relevante pues permite identificar y analizar los posibles cambios que se generan en las estimaciones de los efectos fijos al incrementar la cantidad de efectos aleatorios, tamaños de muestra, así como la varianza de los datos, y realizando la comparación entre los métodos frecuentista y bayesiano. Para responder al objetivo de análisis planteado se analizarán los datos correspondientes a varios experimentos neurocientíficos realizados en el Centro de Investigación en Neurociencias de la Universidad de Costa Rica, y aplicado a 31 ratas macho de cepa Wistar, en donde se mide el número de vocalizaciones emitidas en momentos específicos durante un periodo de 10 minutos, por tres días consecutivos. Por la estructura de los datos, se plantea utilizar un modelo de medidas repetidas de tres niveles incorporando métodos de análisis con enfoques frecuentista y bayesiano. 9 Los objetivos planteados en la presente investigación se presentan a continuación. 2.1 Objetivo General Comparar los enfoques bayesiano y frecuentista para determinar el modelo multinivel más adecuado enfrentando el problema de shrinkage o fenómeno de encogimiento. 2.2 Objetivos específicos 1. Presentar estrategias estadísticas y no estadísticas para enfrentar el problema de shrinkage o fenómeno de encogimiento en los enfoques frecuentista y bayesiano en los distintos modelos, y seleccionar aquellas que se consideren más pertinentes, más allá de la discusión en torno a la significancia estadística. 2. Aplicar los métodos de estimación LASSO, Laplace y matriz de dispersión, para enfrentar el problema de shrinkage en modelos mixtos de datos reales provenientes de los experimentos 723-B7-610 “Motivación y plasticidad neuroconductual: Efectos cognitivos, emocionales y sociales del enriquecimiento ambiental aleatorio e impredecible en ratas” y 742-B6-291 "Cambios en la plasticidad neuronal relacionados con la incubación de la sensibilidad por la anfetamina: evaluación molecular y estructural". 3. Comparar los métodos de estimación propuestos, así como la frecuentista y bayesiana para enfrentar el shrinkage a partir de escenarios simulados, e identificar el mejor método de estimación. 10 3. Marco teórico A continuación, se presenta el análisis bibliográfico, el cual describe técnicas novedosas a considerar durante el desarrollo de la investigación, así como aplicaciones de los enfoques frecuentista y bayesiano. Yi & Ma (2012) realizan un estudio para identificar los marcadores genéticos significativamente asociados con el estado de supervivencia de los ratones, así como estimar los efectos genéticos de estos marcadores. Los datos utilizados corresponden a 116 ratones de un cruce entre cepas (BALB/cByJ y C57BL/6ByJ), además, los sujetos fueron genotipificados en 133 marcadores genéticos que abarcan 20 cromosomas. Cada uno fue infectado con “Listeria monocytogenes”, se tiene que el 30% de la población sobrevive hasta el final del experimento. La particularidad de estos análisis radica en que se cuenta con muchas variables predictoras, las cuales pueden ser agrupadas de manera natural para facilitar el análisis, así mismo, esto puede provocar grupos de variables altamente correlacionadas; por lo que en el análisis es ideal incorporar una estructura jerárquica en las variables de predicción. En este caso se propusieron dos métodos de penalización: LASSO adaptativo bayesiano y modelos jerárquicos bayesianos con dos distribuciones previas para los parámetros de varianza considerando una distribución Cauchy – que incorpora un parámetro que controla la contracción en las estimaciones de varianza –, y la distribución doble exponencial – esta cuenta con un factor que controla el encogimiento de los coeficientes –. Los principales resultados muestran que incorporar parámetros de ajuste en las distribuciones previas para las variables y grupos permite estimaciones confiables de los parámetros y aumentan la capacidad de determinar cuáles son las variables importantes. De ambas distribuciones sobresale la doble distribución exponencial, ya que tienen mayor capacidad de encontrar las variables importantes o significativas, además su tiempo de convergencia es menor. 11 En relación con los modelos multinivel, López et al. (2007) realizaron un análisis con datos longitudinales de registros del peso al nacimiento de ganado criollo sanmartiniano, las variables predictoras son el sexo de la cría (hembra o macho), número de partos de la madre, época de nacimiento (invierno o verano), edad y código del padre de la cría. El modelo planteado tiene una estructura jerárquica de cuatro niveles: (1) categoría de la edad de la madre, (2) familia a la cual pertenece, (3) padre y (4) peso del individuo. Son las variables sexo, época y número de partos los parámetros fijos, además, la variable padre se considera como factor aleatorio, y las respectivas interacciones entre estas variables también son consideradas aleatorias. Los resultados obtenidos en este estudio muestran que utilizar el Mejor Predictor Lineal Insesgado (BLUP – por sus siglas en inglés) para las estimaciones es muy útil en modelos mixtos jerárquicos, mas no así en modelos con estructuras longitudinales. Por otra parte, Barnard et al. (2000) plantean la modelación de una matriz de covarianza en términos de sus desviaciones estándar 𝑆, donde 𝑆 es un vector 𝑘 ∗ 1 y la matriz de correlación 𝑅, donde 𝑅 es una matriz 𝑘 ∗ 𝑘. Los autores suponen dos escenarios: (1) usar estimadores de coeficientes shrinkage, (2) usar un modelo de localización general (GLOM). La alternativa de analizar esta matriz resulta relevante, ya que la estimación de la matriz de covarianza puede ser complicada debido a su dimensionalidad y restricción de ser definida no negativa. Los supuestos considerados en el análisis son: (1) la matriz de desviaciones 𝑆 y la matriz de correlaciones 𝑅 son independientes, (2) no se tiene mucha información previa, (3) la información previa con la que se cuenta intencionalmente refleja débil conocimiento de la matriz 𝑅. En relación con el uso de coeficientes shrinkage, se parte de que no se cuenta con información previa para los parámetros 𝛽𝑗. Además, 𝛽𝑗~𝑖. 𝑖. 𝑑. 𝑁(𝛽, ∑) donde 𝛽 representa el vector de distribuciones previas que puede ser ajustado y ∑ corresponde a la varianza, ambos fijos y conocidos. 12 Los autores utilizaron poca información previa considerando un registro independiente normal de cada una de las desviaciones estándar. Para el primer coeficiente, su valor verdadero en el modelo de regresión oscila entre 1 a 2.8 y la desviación estándar de la muestra es de 0.6, a partir de ello, escogieron las previas de la siguiente manera: Entonces, para 𝑠1 consideraron: 𝑙𝑜𝑔 (𝑠1) ~ 𝑁 (−4, 0.62), en donde se tiene una creencia previa de que los valores del primer coeficiente en las diez regresiones están mucho más juntos de lo que realmente están. Para 𝑠2 se considera: 𝑙𝑜𝑔 (𝑠2) ~ 𝑁 (0, 1), por lo tanto, los valores están más dispersos que en 𝑠1, finalmente para 𝑠3 se replica la distribución de 𝑠2. Se espera que los resultados muestren una mayor contracción en las estimaciones para 𝑠1 y menor contracción en 𝑠2 y 𝑠3. Para el segundo escenario, Barnard et al. (2000) proponen una generalización del GLOM mediante la descomposición de la matriz de covarianza en sus correlaciones y desviación estándar. Este mecanismo permite una mayor flexibilidad al especificar la estructura de covarianza condicional de las variables continuas, y la estimación de contracción de las covarianzas bajo un modelo jerárquico bayesiano. El GLOM utilizado va acorde a la propuesta de Lui y Rubin7 (1998) citado por Barnard et al. (2000). En este caso se parte de la siguiente matriz de covarianzas: 𝑉𝑖 = 𝑑𝑖𝑎𝑔(𝑆𝑖) 𝑅𝑖 𝑑𝑖𝑎𝑔(𝑆𝑖) 7 Estos autores presentan una extensión del modelo GLOM que permite la estimación de una matriz de covarianza de la siguiente manera: ∑ 𝑖 = 𝜆𝑖∑ en la cual, la media geométrica de las 𝛌 se establece en 1. Este modelo supone que las normales multivariantes dentro de cada celda tienen la misma forma elipsoidal, pero posiblemente diferentes tamaños. También permiten restricciones logarítmicas adicionales sobre los coeficientes de proporcionalidad 𝜆𝑖, lo que permite poner más estructura en términos de covariables. Esta estimación presenta como ventaja una mayor flexibilidad en las estimaciones, manteniendo la cantidad de parámetros desconocidos en un tamaño manejable. 13 Donde 𝑆𝑖 es el vector de desviaciones estándar y 𝑅𝑖 la matriz de correlaciones, además, 𝑖 = 1, … , 𝑐 donde 𝑐 es cada celda de la matriz. Para este caso se asumió que 𝑆𝑖 no tiene restricciones y 𝑅𝑖 = 𝑅 para todo 𝑖. Para los cálculos y debido a la estructura de los datos se utilizó el muestreo de Gibbs. Los autores concluyen que la técnica de modelar la matriz de covarianza en términos de las desviaciones estándar y la matriz de correlación es una estrategia común cuando la restricción para la matriz de correlación, de ser definida positiva, es fácil de tratar utilizando el muestreo de Gibbs. Otra propuesta de análisis para las matrices es la descomposición de Cholesky8, mencionada por varios autores, entre ellos Pan y Huang (2014). En este caso se realiza la selección de efectos aleatorios en modelos lineales generalizados mixtos mediante la función de penalización por encogimiento o contracción, para ello plantean una reparametrización de la matriz de covarianza mediante la descomposición Cholesky, en la cual se agrega un término de penalización por contracción a la función de cuasi verosimilitud penalizada (PQL9) en los componentes de varianza. Los métodos de penalización por contracción tienen la ventaja de seleccionar los efectos aleatorios “efectivos”, es decir los que contribuyen con información al modelo, eliminando las covariables que resultan innecesarias en términos informativos; y, para estimar los parámetros asociados, esta penalización es considerada una función de las varianzas de los efectos aleatorios. La principal ventaja de los métodos de penalización por contracción es que, para las variables 8 La factorización o descomposición de Cholesky es una matriz simétrica definida positiva, la cual se puede desagregar como el producto de una matriz triangular inferior (denominada el triángulo de Cholesky) y la traspuesta de la matriz triangular inferior. (Mora-Escobar, H.,2011) 9 La cuasiverosimilitud penalizada (PQL) es propuesta por Breslow y Clayton 1993 en “Approximate inference in generalized linear mixed models”. (Pan y Huang, 2014) https://es.wikipedia.org/wiki/Matriz_sim%C3%A9trica https://es.wikipedia.org/wiki/Matriz_definida_positiva https://es.wikipedia.org/wiki/Matriz_triangular https://es.wikipedia.org/wiki/Matriz_traspuesta 14 significativas al modelo, la penalización es baja o casi nula. Simultáneamente, ahorran tiempos en el cálculo, y son más estables en la estimación del modelo. En relación con el cálculo, se mencionan las ventajas en la estimación PQL, además de funciones de encogimiento como la LASSO10 y SCAD11. En ambos casos estas se encargan de penalizar a las variables no informativas del modelo, de modo que reducen a cero los coeficientes de regresión asociados a dichas variables, por ello las variables significativas se mantienen en el modelo con una leve contracción en la estimación. Los autores realizan simulaciones para evaluar el método propuesto considerando cuatro distribuciones: normal12, bernoulli13, binomial14 y poisson15, 10 Este método de estimación es denominado “Least Absolute Shrinkage and Selection Operator” (LASSO por sus siglas en inglés) y fue propuesto por Robert Tibshirani (1996). 11 El método de estimación es denominado “Smoothly Clipped Absolute Deviation” (SCAD – por sus siglas en inglés) y fue propuesto por Fan y Li (2001). Por otra parte, puede considerarse como una modificación de LASSO, y tiene las propiedades como la continuidad, dispersión e imparcialidad. 12 Este modelo supone que las variables independientes para los efectos fijos y aleatorios siguen un orden en el tiempo. Para cada estimación se hacen 100 simulaciones. Como principales resultados se tiene que un pequeño número de efectos aleatorios es seleccionado de manera incorrecta, en las simulaciones se da una selección efectiva de los efectos entre el 97% y 98%, y as estimaciones de los parámetros tienen un sesgo pequeño. 13 Consideró un modelo similar al caso de la distribución normal, con n = 100 y 10 covariables dependientes en el tiempo, las respuestas son binarias distribuidas con el correspondiente logit de probabilidad de éxito. En este caso la selección de componentes aleatorios baja en comparación con la distribución normal, lo que se interpreta como un ajuste insuficiente o excesivo. 14 Los ajustes de los modelos son semejantes al caso Bernoulli, excepto que esta vez la variable de respuesta provienen de distribuciones binomiales 𝐵(𝑚, 𝑝𝑖𝑗) donde 𝑚 es el número de ensayos y 𝑝𝑖𝑗 es la probabilidad de éxito de la prueba. En este caso, el método propuesto para la selección de efectos aleatorios funciona muy bien en el caso de covariables dependientes en el tiempo. En términos de comparación entre dos casos discretos de datos, el análisis con distribución binomial supera al análisis de Bernoulli, lo cual era esperado por los autores, ya que los datos binomiales proporcionan más información. 15 La configuración de parámetros generales y el mecanismo de generación de datos es similar a los casos de datos normales, bernoulli y binomiales con covariables dependientes en el tiempo. Sin embargo, las respuestas son generadas según las distribuciones de Poisson, y las covariables fueron generadas a partir de una distribución normal multivariante. En este caso los autores concluyen que el método de selección de efectos 15 considerando las variables independientes y dependientes en el tiempo, en este caso solo se presentan los casos de variables dependientes en el tiempo ya que es similar a la estructura de datos que se desea analizar. Cai and Dunson (2006) también utilizan la descomposición de Cholesky al realizar un análisis para conocer cuáles son los factores que afectan la fecundidad de las mujeres, además determinar si estos varían de mujer a mujer. Para ello proponen un modelo lineal generalizado mixto (GLMM) de enfoque bayesiano, en el cual se incorpora la descomposición de Cholesky a la problemática de selección de la matriz de covarianza de los efectos aleatorios, además calculan las probabilidades marginales posteriores de incluir cada predictor en los componentes de efectos fijos y aleatorios de manera separada. Los datos utilizados por los autores, corresponden al tiempo de embarazo de un grupo de 427 mujeres entre 19 y 39 años de edad que estaban registradas en California Department of Consumer Affairs y que proporcionaron la información. La variable de respuesta es medida como un evento discreto de sobrevivencia16. Se cuenta también con información demográfica general de la mujer, datos reproductivos, historial de anticonceptivos y otros factores relacionados con la fertilidad. Para evaluar el comportamiento del procedimiento se realizaron simulaciones mediante la integración numérica y aproximaciones cadenas de Markov vía Monte Carlo (MCMC – por sus siglas en inglés) y el muestreo de Gibbs, para ello consideraron un total de 100 individuos con seis observaciones cada uno, y utilizaron diferentes estructuras de covarianza de los efectos aleatorios en un modelo lineal generalizado mixto con funciones de enlace: identidad, logística y aleatorios propuesto funciona muy bien para los datos de distribución de Poisson con covariables dependientes del tiempo. 16 En el caso de que la mujer registre el 13avo ciclo menstrual, ese registro se censura en el análisis. Cai and Dunson (2006) 16 logarítmica. Un burn-in de 2000, y 20 000 iteraciones para las estimaciones, además de un análisis de sensibilidad con cinco diferentes “semillas” o valores de inicio. De los principales resultados obtenidos en dicho estudio, los autores concluyen que los modelos bayesianos son una buena alternativa en el análisis, la varianza de efectos aleatorios pequeña da buenas aproximaciones y los métodos resultaron robustos en el análisis de sensibilidad. Sin embargo, se hace énfasis en que el enfoque expuesto puede resultar más útil en casos con menos de 20 predictores. Chen (2003) realiza una reparametrización para obtener una estructura condicionalmente lineal al transformar los coeficientes aleatorios en coeficientes de regresión lineal para facilitar el uso de distribuciones previas conjugadas normales. La reparametrización propuesta consiste en seleccionar un subgrupo de elementos del vector de componentes aleatorio, además suponiendo que el componente de efectos fijos es conocido. Para los elementos del vector aleatorio no seleccionados se asume que tienen varianza 0. Se elige también una matriz diagonal no negativa y otra matriz triangular con 1´s en las entradas diagonales, logrando la descomposición de la matriz de covarianza de los efectos aleatorios. Los datos utilizados por Chen (2003) son de 12 centros de salud que participan en el Proyecto Perinatal Colaborativo (CPP) entre 1959 y 1966. Se desea conocer el desarrollo psicomotor de un niño de 8 meses de edad ante la exposición prenatal a bifenilos policlorados (PCB) y el desarrollo motor en niños pequeños, además de investigar si existen diferencias entre los centros de salud. La variable de respuesta es el nivel de exposición que se mide del suero materno tomado en el tercer trimestre del embarazo, las variables independientes son triglicéridos del tercer trimestre y nivel de colesterol sérico, raza de la madre (blanca, negra), educación de la madre (primaria o menos, escuela secundaria, más de secundaria), el orden de nacimiento del niño (primero o no) y un indicador de si el niño alguna vez fue amamantado. 17 Para la estimación se utiliza un modelo jerárquico bayesiano, y son seleccionadas distribuciones previas que facilitan las especificaciones computacionales, además, se establecieron tres cadenas con valores iniciales distintos, considerando 50 000 iteraciones iniciales y 100 000 iteraciones para las estimaciones posteriores, y se realizó también un análisis de sensibilidad. Los principales resultados indican que la parametrización en el modelo mixto facilitó la especificación y cálculo mediante el muestreo de Gibbs, obteniendo resultados robustos. Por otra parte, se unificó el problema de selección de variables y efectos aleatorios. Finalmente, el uso de distribuciones no normales, como la distribución t de Student en los términos de error, o efectos aleatorios, permitió obtener resultados robustos. Huang et al. (2006) analiza la efectividad de cierto fármaco que reciben algunos sujetos para medir el comportamiento del virus VIH. El objetivo de la investigación corresponde por una parte a la estimación de ciertos parámetros a nivel del individuo, y a un análisis en datos longitudinales a largo plazo con datos desbalanceados. Los datos utilizados corresponden a la simulación de un ensayo clínico con 20 pacientes portadores del VIH que reciben tratamiento antiviral a largo plazo, las mediciones se realizan cada 25 días desde el día 0 hasta el día 200 del seguimiento. Además, Huang et al. (2006) mencionan que los modelos comúnmente utilizados para este tipo de análisis tienen una limitante importante: la relación que se puede establecer entre la variable de respuesta y los factores independientes es a corto plazo. Además, se sabe qué variable dependiente tiene algunos patrones similares entre los pacientes aun conservando sus características individuales, por ello los modelos jerárquicos de efectos mixtos no lineales (NLME) parecen ser los más indicados. Para el análisis se utiliza un modelo de efectos mixtos no lineal bayesiano (BNLME) debido a que los programas computacionales no eran capaces de realizar las estimaciones de NLME. 18 Dentro de los parámetros de este modelo se incorpora la información de marcadores de fenotipo que permiten cuantificar ciertos fármacos en el individuo a través del tiempo, los parámetros de un modelo de exposición a fármacos, la resistencia del virus, la eficacia antiviral del fármaco y la interacción entre la información de las células susceptibles a la infección del VIH, de células infectadas y del virus libre. La especificación del modelo se realiza en tres etapas: (1) se definió la variación dentro del conglomerado, (2) se definió la variación entre conglomerados, (3) se especificaron las distribuciones previas. Por otra parte, los valores de los hiperparametros se obtienen de estudios anteriores o de literatura de referencia, en los casos donde se carece de información son asignadas distribuciones previas no informativas con varianzas grandes. Las estimaciones son realizadas mediante las Cadenas de Markov Vía Monte Carlo (MCMC) combinando los métodos de Gibbs y Metropolis-Hastings, se consideran 30 000 iteraciones iniciales y 120 000 para la estimación, además se realiza un análisis de sensibilidad. Dentro de los principales resultados se tiene que los modelos proporcionaban un buen ajuste a los datos, las estimaciones de sesgo para la población y los parámetros individuales fueron muy pequeñas, por su parte las correspondientes desviaciones estándar de las estimaciones se consideraron razonables. Así mismo, el modelo propuesto ajustó los datos moderadamente bien para la mayoría de los sujetos en el estudio, en menos del 10% se tuvo patrones de respuesta inusuales a la carga viral, mediciones inexactas de la exposición a fármacos y/o adherencia para estos sujetos. Por otra parte, Kinney y Dunson (2006), realizan un análisis para la selección de efectos aleatorios en un modelo mixto logístico mediante la implementación del enfoque bayesiano con el muestreo de Gibbs, abordan principalmente dos problemáticas: definición de las distribuciones previas subjetivas, así como la ineficiencia computacional pues se tiene lentitud en las estimaciones. 19 Los autores se centran en “Stochastic Search Variable Selection” (SSVS) iniciando con el modelo que contiene todas las variables, buscando los modelos con alta probabilidad posterior, y de esta manera se logra identificar el subconjunto de predictores que tienen coeficientes no nulos o varianzas diferentes de cero para aplicar el muestreador de Gibbs. También se considera la descomposición de Cholesky para la especificación de las matrices a utilizar en las estimaciones. Las distribuciones previas para los componentes de los efectos fijos son Bernoulli y Gamma, y en los componentes aleatorios las distribuciones son similares. Por otra parte, se menciona como posibles alternativas para el análisis los modelos logísticos de datos ordinales y probit. La técnica se ilustra con el análisis de datos de epidemiología y datos simulados. Los datos pertenecen al Proyecto Perinatal Colaborativo (CPP), la recolección de la información se realizó entre 1959 y 1966 con un total de 5389 participantes. Se mide el efecto del DDE (Dicloro Difenil Dicloroetileno), un metabolito del DDT (Dicloro Difenil Tricloroetano), el cual es medido en el suero materno luego de una pérdida de embarazo y es considerada una variable binaria. Las covariables son la edad de la madre, el índice de masa corporal, el tabaquismo, los niveles de colesterol y triglicéridos. Para la estimación de efectos aleatorios se considera que los 12 centros de los cuales se obtuvo la información son heterogéneos entre sí. El objetivo es seleccionar un modelo logístico de efectos mixtos que relacione los niveles de DDE y la pérdida de embarazo, en el cálculo de las estimaciones se consideran 30 000 iteraciones y un burn-in de 5000 iteraciones. En la simulación se considera un modelo logístico de respuesta binaria con tres covariables de distribución uniforme para 30 observaciones de 200 sujetos, además se eligió un rango de valores posibles o realistas para las varianzas de efectos aleatorios y se realizó un análisis de sensibilidad. Los resultados obtenidos son comparados con el enfoque de cuasi verosimilitud penalizada, ya que es muy utilizado en las estimaciones de GLMM. El enfoque 20 propuesto da ventajas en la selección de los efectos fijos y aleatorios pues lo realiza al mismo tiempo. Por su parte el modelo bayesiano permitió el desarrollo computacional para el cálculo de las probabilidades marginales, lo cual es de utilidad en el muestreo de Gibbs. En el análisis de sensibilidad, se obtuvo resultados robustos con poca diferencia en las estimaciones de los parámetros. También, Bondell et al. (2010) realizan una simulación para identificar de forma simultánea los principales predictores para los componentes fijos y aleatorios en un modelo mixto. Para el análisis se usan los datos de la Red de Tendencias y Estado del Aire Limpio (CASTNet) que contiene el registro de los lugares de Estados Unidos con medidas repetidas de contaminación y variables meteorológicas. Estos datos permiten realizar simulaciones más precisas de la calidad del aire. La propuesta consiste en realizar la selección simultánea de factores fijos y aleatorios utilizando una descomposición de Cholesky modificada, mediante una reparametrización de los modelos lineales mixtos. Para el análisis se parte del supuesto de distribución normal en la distribución condicional y la distribución de los efectos aleatorios, en caso de que este supuesto no se cumpla, se puede ocurrir falta de robustez. Para las penalizaciones se propone el método LASSO ajustado. En la selección del parámetro de ajuste se considera el algoritmo de Expectation- Maximization17 (EM) en el cual se supone que los efectos aleatorios no son observados y se le asigna un valor fijo al parámetro de regularización no negativo utilizado para estimar el modelo de regresión lineal LASSO ajustada. Los resultados obtenidos de la simulación muestran que el rendimiento de los métodos típicos que seleccionan separadamente los componentes fijo y aleatorio, 17 Utilizado por Laird y Ware en “Random-effects models for longitudinal data”, 1982 y Laird, Lange y Stram en “Maximum likelihood computations with repeated measures: application of the EM algorithm”, 1987. (Bondell et al., 2010) 21 como el criterio BIC o GIC18 extendido, no es tan bueno como el método propuesto que selecciona simultáneamente los componentes fijos y aleatorios. Se tiene también el caso de un ensayo clínico con datos hipotéticos, el cual tiene como objetivo analizar un nuevo medicamento para combatir la depresión. Los intereses de la investigación son enunciados a continuación: 1) Lograr evidencia de la efectividad del medicamento. 2) Establecer si el efecto del medicamento se maximiza o minimiza con el paso del tiempo. 3) Establecer si hay evidencias para afirmar que el efecto del medicamento tenderá a ser el mismo en todos los individuos con estos mismos criterios de inclusión. Para ello se cuenta con un grupo de 24 sujetos, los cuales son asignados aleatoriamente a un grupo de control (12) y un grupo de tratamiento (12). Los sujetos del grupo de control reciben un placebo y los demás reciben el medicamento. Se realizan cuatro mediciones, utilizando una escala, la primera medición corresponde a la línea base, es decir antes de administrar el tratamiento, las demás mediciones se realizan 1 mes, 3 meses y 6 meses después de la línea base. En análisis se realiza en dos escenarios, en el primero se considera la base de datos completa y en el segundo son eliminados algunos de los valores para realizar las estimaciones con valores faltantes, en este último caso se tienen distintas razones por las cuales los valores fueron eliminados. Se utilizan tres modelos multinivel, el primero corresponde a un modelo de efectos fijos, el segundo modelo incorpora efectos fijos y aleatorios, y el tercer modelo incorpora además de los efectos fijos y aleatorios, las correlaciones entre los efectos aleatorios. Los principales resultados indican que para la variable grupo de tratamiento la estimación del parámetro cambia a medida que se incorporar efectos aleatorios en 18 El GIC es un criterio de selección similar a los criterios AIC, BIC y AIC condicional, es mencionado por Rao y Wu en “A strongly consistent procedure for model selection in regression problems”,1989. (Bondell et al., 2010) 22 los modelos, la estimación de la constante se mantiene similar en los tres modelos. Además, entre el segundo y tercer modelo estimado, se presentan cambios en los parámetros correspondientes a la variable tiempo y en los errores estándar de los parámetros fijos del modelo; en el caso del tercer modelo que incorpora la estimación de las covarianzas que pierde significancia estadística. Ntzoufras y Lykou (2013) proponen un método bayesiano para lograr la reducción y selección de variables con la regresión LASSO, el problema que se trata de resolver inicialmente es la especificación del valor del parámetro 𝛌 a través de los factores de Bayes y utilizando el método de MCMC para realizar las estimaciones. Se le asocia a este parámetro los valores de Pearson y la correlación parcial en los límites entre significación e insignificancia según lo definido por los factores de bayes, y se realiza una comparación entre la regresión LASSO ordinaria y la bayesiana. Se realizan dos estudios de simulación y la aplicación del modelo a un conjunto de datos reales, además se considera el escenario en el cual se tiene un número de predictores menor al número de observaciones. Para identificar fácil y de manera clara el valor del parámetro de contracción, los autores parten del nivel de correlación identificado entre en el límite entre las variables significativas e insignificantes para el factor Bayes, de la siguiente manera: (1) Se traza un rango de valores de correlación entre las covariables y variable de respuesta, independientemente del valor de 𝛌. Así a partir de las medidas de correlación de Pearson se tendrán factores de Bayes que medirán la evidencia a favor de la adición de esta covariable al modelo (inferior a uno para todos los valores de 𝛌). Se identifica también el rango de valores de correlación donde la inclusión de esta covariable nunca es apoyada a posteriori con la suficiente firmeza; es decir, el factor de Bayes correspondiente nunca llega a ser más alto que un nivel específico α para todos los valores de 𝛌. (2). Se especifica el valor de 𝛌, partiendo de los siguientes escenarios: (2.1) considerando a 𝛌 como constante, en este caso se especifica 𝛌 definiendo los 23 niveles de importancia práctica (es decir, cuando los factores de bayes son iguales a uno) según las medidas de correlación producidas en el paso 1 – o variable aleatoria –, (2.2) considerando a 𝛌 como aleatorio: se toma el rango de correlaciones producido especificando un hiperparametro sensible para 𝛌. Se tiene que, a partir de la especificación de los hiperparametros del paso inicial, se evita el uso de distribuciones previas no informativas para el parámetro de contracción 𝛌. Los autores mencionan que conforme aumenta el 𝛌 la distribución se vuelve más informativa. Y, para valores pequeños de 𝛌 no se implementa la contracción, mientras que 𝛌 aumenta todos los coeficientes se reducen a cero. Dentro de las distribuciones previas contempladas se tiene la doble exponencial, para cada coeficiente 𝛽𝑗 del modelo, con ventajas en la selección de variables ya conocidas en el enfoque bayesiano. Se realiza un análisis de sensibilidad en los factores Bayes en la elección de los parámetros de contracción cuando se realiza la regresión de LASSO múltiple. Se exploran las ventajas de los métodos de selección de variables y encogimiento. La contracción se logra mediante el uso de un producto de distribuciones previas dobles independientes, e independientes para los coeficientes de regresión. La selección de variables se logra a través de los indicadores de inclusión de variables binarias habituales incluidos en el predictor lineal. La estimación de las distribuciones posteriores (incluido el modelo posterior y las probabilidades de inclusión variable) se logra a través de un simple esquema MCMC. Finalmente, los resultados se comparan con una variedad de métodos: LASSO bayesiano partiendo de la correlación, LASSO bayesiano con hiperparametros gamma, logaritmo LASSO con selección del parámetro de contracción mediante validación cruzada, y otras como normal mixto con inversa gamma. Se concluyó que los métodos LASSO bayesianos funcionan eficazmente en todos los ejemplos, rastreando los efectos importantes con altas probabilidades posteriores y eliminando las covariables no importantes con bajas probabilidades posteriores. 24 4. Marco metodológico 4.1 Población bajo estudio Los datos para el análisis son fuente primaria de información, correspondiente a un experimento en Neurociencias, asociado a los proyectos de investigación: 723-B7- 610 “Motivación y plasticidad neuroconductual: Efectos cognitivos, emocionales y sociales del enriquecimiento ambiental aleatorio e impredecible en ratas” y 742-B6- 291 "Cambios en la plasticidad neuronal relacionados con la incubación de la sensibilidad por la anfetamina: evaluación molecular y estructural", estos se realizaron en el Centro de Investigación en Neurociencias de la Universidad de Costa Rica, en ambos casos es el señor Juan Carlos Brenes Sáenz el investigador principal, por lo cual en caso de alguna consulta referida a los datos, él es a quien se debe dirigir. La población bajo análisis comprende 31 ratas machos de la cepa Wistar. Cada rata tiene aproximadamente dos meses de edad y un peso entre 200 y 210 gramos. Los animales fueron albergados en un espacio previamente adaptado con alimento y agua, tienen ciclos de luz oscura entre las 7a.m. y 7p.m., así como una temperatura que oscila entre los 23°C y los 29°C. 4.2 Análisis de datos originales Los objetivos generales asociados a los experimentos de los cuales se desprende la información utilizada en esta investigación corresponden a: - Estudiar procesos relaci0onados con la plasticidad neuronal subyacentes a la incubación de la sensibilidad conductual por exposición a anfetamina, para contribuir al conocimiento sobre los cambios plásticos inducidos por drogas psicoestimulantes. - Determinar si el incremento en el valor motivacional del enriquecimiento Ambiental (EA) en ratas potencia los efectos estimulatorios del EA sobre el aprendizaje y la memoria, la interacción y comunicación social, la sensibilidad a la anfetamina y la expresión de genes asociados a la plasticidad neural. 25 Este modelo se desarrolla a partir de la evidencia de que la anfetamina puede producir síntomas similares a la manía (euforia e hiperactividad psicomotora) en individuos sanos y precipitar episodios o agravar la sintomatología en pacientes con trastorno bipolar (Smith y Daves 1977, Willson et al., 2005, citado por Pereira et al. 2013); y debido a que la hiperactividad, la euforia y verbosidad excesiva son tres de los síntomas típicos de la anfetamina en humanos y de los episodios de manía dentro del trastorno bipolar (Van Kammen y Murphy, 1975; Willson et al., 2005, citado por Pereira et al. 2013), se partió del supuesto de que las vocalizaciones ultrasónicas (USVs) inducidas por la anfetamina con una frecuencia promedio superior o igual a los 50 kilohertzios (kHz) pueden ser consideradas como un marcador conductual para modelar la manía en animales de laboratorio. Los fármacos que se desean analizar corresponden al tamoxifeno, el litio y la miricitrina, de los cuales solo el litio es un medicamento clínicamente aprobado para tratar el trastorno bipolar, los otros dos están en proceso de validación y gran parte del estudio era mostrar los potenciales efectos antimaniáticos de dichos fármacos. Durante el experimento se realizó un tamizaje midiendo la cantidad de USVs espontáneas emitidas por las ratas al ser colocadas individualmente en una caja con borucha fresca. A partir de esta información, los sujetos fueron asignados aleatoriamente a dos grupos, pero de manera balanceada para que la variabilidad poblacional de la tasa de USVs estuviera representada lo más homogéneamente posible entre de los grupos. El análisis de datos se realizó mediante un análisis de varianza múltiple (MANOVA), considerando una significancia de p<0.05. Se comparó cada sustancia suministrada con el respectivo grupo de control, es decir: SAL+SAL vs. SAL+AMP, LI+SAL vs. LI+AMP y TAM+SAL vs. TAM+AMP, y se analizaron las vocalizaciones emitidas, así como la actividad locomotora. El experimento consistió en colocar a cada uno de los sujetos de manera independiente durante tres días consecutivos en un campo abierto (CA) de acrílico 26 de 40 cm x 40cm x 40 cm, con el piso cubierto de borucha, bajo luz roja y con un micrófono ultrasónico colocado a 45 cm sobre el suelo de la caja por un periodo de 10 minutos cada día. Las USVs se monitorearon con un micrófono ultrasónico (UltraSoundGate Condenser MicrophoneCM16) y gravados con el software Avisoft Recorder 2.7. Al finalizar el tiempo del test, posterior a la extracción de cada rata y antes de ingresar a la siguiente, se cambió la borucha del CA. En la ejecución del experimento (considerando únicamente la aplicación de las sustancias salina y anfetamina): día 1, se realizó una línea base para tener la actividad locomotora y las USVs espontánea de los animales previo a los tratamientos. En el día 2 se les suministró 0.9% de solución salina (sub cutáneo) a todos los sujetos 15 min antes de colocarlos en el CA. En el día 3 los sujetos reciben dos inyecciones previo al test: la primera inyección se aplica 25 min antes del CA, en la que se les suministró a todos los sujetos solución salina; diez min después (15 minutos antes del CA) se aplicó la segunda inyección, a un subgrupo de los sujetos se les suministró salina (16 ratas) y al otro subgrupo anfetamina (2.5mg/kg) (15 ratas). Además, se midió la cantidad de vocalizaciones emitidas según los cuatro subtipos19 en cada minuto del experimento. 19 Los subtipos de vocalizaciones son definidos por Pereira et al. (2014), de acuerdo al tono, frecuencia y forma. Dicha clasificación se realizó manualmente a partir de espectrogramas. Los tipos de vocalización son los siguientes: a. Flat calls (planas): cuando los cambios de frecuencia dentro de un único elemento de cada USV eran iguales o inferiores a 5-kHz, pero la diferencia en la frecuencia de inicio y final sí pudo ser superior a 5-kHz. b. Step-calls. (escalonadas): son USV de frecuencia modulada (FM) y ocurren cuando la vocalización flat tiene al menos un elemento plano corto superpuesto al inicio y/o al final de la misma. c. Trills (zigzag): son USV de FM con cambios de frecuencia superiores a 5-kHz o con dos o más cambios de frecuencia en direcciones opuestas d. Mixed calls (mixtas): son USV de FM que no se incluyeron en las categorías anteriores y que pueden tener uno o más componentes planos y/o zigzag superpuestos. 27 Finalmente, por medio del programa Avisoft Bioacustics se midió la frecuencia (en Hercios), intensidad (en decibeles) y duración promedio (en segundos) para cada vocalización. 4.3 Análisis de datos en esta investigación Para el presente análisis se debe aclarar que, por un tema de confidencialidad de los proyectos y la información recolectada en los mismos, solo se tuvo acceso a los datos provenientes de animales tratados con solución salina-salina y salina- anfetamina. Además, no se incluye la variable tipo de vocalización, densidad ni duración de la misma, solamente la información correspondiente a la cantidad de vocalizaciones. Debido al interés del estudio, no se consideró la información de cada registro de vocalización, sino que se hizo un conteo de las vocalizaciones emitidas para cada uno de los sujetos en cada minuto de los diez que duró el experimento, por otra parte, por criterio y recomendación del experto, se utilizó la información de los minutos 1, 4, 7 y 10 del CA en los tres días. Por lo tanto, para cada sujeto se tienen mediciones cada día en el que se desarrolló el experimento, y a su vez, para cada día del experimento, se tiene la información de los cuatro minutos seleccionados, es decir, cada rata es un conglomerado. Así, se tiene una estructura de tres niveles, en la cual el primer nivel corresponde a las mediciones de las vocalizaciones en los minutos indicados anteriormente y se cuenta con un total de 372 registros. El segundo nivel corresponde a las mediciones en los tratamientos recibidos por los sujetos, en los tres días, para un total de 93 registros; finalmente, el tercer nivel corresponde a las 31 ratas o conglomerados. La estructura de los datos se presenta en la imagen 1. 28 Imagen 1 Diagrama de la estructura de los datos Día 1: sin tratamiento Día 2: solución salina Día 3: solución salina o anfetamina M1: medición de vocalizaciones en el minuto 1 – día 1 M2: medición de vocalizaciones en el minuto 4 – día 1 M3: medición de vocalizaciones en el minuto 7 – día 1 M4: medición de vocalizaciones en el minuto 10 – día 1 M5: medición de vocalizaciones en el minuto 1 – día 2 M6: medición de vocalizaciones en el minuto 4 – día 2 M7: medición de vocalizaciones en el minuto 7 – día 2 M8: medición de vocalizaciones en el minuto 10 – día 2 M9: medición de vocalizaciones en el minuto 1 – día 3 M10: medición de vocalizaciones en el minuto 4 – día 3 M11: medición de vocalizaciones en el minuto 7 – día 3 M12: medición de vocalizaciones en el minuto 10 – día 3 Se debe considerar que debido a que en el día tres del experimento los sujetos son asignados a distintos tratamientos – 16 reciben salina y 15 anfetamina – el diseño del conjunto de datos es desbalanceado en el tercer día. Para las variables categóricas tratamiento y minuto se utiliza la notación de Gaviria & Castro (2005) para especificar el modelo mediante la codificación completa que emplea tantas variables dummies como categorías de la variable predictora. A continuación, son detallas las variables a utilizar en el análisis: Rata Día 1 M1 M2 M3 M4 Día 2 M5 M6 M7 M8 Día 3 M9 M10 M11 M12 29 Variable de respuesta: - Número de vocalizaciones: cantidad de vocalizaciones emitidas. Variables predictoras: - Día del tratamiento: día 1 – sin tratamiento, día 2 – salina, día 3 – salina o anfetamina. - Medición: en total se tienen 12 mediciones y corresponden al minuto en el cual se registró la variable de respuesta durante la ejecución el experimento 1, 4, 7 y 10, cada día. 4.4 Técnicas a emplear Para predecir la cantidad de vocalizaciones emitidas por las ratas, se utilizaron los modelos multinivel considerando estimaciones frecuentistas y bayesianas, en este caso puntual, el empleo de los modelos mixtos responde a la misma naturaleza jerárquica de los datos, pues las vocalizaciones son mediciones en distintos y determinados momentos del experimento, que a su vez pertenecen a una misma rata; en este caso específico corresponde a un modelo de medidas repetidas. En relación con los modelos mixtos, estos son modelos estadísticos para variables de respuesta continua, y en los cuales los residuos se distribuyen normalmente. A diferencia de los modelos de regresión clásica pueden no tener varianza constante o no ser independientes. Se tiene que, en casos donde los estudios tienen datos agrupados como clúster, diseños experimentales, estudios longitudinales o de medidas repetidas, son recomendados estos modelos pues le brindan al investigador una herramienta más flexible y poderosa para el análisis de los datos (West et al., 2007). En los modelos mixtos se pueden diferenciar dos tipos de parámetros, por una parte, están los fijos que corresponden a los efectos medios de la población y los aleatorios que corresponde a las varianzas y covarianzas de todos los niveles del modelo 30 (Gaviria & Castro, 2005, pág. 59). Estos permiten conocer diferencias entre los sujetos, ya sea considerando variables del conglomerado o de la unidad de análisis que conforma el conglomerado, por ende, omitir el uso de modelos mixtos bajo estas condiciones provocaría estimaciones erradas, tal como lo mencionan Gaviria & Castro (2005): Al ignorar la estructura de los datos, el problema que se produce es que eliminamos toda la varianza interna de los grupos, que puede llegar a ser el 80% o del 90%. Así las relaciones parecerán como muy fuertes, y pueden ser de hecho muy distintas de los resultados con las variables desagregadas. (pág. 12) También, los modelos multinivel resuelven el dilema de no independencia de los datos, que, al ignorar esta condición, puede mostrar resultados significativos cuando no lo son, por ende, las conclusiones estarían erradas. Los modelos mixtos: “constituyen una estrategia analítica que permite la formulación jerárquica de las fuentes de variación y con capacidad para dar cuenta de esta estructura” (Gaviria & Castro, 2005, pág. 14). Finch et al. (2014, pág. 37) señalan como supuestos de los modelos multinivel de dos niveles los siguientes: 1. En el segundo nivel los residuales son independientes entre clústers. 2. En el segundo nivel los interceptos y coeficientes son asumidos como independientes de los errores de primer nivel. 3. En el primer nivel los residuales están distribuidos normalmente y tienen varianza constante. 4. En el nivel dos el intercepto y las pendientes deben tener una distribución normal multivariada con una matriz de covarianza constante. Sin embargo, los modelos multinivel suponen que la distribución de la variable de respuesta es normal, por ende, para este caso de estudio al tener un variable de respuesta Poisson deben utilizarse los modelos lineales generalizados, los cuales permiten la estimación de modelos cuando en estos la variable de respuesta no tiene distribución normal. 31 Para la estimación frecuentista del modelo, se consideró la máxima verosimilitud restringida (RELM – por sus siglas en inglés) y no la conocida máxima verosimilitud (MLE – por sus siglas en inglés). Ello porque la RELM en el contexto multinivel incorpora un ajuste en los grados de libertad, ya que los considera en la estimación de los componentes de varianza, por ende, las estimaciones serán más apropiadas (Finch, 2014). Por otra parte, en relación con el enfoque bayesiano, su principal característica es incorporar en las estimaciones información previa que involucra el conocimiento del experto, o puede obtenerse también de estudios preliminares que tengan relación con el tema de interés. Esta información en conjunto con la información que generan los datos, conforma la denominada distribución a posteriori de la estimación bayesiana. La distribución previa puede ser informativa o no informativa: la primera hace referencia a información sustentada en resultados de investigaciones que logran identificar atributos de variables de interés, la no informativa que hace referencia a las situaciones en las cuales no se conoce del fenómeno lo suficiente. (Carlin & Louis, 2000) Las estimaciones del método bayesiano representan una alternativa a los métodos frecuentistas, se realizan mediante MCMC, el cual es un proceso iterativo donde la distribución a previa se combina con la información de la muestra actual y se genera la distribución a posteriori para cada uno de los parámetros del modelo. (Finch et al., 2014, pág. 169). Mediante un proceso repetitivo en el cual se actualiza la información en las cadenas de Markov, cuando cada cadena se estabiliza en un valor indicado se dice que la cadena converge y dicho valor corresponde a la estimación del parámetro en el modelo. Además, las aproximaciones de las distribuciones a posteriori resultan más precisas en la implementación del método de MCMC (Barnard, McCulloch & Meng, 2000), ya que la estimación mejora con el número de iteraciones; también hay una mayor flexibilidad al comparar modelos con una cantidad diferente de componentes de 32 varianza. Pero, este enfoque implica una mayor carga computacional para el investigador (Kinney y Dunson, 2006). Algunas de las ventajas de utilizar estadística bayesiana en la estimación de modelos multinivel son: (1). No supone una distribución de los datos, la estimación de los intervalos de credibilidad no se afecta por distribuciones asimétricas de los datos. (2). La estadística bayesiana puede ser de mucha utilidad cuando el modelo a estimar es muy complejo, pues en algunos casos las estimaciones frecuentistas pueden no lograr la convergencia. (3). Pueden estimar de manera precisa los parámetros con muestras pequeñas. (4). Las estimaciones bayesianas pueden ser utilizadas en los casos donde las estimaciones frecuentistas también funcionan. (Finch et al., 2014, pág. 168) Por otra parte, se debe indicar que en el contexto de experimentos neurocientíficos los tamaños de muestra tienden a ser pequeños, lo cual podría desfavorecer las estimaciones en los modelos estadísticos frecuentistas, no así en los modelos de enfoque bayesiano. Algunos ejemplos de estudios con tamaños de muestra reducidos son: Brenes at al. (2016) en “Differencial Effects of Social and Physical Environnmental Enrichment on Brain Plasticity, Cognition, and Ultrasonic Communication in Rats” realizan los análisis en 48 sujetos, por otra parte Brenes & Schwarting (2014) en su investigación “Attribution and Expression of Incentive Salience Are Differentially Signaled by Ultrasonic Vocalizations in Rats” tienen una población de estudio de 30 sujetos, otros casos aplicados que se relacionan con el tema son: “The Rat Pup Study” en Mixed-Effects Models in S and S-PLUS realizado por Pinheiro & Bates (2000) en el cual los datos fueron recolectados de 30 ratas, finalmente, Douglas et all (2004) en “Pontine and basal forebrain cholihergic interaction: implications for sleep and breathing” cuentan únicamente con cinco sujetos. Todos los anteriores relacionados en al área médica y con sujetos sometidos a análisis. 33 Respecto a la especificación de los modelos multinivel, el enfoque bayesiano es muy similar al enfoque frecuentista, los valores de las covarianzas y los errores pueden ser especificados en cada nivel y ser una fuente potencial acumulativa para la explicación del error en la variable de respuesta (Congdon, 2003, pág. 138). Para la estimación es importante identificar el tipo de variable dependiente en el análisis (continua, discreta, dicotómica, ordinal), pues a partir de ello se establecerán los lineamientos necesarios en la especificación del modelo que será utilizado. En este análisis, específicamente en la estimación bayesiana del modelo multinivel, se utilizaron previas no informativas. En relación con el fenómeno de encogimiento en los coeficientes fijos de los modelos mixtos, en la literatura se conocen como estimadores shrinkage a los ponderadores o valores utilizados para ajustar las estimaciones de los coeficientes. Así mismo, Royle & Link (2002) indican que los estimadores de contracción que surgen como consecuencia de la estimación (o predicción) de un mayor número de coeficientes aleatorios ya sea bajo consideraciones bayesianas o frecuentistas, lo que buscan es minimizar el error cuadrático medio de la predicción. También mencionan que, los coeficientes de contracción o encogimiento pueden mejorar la estimación de los coeficientes del modelo en los siguientes casos: donde se tenga un gran número de parámetros en el modelo, o donde se tiene tamaños de muestra pequeños. Gaviria & Castro (2005) señalan que el tamaño del shrinkage – considerado como un ponderador – es determinado por la cantidad de elementos en los niveles, de esta manera “conforme aumenta el número de unidades del primer nivel (𝑛𝑗), este término tiende a uno. Del mismo modo, cuando el número de unidades individuales en cada contexto desciende, este factor tiende a cero.” (pp. 71-72). Gaviria & Castro (2005) también mencionan que la varianza también es determinante para el factor shrinkage, pues si dentro de los conglomerados ésta es pequeña, es decir los sujetos son muy similares o hay muchos sujetos en el 34 conglomerado el factor shrinkage tiende a uno. Sin embargo, cuando los conglomerados tienen medias similares, las varianzas de los conglomerados son mayores que entre los conglomerados, es decir los grupos no son homogéneos, o la cantidad de sujetos es pequeña el factor de encogimiento tiende a cero. 4.5 Modelos analizados Inicialmente se realizó un análisis descriptivo de los datos, así como un análisis preliminar del fenómeno de shrinkage considerando distintos escenarios mediante un estudio de simulaciones. En el análisis descriptivo inicial, para conocer el comportamiento de los datos (anexo 1), se desagrega la variable cantidad de vocalizaciones por día y sustancia suministrada en el experimento, así como para cada sujeto. También, se creó una variable denominada “grupo” para identificar el comportamiento de los sujetos que recibieron anfetamina y los que no, y con ello generar insumos que permitieron determinar el mejor modelo que permita analizar los datos. A partir de las estadísticas descriptivas y, considerando que la variable de respuesta es un conteo mayor o igual a cero con distribución asimétrica positiva, aunado a ello, debido a que la varianza es mayor a la media, es decir, se tiene sobre dispersión se decidió utilizar un modelo lineal generalizado con distribución QuasiPoisson y función de enlace logaritmo natural. Los principales componentes de un modelo lineal generalizado son: (1) el componente sistemático es expresado como una relación lineal, (2) una función denominada función de enlace o link que indica la relación entre el valor esperado y el vector de covariables, (3) cuenta con un componente aleatorio que corresponde a la distribución condicional entre la variable de respuesta y el vector de covariables. En este caso puntual se tiene que la función de verosimilitud de la Poisson corresponde a: 35 𝐿(𝜇) = 𝑓(𝜇|𝑋) = 𝑒−𝜇 ∗ 𝜇𝒙 𝒙! Aplicando el logaritmo y la función se tiene: 𝐿𝐿(𝜇) = −𝑙𝑛(𝒙!) + 𝑙𝑛(𝑒−𝜇) + 𝑙𝑛(𝜇𝒙) 𝐿𝐿(𝜇) = −𝑙𝑛(𝒙!) − 𝜇 + 𝒙 ∗ 𝑙𝑛(𝜇) La función de enlace en el caso Poisson es 𝑙𝑛(𝜇). Es decir: 𝑙𝑛(𝜇) = 𝜷′𝒙 𝜇 = 𝑒𝑥𝑝(𝜷′𝒙) Finamente se tendría que la log-verosimilitud del modelo es: 𝐿𝐿(𝜇) = −𝑙𝑛(𝒙!) − 𝑒𝑥𝑝(𝜷′𝒙) + 𝑥 ∗ (𝜷′𝒙) Considerando todos los parámetros del modelo y de manera general, el predictor lineal del modelo puede ser detallado como sigue: 𝜂𝑖 = 𝛽0 + ∑ 𝑥𝑖𝑗𝛽𝑗 + ∑ ∑ 𝑧𝑖𝑗𝛽𝑗 𝐽𝑘 𝑗∈𝐺𝑘 𝐾 𝑘=1 𝐽0 𝑗=1 = 𝑿𝑖𝜷 Donde 𝛽0 es el intercepto, 𝑥𝑖𝑗 y 𝑧𝑖𝑗 representan los valores observados para las variables no agrupadas y agrupadas respectivamente, 𝛽𝑗 son los coeficientes del modelo. La notación 𝑗 ∈ 𝐺𝑘 indica el grupo de la variable j. Además, 𝑋𝑖 contiene todas las variables, 𝛽 es el vector que contiene todos los coeficientes, incluyendo el intercepto. De manera resumida, se tiene que 𝑋𝑖 = (1, 𝑥𝑖1, 𝑥12, … , 𝑥𝑖𝐽) y 𝛽 = (𝛽0, 𝛽1, … , 𝛽𝐽) ′ , donde 𝐽 = ∑ 𝐽𝑘 𝐾 𝑔=0 corresponde al total de variables. La media de la variable de respuesta es relatada mediante la función de enlace 𝑔: 36 𝐸(𝑦𝑖|𝑋𝑖) = 𝑔(𝑋𝑖𝛽) Finalmente, la distribución de los datos puede ser representada de la siguiente manera: 𝑝(𝑦|𝑋𝛽, ∅) = ∏ 𝑝(𝑦𝑖|𝑋𝑖𝛽, ∅) 𝑛 𝑖=1 Donde ∅ es el parámetro de dispersión, en el caso de las distribuciones de Poisson y Binomial, dicho parámetro toma valor 1. Para la selección de las variables que se incluyen en el modelo, se parte del conocimiento del experto en el tema de neurociencias así como algunos requerimientos para conservar la estructura jerárquica de los datos (ver anexo 2), para el análisis se utilizó la plataforma de programación R 3.5.2, los paquetes glmer20 Bates et a. (2015) para la estimación frecuentista y el paquete MCMCglmm21 Jarrod D Hadfield (2010) en el caso de la estimación bayesiana, además del programa STATA 15. (Para mayor detalle ver anexo 3.1 y 3.2 respectivamente). Para identificar la presencia de shrinkage en el modelo seleccionado, e iniciando con la estimación frecuentista, se estima un modelo de efectos fijos para las variables día y medición, seguido, se estima un modelo con tres coeficientes aleatorios22 y finalmente otro modelo con los siete efectos aleatorios de las variables 20 El glmer se utiliza en la estimación frecuentista de Modelos Mixtos Lineales Generalizados, utiliza el algoritmo de mínimos cuadrados penalizados. 21 Este paquete realiza la estimación en Modelos Lineales Generalizados Mixtos utilizando las técnicas de Cadenas de Markov Vía Monte Carlo, y utiliza como información previa a la distribución Wishart inversa para las (co)varianzas y una distribución normal para los efectos fijos. Utiliza un burn-in de 3000 iteraciones y 10 000 iteraciones para la estimación de los coeficientes. 22 Los coeficientes aleatorios corresponden las variables M1, M5 y M11. Las variables M1 y M5 corresponden a las primeras mediciones que son registradas al iniciar el experimento el día 1 y 2 respectivamente, estas toman relevancia ya que corresponde al registro de las primeras reacciones de adaptación del sujeto al CA, así las más recientes posterior a la aplicación de la inyección; por otra parte, la variable M11 corresponde a mediciones realizadas el tercer día del experimento, la que también resulta valiosa en el análisis pues es la que 37 día y medición, esto se repite en la estimación bayesiana. La estructura de los datos corresponde a: Imagen 2 Diagrama de la estructura de los datos analizada Donde se tiene que: Día 1: sin tratamiento Día 2: solución salina Día 3: solución salina o anfetamina M1: medición de vocalizaciones en el minuto 1 – día 1 M4: medición de vocalizaciones en el minuto 10 – día 1 M5: medición de vocalizaciones en el minuto 1 – día 2 M6: medición de vocalizaciones en el minuto 4 – día 2 M9: medición de vocalizaciones en el minuto 1 – día 3 M10: medición de vocalizaciones en el minuto 4 – día 3 M11: medición de vocalizaciones en el minuto 7 – día 3 Los modelos analizados son detallados a continuación: Modelo 1: de efectos fijos refleja el comportamiento del sujeto luego del proceso de adaptación al CA así como de las inyecciones, de igual manera corresponde a mediciones en el día que se aplica el tratamiento o sustancia de interés. Rata Día 1 M1 M4 Día 2 M5 M6 Día 3 M9 M10 M11 38 𝑦𝑖𝑗𝑘 = 𝛽0𝑗𝑘 + 𝛽1𝑗𝑘𝑋1𝑖𝑗𝑘 + 𝛽2𝑗𝑘 𝑋2𝑖𝑗𝑘 + 𝛽3𝑗𝑘 𝑋3𝑖𝑗𝑘 + 𝛽4𝑗𝑘 𝑋4𝑖𝑗𝑘 + 𝛽5𝑗𝑘 𝑋5𝑖𝑗𝑘 + 𝛽6𝑗𝑘 𝑋6𝑖𝑗𝑘 + 𝛽7𝑗𝑘 𝑋7𝑖𝑗𝑘 + 𝛽8 𝑋8𝑗𝑘 + 𝛽9 𝑋9𝑗𝑘 + ɛ𝑖𝑗𝑘 Modelo 2: incorpora coeficientes de efectos aleatorios para tres variables 𝑦𝑖𝑗𝑘 = 𝛽0𝑗𝑘 + 𝛽1𝑗𝑘𝑋1𝑖𝑗𝑘 + 𝛽2𝑗𝑘 𝑋2𝑖𝑗𝑘 + 𝛽3𝑗𝑘 𝑋3𝑖𝑗𝑘 + 𝛽4𝑗𝑘 𝑋4𝑖𝑗𝑘 + 𝛽5𝑗𝑘 𝑋5𝑖𝑗𝑘 + 𝛽6𝑗𝑘 𝑋6𝑖𝑗𝑘 + 𝛽7𝑗𝑘 𝑋7𝑖𝑗𝑘 + 𝛽8 𝑋8𝑗𝑘 + 𝛽9 𝑋9𝑗𝑘 + ɛ𝑖𝑗𝑘 𝑦𝑖𝑗𝑘 = (𝛽0 + 𝑉0𝑗 + 𝑉0𝑗𝑘) + (𝛽1 + 𝑉1𝑗𝑘)𝑋1𝑖𝑗𝑘 + 𝛽2𝑋2𝑖𝑗𝑘 + (𝛽3 + 𝑉3𝑗𝑘) 𝑋3𝑖𝑗𝑘 + 𝛽4𝑗𝑘𝑋4𝑖𝑗𝑘 + 𝛽5𝑗𝑘𝑋5𝑖𝑗𝑘 + 𝛽6𝑗𝑘 𝑋6𝑖𝑗𝑘 + (𝛽7 + 𝑉7𝑗𝑘) 𝑋7𝑖𝑗𝑘 + 𝛽8 𝑋8𝑗𝑘 + 𝛽9𝑋9𝑗𝑘 + ɛ𝑖𝑗𝑘 𝑦𝑖𝑗𝑘 = 𝛽0 + 𝛽1𝑋1𝑖𝑗𝑘 + 𝛽2𝑋2𝑖𝑗𝑘 + 𝛽3 𝑋3𝑖𝑗𝑘 + 𝛽4𝑗𝑘𝑋4𝑖𝑗𝑘 + 𝛽5𝑗𝑘𝑋5𝑖𝑗𝑘 + 𝛽6𝑗𝑘 𝑋6𝑖𝑗𝑘 + 𝛽7𝑋7𝑖𝑗𝑘 + 𝛽8𝑋8𝑗𝑘 + 𝛽9𝑋9𝑗𝑘 + (𝑉0𝑗 + 𝑉0𝑗𝑘 + 𝑉1𝑗𝑘𝑋1𝑖𝑗𝑘 + 𝑉3𝑗𝑘 𝑋3𝑖𝑗𝑘 + 𝑉7𝑗𝑘 𝑋7𝑖𝑗𝑘) + ɛ𝑖𝑗𝑘 Igualmente, en el segundo nivel se calcula la varianza correspondiente al intercepto 𝜎𝑉0 2 únicamente, para el tercer nivel se presenta la estructu