Introdución aos Modelos Mixtos

Traballo Fin de Grao INTRODUCIÓN AOS MODELOS MIXTOS Diego Losada González 2021/2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao INTRODUCIÓN AOS MODELOS MIXTOS Diego Losada González Xullo, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Estatística e Investigación Operativa Título: Introdución aos Modelos Mixtos Breve descrición do contido As estruturas xerárquicas de datos (estruturas multinivel) son frecuentes nas Ciencias Sociais, a Medicina ou a Bioloxía. Na formulación de estruturas xerárquicas asúmese que cada individuo pertence a un único grupo e o obxectivo é analizar as relacións a dous niveis: entre os grupos e dentro dos mesmos. A análise, tanto descritiva como inferencial, deste tipo de poboacións con estruturas complexas é o obxecto dos denominados modelos multinivel. Cabe sinalar que os modelos multinivel, terminoloxía que provén da Estatística Educacional, tamén se coñecen como modelos lineais xerárquicos ou modelos mixtos (Estatística, Bioestatística), modelos de efectos aleatorios ou modelos de coeficientes aleatorios (Econometría) ou modelos de compoñentes da varianza (Deseño de experimentos). Breve planificación: A modo de orientación, o traballo podería organizarse nas seguintes seccións: Modelo de análise da varianza: ANOVA Modelo de análise da varianza con efectos aleatorios: RANOVA Introdución aos modelos multinivel con resposta continua. iii iv Ademais, presentaremos diferentes modelos mixtos aplicados tanto a conxuntos de datos ou a datos simulados. Para iso, empregaremos o software estatístico libre R (https://www.r-project.org/). Recomendacións As estruturas xerárquicas de datos (estruturas multinivel) son frecuentes nas Ciencias Sociais, a Medicina ou a Bioloxía. Outras observacións Índice Resumo ix Introdución xi 1. ANOVA e ANCOVA 1 1.1. O modelo ANOVA ................................... 1 1.1.1. ANOVA como modelo linear xeral . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2. Estimación dos parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.3. Análise da varianza e test F . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.4. Comparacións múltiples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2. O modelo ANCOVA .................................. 12 1.2.1. ANCOVA seninteracción ........................... 12 1.2.2. Estimación dos parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.3. ANCOVA coninteracción ........................... 13 2. Introdución aos modelos mixtos 15 2.1. Datosmultinivel .................................... 15 2.2. A necesidade de ter en conta os distintos niveis . . . . . . . . . . . . . . . . . . . 17 3. RANOVA 19 3.1. Análise da varianza con efectos aleatorios . . . . . . . . . . . . . . . . . . . . . . 19 v vi ÍNDICE 3.2. Estimación dos parámetros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2.1. Estimación da media global µ......................... 21 3.2.2. Estimación das compoñentes da varianza . . . . . . . . . . . . . . . . . . . 22 3.2.3. Predición dos efectos aleatorios . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3. Contraste sobre os efectos grupais . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4. Modelos mixtos con covariables relativas ao primeiro nivel 35 4.1. Modelo con intercepto aleatorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2. Modelo con intercepto e pendente aleatorios . . . . . . . . . . . . . . . . . . . . . 40 5. Modelos mixtos con covariables relativas ao segundo nivel 47 5.1. Variables contextuais composicionais . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.2. Variables contextuais globais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.3. Interacción entre niveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6. Conclusións 55 A. Código de R 57 A.1. Creación da base de datos mates ........................... 57 A.2.Introdución ....................................... 58 A.2.1.Figura1..................................... 58 A.3. Creación da base de datos mates7 ........................... 59 A.4.ANOVA......................................... 60 A.4.1.Figura1.1.................................... 60 A.4.2. Función Bonferroni_taboa .......................... 61 A.4.3. Test F e contrastes pareados . . . . . . . . . . . . . . . . . . . . . . . . . . 63 A.5.ANCOVA ........................................ 63 A.5.1.Figura1.2.................................... 63 ÍNDICE vii A.6.RANOVA........................................ 65 A.6.1. VARCOMP e contraste sobre os efectos das escolas . . . . . . . . . . . . . 66 A.6.2.Figura3.5.................................... 66 A.6.3. Conxunto de datos u0df . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 A.6.4. Figura 3.3: Gráfico de eiruga . . . . . . . . . . . . . . . . . . . . . . . . . 67 A.6.5.Figura3.2.................................... 68 A.7. Modelos mixtos con covariables relativas ao nivel 1................. 71 A.7.1. Modelo con intercepto aleatorio e pendente fixa . . . . . . . . . . . . . . . 71 A.7.2. Modelo con intercepto aleatorio e pendente fixa con SocialMin ...... 71 A.7.3. Modelo con intercepto e pendente aleatorios e mais variable categórica . . 72 A.7.4.Figura4.3.................................... 72 A.8. Modelos mixtos con covariables relativas ao nivel 2................. 76 A.8.1. Variable contextual composicional . . . . . . . . . . . . . . . . . . . . . . 76 A.8.2. Variable contextual global . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Bibliografía 79 xiv INTRODUCIÓN εi Yi β0+β1x 0 5 10 15 20 −2 −1 0 1 Status socio−económico de 40 nenas/os Nota en Matemáticas de 40 nenas/os Figura 1: Diagrama de dispersión e recta axustada por mínimos cadrados (en azul) para a nota en Matemáticas de 40 nenas/os de 11 anos escollidas/os aleatoriamente fronte ao status socioeconómico destes. O segmento vertical vermello representa o residuo do dato i-ésimo (en violeta). Pódese dar un paso máis e extender o modelo linear simple mediante a consideración de máis dunha variable explicativa continua que tamén sirva para explicar a variable resposta Y. Este modelo de regresión denomínase modelo linear múltiple e pode expresarse como Y=X′β+ϵ,(2) onde o erro εsegue a satisfacer as hipóteses de homocedasticidade, normalidade e independencia supostas no modelo linear simple. De novo por mínimos cadrados pódese estimar o vector de parámetros β, obtendo b β= (X′X)−1X′Y, que segue unha distribución b β∈Np(β, σ2(X′X)−1)3, onde pé o número de variables explicativas a considerar. Finalmente, o modelo linear xeral non é máis que a xeneralización de (2) e abrangue a consideración tanto de variables explicativas continuas como discretas. Pode atoparse máis información sobre estes modelos na obra de Faraway [8]. Ao longo deste traballo tratarase de extender este tipo de modelos a situacións nas cales a variable resposta Yesté medida en diferentes grupos, como no Exemplo 0.1, onde se fornecen os resultados en Matemáticas de alumnado de diferentes centros educativos e non é estraño pensar que as diferentes escolas teñan un “efecto” nas notas acadadas polas/os estudantes do centro. Ademais, todo o código de empregado ao longo da totalidade do traballo atoparase dispoñible no repositorio Modelos_Mixtos_con_R de GitHub e será de carácter enteiramente reproducible. En gran parte atoparase tamén dispoñible no Anexo A. 3Z= (Z1, ... , Zm)∈Nm(µ, σ2Im)presenta unha distribución normal estándar multivariante se cada unha das súas compoñentes Zjten distribución normal estándar univariante e son mutuamente independentes. Capítulo 1 ANOVA e ANCOVA Ao longo deste capítulo analizarase o modelo de análise da varianza, tamén coñecido como modelo ANOVA (acrónimo de ANalysis Of VAriance), que non é máis que un modelo de regresión no cal hai unha única variable explicativa, que ademais é discreta ou categórica. Logo incluiranse variables explicativas discretas e continuas asemade para construír o modelo de análise da covarianza ou modelo ANCOVA (segundo terminoloxía inglesa, ANalysis of COVAriance). 1.1. O modelo ANOVA O modelo de análise da varianza ou ANOVA é semellante ao modelo de regresión linear (2), pero difire deste dado que considera unha variable explicativa discreta ou categórica. 1.1.1. ANOVA como modelo linear xeral Supóñase que se ten unha variable continua Ymedida en Jgrupos, isto é, Jmostras independentes onde cada mostra está formada por variables independentes entre si e con idéntica distribución N(µj, σ2). Denótase por µjá media da mostra da variable resposta relativa ao grupo j-ésimo, para cada j∈ {1, ... , J}. Ademais, o índice iserá empregado para denotar o valor da variable resposta Yna i-ésima observación para o grupo j. Así disporase dunha mostra {Yij} con i= 1, ... , njpara cada grupo j= 1, ... , J. Denótase tamén por n=PJ j=1 njao número de observacións totais da mostra. Nótese que por supoñerse todas as varianzas σ2iguais, o modelo de análise da varianza é homocedástico por definición. 1 2 1. ANOVA e ANCOVA O modelo ANOVA pode escribirse entón como un modelo de regresión do seguinte xeito: Yij =µj+εij,con i= 1, ... njej= 1, ... , J;(1.1) onde os εij ∈N(0, σ2)son independentes para cada j= 1, ... , J e para cada i= 1, ... , nj. Como xa se intúe en (1.1), o modelo ANOVA é un modelo linear xeral xa que pode ser expresado en notación matricial da forma Y=Xβ+ϵ, mediante as seguintes expresións: Y=                   Y11 . . . Yn11 . . . . . . Y1J . . . YnJJ                  n×1 ,X=                         1 0 ··· 0 . . . . . ..... . . 1 0 ··· 0 0 1 ··· 0 . . .. . ..... . . 0 1 ··· 0 . . .. . ..... . . 0 0 ··· 1 . . .. . ..... . . 0 0 ··· 1                        n×J , β =µ=       µ1 µ2 . . . µJ       J×1 ,ϵ=                   ε11 . . . εn11 . . . . . . ε1J . . . εnJJ                  n×1 ;(1.2) de onde se deduce que E(Y)=(µ1,(n1) ... , µ1, µ2,(n2) ... , µ2, ... , µJ,(nJ) ... , µJ)e que a matriz de covarianzas do modelo é V ar(Y) = V ar(Xβ+ϵ) = V ar(ϵ) = σ2In, con Indenotando a matriz identidade de dimensión n×n. 1.1.2. Estimación dos parámetros Os parámetros a estimar son as medias de cada grupo j= 1, ... , J; é dicir, as compoñentes do vector µdefinido en (1.2). Notación 1.1.Denotarase por Y•j=Pnj i=1 Yij á suma dos valores da variable resposta para o grupo j. Esta notación sustitúe un índice efectuando a adición de todas as observacións obtidas ao recorrer ese índice e fixando os restantes. Empregarase esta notación ao longo do traballo. Nótese que, deste mesmo xeito, Y•j=1 njPnj i=1 Yij =Y•j/njdenotaría a media da mostra dos datos do grupo j,Y•• =PJ j=1 Pnj i=1 Yij denotaría a suma de todos os valores da variable resposta eY•• =PJ j=1 nj nY•jsería a súa media, así pois, a media na mostraxe de todos os datos. Ao igual que para calquera modelo de regresión linear xeral, nos modelos de regresión con erro normal, os métodos de mínimos cadrados ou de máxima verosimilitude proporcionan o mesmo estimador para cada µj, daquela é indiferente que método se empregue para estimar os coeficientes do modelo ANOVA. 1.1. O modelo ANOVA 3 Conforme ao criterio de mínimos cadrados, a expresión a minimizar é: Q= J X j=1 nj X i=1 (Yij −µj)2= n1 X i=1 (Yi1−µ1)2+ n2 X i=1 (Yi2−µ2)2+··· + nJ X i=1 (YiJ −µJ)2. Como cada parámetro só aparece nun sumatorio, para minimizar Qpoderíase pensar en minimizar cada sumatorio por separado. Así, denotando Qj=Pnj i=1 (Yij −µj)2; dQj dµj =−2 nj X i=1 (Yij −µj)=0⇒ nj X i=1 Yij =njµj⇒bµj=1 nj nj X i=1 Yij =Y•j. Deste xeito, o valor estimado para a media poboacional de cada grupo resultará ser a media da mostra de valores Yij asociados ao grupo j, como se podía intuír. Para ver que efectivamente os métodos de máxima verosimilitude e de mínimos cadrados proporcionan os mesmos estimadores para todos os parámetros µj, con j= 1, ... , J; basta con percatarse de que, dado que Yij ∈N(µj, σ2), a función de verosimilitude resulta ser L(µ1, ... , µJ;σ2) = 1 (2πσ2)n/2·exp  −1 2 J X j=1 nj X i=1 Yij −µj σ2 , e efectivamente, maximizar Lcon respecto aos parámetros µjé equivalente a minimizar Qou exp(Q). Deste xeito, os residuos do modelo ANOVA, que seguen a ser a diferenza entre os valores observados e os axustados, resultan ser bεij =Yij −b Yij =Yij −Y•j. Así, os residuos representan a desviación dunha observación con respecto á media estimada do seu grupo. Unha propiedade importante é que os residuos bεij suman cero para cada grupo j, posto que Pnj i=1 bεij = Pnj i=1 (Yij −Y•j) = Pnj i=1 Yij −njY•j= 0 para cada j= 1, ... , J. Exemplo 1.2. Para ilustrar o modelo ANOVA, seguirase empregando a base de datos presentada no Exemplo 0.1. Construirase un modelo en que explique a nota acadada en Matemáticas (variable resposta) en función dunha única variable explicativa, que ademais sexa discreta. Unha situación interesante sería considerar unha variable auxiliar que permita determinar a pertenza ou non do alumnado considerado ás únicas 7 escolas cuxo 100 % do alumnado participou no estudo presentado en [4], coa finalidade de evitar nesgos de selección. Tales escolas denotaranse por E1, E2, E3, E4, E5, E6 e E7. A base de datos creada mediante a consideración de unicamente estas sete escolas denominarase mates7 e será a que se empregue de aquí en diante. A continuación amósase o código elaborado para definir o modelo, xunto co resumo do mesmo obtido coa función summary de . 4 1. ANOVA e ANCOVA anova_mates7 <- lm(NotaMates ~Escola -1) summary(anova_mates7) ## ## Call: ## lm(formula = NotaMates ~ Escola - 1) ## ## Residuals: ## Min 1Q Median 3Q Max ## -17.6957 -3.2007 0.7324 3.6733 10.1293 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## EscolaE1 18.1116 0.7533 24.04 <2e-16 *** ## EscolaE2 16.9639 1.0904 15.56 <2e-16 *** ## EscolaE3 19.7156 0.7138 27.62 <2e-16 *** ## EscolaE4 12.3102 0.8570 14.37 <2e-16 *** ## EscolaE5 18.4557 0.6619 27.89 <2e-16 *** ## EscolaE6 16.2323 0.7620 21.30 <2e-16 *** ## EscolaE7 14.8637 0.6505 22.85 <2e-16 *** ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## ## Residual standard error: 4.997 on 300 degrees of freedom ## Multiple R-squared: 0.9219,Adjusted R-squared: 0.9201 ## F-statistic: 506.1 on 7 and 300 DF, p-value: < 2.2e-16 No apartado concernente á estimación dos coeficientes proporciónase, xunto a cada escola, a estimación bµjda media das notas de Matemáticas acadadas nese mesmo centro educativo “Ej” para un j entre 1e7; que como se viu, non é máis que a media da mostra das notas do correspondente colexio, Y•j. Por outra banda, o nivel crítico asociado a cada unha destas estimacións non ten ningunha utilidade; pois é o p-valor asociado ao contraste de que tales estimacións son nulas, o que non ten sentido tratándose de variables enteiramente positivas (salvo datos illados). Na Figura 1.1 móstrase o comportamento da nota acadada en Matemáticas para as diferentes escolas consideradas a través de diagramas de caixas. Ademais, engádese a media da variable resposta en cada grupo (liña punteada dentro do diagrama de caixa) e a media global (liña 1.1. O modelo ANOVA 5 descontinua azul). Unha primeira ollada amosa que as medias (liña punteada presente en cada caixa) concernentes á nota en Matemáticas quizáis difiran entre as diferentes escolas máis do que deberían se a causa fose soamente o azar. 0 5 10 15 20 25 E1 E2 E3 E4 E5 E6 E7 Escolas Notas en Matemáticas Figura 1.1: Diagramas de caixas da nota de Matemáticas nas 7 escolas con participación absoluta. A liña horizontal descontinua representa a media da mostra global, mentres que en cada caixa a raia punteada é a media da mostra de cada escola. 1.1.3. Análise da varianza e test F Dado o modelo (1.1) con medias grupais µjonde j= 1, ... , J; poderíase pensar se estas medias son todas iguais ou hai algunha que difire do resto. O caso particular de dous grupos, J= 2, non é máis ca un contraste de comparación de medias en poboacións normais, é dicir, o ben coñecido t-test1. Estudaranse máis detalladamente os casos con J≥3. Para contrastar medias de máis de dous grupos de datos, poderíase pensar en facer os contrastes por pares, pero esta estratexia pode ser perigosa. Se temos varios grupos de datos e facemos varias comparacións; é probable que, co tempo, atopemos unha diferenza só por azar. En lugar disto, deberíase aplicar un test integral, e é aquí onde xorde o denominado test F. 1Dada unha variable X∈N(µX, σ2 X)é coñecido que X−µX ScX /√n∈Tn−1, onde ScX é a cuasivarianza da mostra X1, ... , Xn, e tal pivote permitirá realizar intervalos de confianza e contrastes de hipóteses para a media poboacional µX. Así, considerando a variable X−Y∈N(µX−µY, σ2), o pivote convértese en X−Y Sc/√n∈Tn−1; onde S2 cé a cuasivarianza da mostra X1−Y1, ... , Xn−Yn, e permite, en particular, contrastar a diferenza das medias poboacionais. Este tipo de contrastes coñécense como t-test posto que a distribución do pivote é unha t de Student. 6 1. ANOVA e ANCOVA OANOVA emprega un test de hipóteses simples para comprobar se as medias de varios grupos son iguais; isto é, se µ1=··· =µJ, ou se algunha das medias é distinta. O contraste a estudar é o seguinte:    H0:µ1=··· =µJ, Ha:polo menos unha das medias difire do resto. (1.3) Como cada εij ten distribución normal, tamén Yij ∈N(µj, σ2); dado que Yij é unha función linear de εij. Ademais, a independencia dos εij implica a independencia dos Yij. Nótese que, á vista das hipóteses do modelo ANOVA e baixo a hipótese nula H0, estaríase a dicir que a distribución da variable resposta é a mesma nos diferentes grupos considerados. Ademais, a variabilidade total das observacións Yij sen ter en conta os distintos grupos (baixo H0) vén dada polas desviacións de cada valor observado respecto da media global Yij −Y••. Pola contra, se temos en conta os grupos (baixo a hipótese alternativa Ha), a variabilidade debida ao erro vén dada polas desviacións de cada valor observado respecto da media estimada do seu respectivo grupo: Yij −Y•j(= bεij). A diferenza entre ambas expresións coincide coa diferenza entre a media estimada de cada grupo e a media global, isto é, (Yij −Y••)−(Yij −Y•j) = Y•j−Y••. Así, pódese descompoñer a variabilidade total da variable resposta en dúas compoñentes: Yij −Y•• = (Y•j−Y••)+(Yij −Y•j), e en consecuencia, efectuando o cadrado da expresión anterior teríase que J X j=1 nj X i=1 (Yij −Y••)2= J X j=1 nj X i=1 (Y•j−Y••)2+ J X j=1 nj X i=1 (Yij −Y•j)2+ + J X j=1 nj X i=1 2(Y•j−Y••)(Yij −Y•j). Agora ben, J X j=1 nj X i=1 2(Y•j−Y••)(Yij −Y•j) = 2 J X j=1 (Y•j−Y••)· nj X i=1 (Yij −Y•j)!= 0; xa que Pnj i=1 (Yij −Y•j) = Pnj i=1 Yij −njY•j=Y•j−Y•j= 0, e entón dedúcese que J X j=1 nj X i=1 (Yij −Y••)2 | {z } variabilidade total = J X j=1 nj(Y•j−Y••)2 | {z } variabilidade entre grupos + J X j=1 nj X i=1 (Yij −Y•j)2 | {z } variabilidade dentro dos grupos .(1.4) 1.1. O modelo ANOVA 7 O termo da esquerda serve como medida da variabilidade total das observacións Yij e denótase por V T. O primeiro termo na dereita de (1.4) denótase por V EG (variabilidade entre grupos), e o segundo por V DG (variabilidade dentro dos grupos). Deste xeito, (1.4) pode escribirse de forma equivalente como V T =V EG +V DG. Correspondente a esta descomposición da suma total de cadrados, tamén podemos obter unha descomposición dos graos de liberdade asociados a cada sumando: V T ten n−1graos de liberdade asociados. Hai en total nobservacións Yij −Y••, pero un grao de liberdade pérdese porque as desviacións non son independentes no sentido de que suman cero: PJ j=1 Pnj i=1 (Yij −Y••) = Y•• −JnjY•• =Y•• −Y•• = 0. V EG ten J−1graos de liberdade asociados. Hai Jdesviacións da media estimada Y•j−Y•• de cada grupo j, pero un grao de liberdade pérdese igual que na V T, pois PJ j=1 nj(Y•j−Y••) = PJ j=1 njY•j−Y•• PJ j=1 nj=PJ j=1 Y•j−nY •• =Y•• −Y•• = 0 V DG ten n−Jgraos de liberdade asociados. Basta considerar a compoñente de V DG para cada grupo j:Pnj i=1 (Yij −Y•j)2, que é o equivalente á suma total de cadrados considerando só o grupo j, polo que ten nj−1graos de liberdade asociados. Así, os graos de liberdade asociados a V DG son (n1−1) + (n2−1) + ··· + (nJ−1) = n−J. Nótese que J X j=1 nj X i=1 (Yij −Y•j)2= J X j=1 (nj−1)(Yij −Y•j)2 nj−1= J X j=1 (nj−1)S2 j; onde S2 jdenota a cuasivarianza de Yno grupo j. Recompílase esta información na Táboa 1.1, que posteriormente verase que é un bosquexo da verdadeira saída de ao realizar o contraste (1.3). Graos liberdade Sumas de cadrados Grupo J−1V EG =PJ j=1 nj(Y•j−Y••)2 Residuos n−J V DG =PJ j=1 (nj−1)S2 j Total n−1V T =PJ j=1 Pnj i=1 (Yij −Y••)2 Táboa 1.1: Táboa de descomposición da variabilidade asociada a un modelo ANOVA. A análise da varianza ANOVA céntrase en comparar as medias de cada grupo a través da análise da varianza entre grupos e dentro de cada grupo. Esto é, ANOVA considera simultaneamente moitos grupos e evalúa se as súas medias na mostraxe difiren máis do que se esperaría da 8 1. ANOVA e ANCOVA variación natural. Esta variabilidade é a media cuadrática entre grupos e que se denotará por MEG. Por outra banda, para conseguir un valor de referencia sobre canta variabilidade debe esperarse entre as medias da mostra emprégase unha estimación da varianza dentro dos grupos, o que se chamará erro cuadrático medio ou media cuadrática dentro dos grupos MDG. As medias cuadráticas obtéñense dividindo cada suma de cadrados polos seus graos de liberdade asociados: MEG =V EG J−1eMDG =V DG n−J. Para efectuar o contraste de igualdade entre as medias de todos os grupos simultaneamente, empregamos o test F, cuxo estatístico Fadoptará a forma F=MEG MDG =V EG/(J−1) V DG/(n−J), que non é máis que o cociente entre a variabilidade entre grupos e a variabilidade dentro dos grupos. Valores grandes de Faportan indicios a favor de Ha, xa que MEG tenderá a exceder MDG cando Haé certa, xa que a variabilidade entre grupos é maior que a variabilidade dentro de cada grupo. É dicir, teríanse os grupos ”separados“, o cal é un indicativo de que non todas as medias son iguais. Valores pequenos de Faportan indicios a favor de H0, xa que MEG eMDG teñen o mesmo valor esperado baixo H0(ver [20, pp. 538–542]), isto é, baixo a hipótese nula de que todas as medias dos grupos son iguais calquera diferenza entre as medias da mostra débese unicamente ao azar. Para poder construír unha regla de decisión e a rexión crítica desta, necesítase coñecer a distribución do estatístico F. Baixo H0, tense que todas as medias grupais µjson iguais e polo tanto todas as respostas Yij teñen a mesma distribución, e como consecuencia da hipótese de normalidade do modelo ANOVA, aplicando o Teorema de Cochran (pódese consultar en [20, pp. 92–93]) teríase que: Baixo H0,V EG σ2∈χ2 J−1,V DG σ2∈χ2 n−J;e son independentes2. Por tanto, F=MEG MDG =V EG/(J−1) V DG/(n−J)= V EG/σ2 J−1 V DG/σ2 n−J ∈FJ−1,n−J(Baixo H0). Esto é, se H0é certa e as condicións do modelo se verifican, entón o estatístico Fsegue unha distribución F de Snedecor3con graos de liberdade J−1en−J. Ademais, a rexión crítica 2Se Z1, ... , Zmson variables aleatorias normais estándar e independentes; entón X=Z2 1+··· +Z2 m∈χ2 m segue unha distribución chi-cadrado con mgraos de liberdade. É unha distribución non negativa, con media me varianza 2m. 3Se X1∈χ2 m1, X2∈χ2 m2e son independentes, entón F=X1/m1 X2/m2∈Fm1,m2e dise que ten distribución Fde Snedecor con m1em2graos de liberdade. 1.1. O modelo ANOVA 9 do test para un nivel de significación αserá: Rexéitase H0:µ1=··· =µJse F > f1−α;J−1,n−J,(1.5) onde f1−α;J−1,n−Jrepresenta o cuantil de orde 1−αda distribución F de Snedecor con J−1e n−Jgraos de liberdade. Finalmente, para ver como se efectúan os cálculos do modelo de análise da varianza pódese seguir a organización da táboa ANOVA de , que aparece reflexada na Táboa 1.2. Graos liberdade Sumas de cadrados Medias cuadráticas F Grupo J−1V EG =PJ j=1 nj(Y•j−Y••)2MEG =V EG J−1 MEG MDG Residuos n−J V DG =PJ j=1 (nj−1)S2 jMDG =V DG n−J Total n−1V T =PJ j=1 Pnj i=1 (Yij −Y••)2 Táboa 1.2: Táboa de descomposición de variabilidade asociada a un modelo ANOVA obtida grazas ao entorno estatístico . Exemplo 1.3. Aplicando a función anova de ao modelo presentado no Exemplo 1.2, onde recórdese que a variable resposta é a nota acadada en Matemáticas e está medida en 7escolas diferentes, pódese obter a táboa de descomposición da varianza que se amosa na Táboa 1.2. anova(lm(NotaMates ~Escola)) ## Analysis of Variance Table ## ## Response: NotaMates ## Df Sum Sq Mean Sq F value Pr(>F) ## Escola 6 1569.3 261.558 10.475 1.49e-10 *** ## Residuals 300 7490.6 24.969 ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 Obsérvase que o valor do estatístico asociado ao test Fé de 10.475. A partir del é sinxelo calcular o nivel crítico empregando (1.5), que resulta ser 1.49 ×10−10. Agora ben, xa que o nivel crítico ou p-valor é menor que os niveis de significación habituais e, en particular, moito menor que α= 0.001; que será o nivel de significación que se empregue ao longo de todo o traballo, séguese que existen evidencias estatisticamente significativas a favor da hipótese alternativa de que polo menos algunha media é diferente das demais, tal e como se intuía na Figura 1.1. 16 2. Introdución aos modelos mixtos para que un conxunto de datos sexa xerárquico, é dicir, que estea composto por datos multinivel, é que un individuo só pode pertencer a un único grupo do nivel 2, así como a un único grupo de cada un dos niveis posteriores. Do mesmo xeito, un grupo do nivel 2só pode pertencer a un único grupo no nivel 3, así como a un único grupo de cada un dos niveis posteriores ao seu; e o mesmo ten que ocorrer para todos os niveis presentes no conxunto de datos en consideración. Exemplo 2.1. Os datos multinivel son bastante frecuentes no eido educativo, de feito, no Capítulo 1 xa se viu unha estrutura de dous niveis na que as/os alumnas/os compoñen o nivel 1e á súa vez están aniñados en escolas, as cales constitúen o nivel 2; aínda que tamén se poden dar estruturas de máis de dous niveis, segundo [15], en Educación é común atoparse con estruturas de 5niveis (estudante, clase, escola, distrito e área xeográfica). No conxunto de datos mates7 xa empregado anteriormente téñense dispoñibles as variables nota en Matemáticas (que se denota por NotaMates), status socio-económico (que se denota por StSE), pertenza a un grupo racial minoritario (que se denota por SocialMin) e sexo da/o alumna/o (que se denota por Sexo); relativas ao alumnado e que serían polo tanto as variables asociadas ao nivel 1, mentres que as variables media dos status socio-económico do alumnado de cada escola (que se denota por MStSE), número de estudantes de cada escola (que se denota por Tamaño), carácter público ou privado de cada escola (que se denota por Sector), proporción de estudantes de cada escola que participan no estudo [4] (que se denota por Particip), medida do ambiente discriminatorio de cada escola (que se denota por AmbDiscrim) e posesión de máis do 40 % das/os matriculadas/os en cada escola que sexan procedentes de grupos raciais minoritarios (que se denota por MaioriaMin); serían as asociadas ao nivel 2, debido a que estas non varían entre todas as observacións que incumben a unha mesma escola. Ademais, como xa se adiantaba ao final da Introdución, as notas acadadas en Matemáticas polas/os estudantes probablemente estén enormemente influenciadas polas diferentes escolas. 2.2. A necesidade de ter en conta os distintos niveis 17 Nesta situación, é importante destacar que a nota acadada na materia de Matemáticas por dúas/dous alumnas/os de distintos institutos que quizáis nin sequera se coñecen é totalmente independente. Agora ben, isto non ocorre con nenas e nenos que acoden diariamente á mesma escola, posto que terán bastante en común; dende o seu estilo de vida até a organización das clases de Matemáticas; pode ser incluso que lles imparta a materia o mesmo docente. Por conseguinte, a asunción de independencia feita nos modelos clásicos estudados ao longo do Grao non é asumible nesta situación; e resulta naturalmente necesario construír novos modelos que teñan en conta a característica inherente que é a xerarquía do conxunto de datos que se está a tratar. 2.2. A necesidade de ter en conta os distintos niveis Cando os individuos forman grupos ou clusters, é obvio pensar que o máis probable é que individuos clasificados nun mesmo grupo resulten ter un comportamento máis semellante que uns individuos calquesquera de grupos diferentes e polo tanto con menos información en común. Tal e como afirma Goldstein [13], incluso cando os grupos se asignan aleatoriamente (no peor dos casos) a miúdo tende a haber diferenzas entre estes. Unha asunción típica na Estatística clásica, en particular no modelo de regresión dun só nivel presentado en (2), é que as observacións son independentes e identicamente distribuídas. Agora ben, se o conxunto de datos está formado por datos multinivel e aínda que non se teña en conta o efecto dos grupos á hora de construír un modelo de regresión, entón a hipótese de independencia non se verifica. Un xeito de ter en conta os efectos grupais sería incluír no modelo variables dummy2como variables explicativas, tal e como se fixo previamente para os modelos ANOVA eANCOVA (modelos de efectos fixos) en (1.1) e (1.6), respectivamente; os cales poderían ser apropiados se o interese principal fose facer Inferencia sobre precisamente os grupos presentes na mostra. Así e todo xorden impedimentos cando o número de grupos é grande, xa que o número de parámetros a estimar medra considerablemente e o modelo de regresión pode non ser o suficientemente eficiente. Ademais, se o interesante non son precisamente eses grupos, senón que se consideran como unha mostra (aleatoria) dunha poboación máis grande e o que se quere é facer Inferencia sobre todos os grupos en xeral; por exemplo, se en lugar das 7escolas o realmente importante é sacar conclusións arredor de tódalas escolas estadounidenses, entón os modelos de efectos fixos ANOVA eANCOVA deixan de ser válidos. É aquí onde xorden os modelos mixtos ou modelos multinivel, en particular; os modelos mixtos con variable resposta continua que se abordarán nas seguintes seccións. Ademais, de aquí en diante, consideraranse unicamente datos xerárquicos cunha estrutura de dous niveis, onde os individuos se atopan no nivel 1e están aniñados en grupos no nivel 2. As observacións 2As variables dummy son unhas variables ficticias que serven para indicar a posible pertenza dos individuos a cada grupo, tomando o valor 1no caso de que un individuo pertenza a un determinado grupo, e o 0en caso contrario. 18 2. Introdución aos modelos mixtos entre niveis ou grupos distintos serán independentes, mentres que as observacións dentro dun mesmo grupo resultarán dependentes entre si posto que pertencen á mesma subpoboación. Por conseguinte, falarase de dúas fontes de variación:entre grupos edentro dos grupos ou intra-grupos. Para o axuste con de modelos mixtos empregarase o paquete lme4 (acrónimo de Linear Mixed-Effects Models using ’Eigen’ and S4) que pode consultarse en Bates et al. [3], que é unha versión máis moderna do antigo nlme que contiña o estudo de [4] de onde se extraeu a información para crear o conxunto de datos mates7. O paquete lme4 emprega métodos de álxebra linear máis eficientes (os do paquete Eigen), ademais de ser máis rápido computacionalmente e empregar menos memoria. Resultará de especial interese a función lmer, que serve para axustar modelos mixtos lineares. No seguinte capítulo, comezarase explicando o modelo mixto máis sinxelo (o coñecido como modelo RANOVA), que non é máis que a extensión natural do modelo ANOVA; e deseguido engadiranse ao modelo covariables medidas no nivel 1(Capítulo 4) e concernentes ao nivel 2 (Capítulo 5); asemade irase ilustrando o comportamento destes novos métodos na práctica empregando a base de datos mates7 a prol de construír modelos máis sofisticados que permitan explicar mellor as notas acadadas polo alumnado na materia de Matemáticas e coñecer que parte das notas é debida ao propio alumnado e cal foi debido ao “efecto” da escola, no cal van implicitamente incluídas as modalidades de ensino, a formación do profesorado, etc. No Anexo A reproducirase parte do código de empregado. Lémbrese que a totalidade do código de atoparase ademais no repositorio Modelos_Mixtos_con_R de GitHub e será de carácter enteiramente reproducible, tanto o relativo á construción dos modelos que se estudarán a continuación, como o atinente a todas as figuras ilustradas ao longo da totalidade deste TFG. Capítulo 3 RANOVA Tal e como se adiantaba no capítulo anterior, as seguintes páxinas adicaranse a explicar a natureza do modelo mixto máis sinxelo de todos, que resulta ser a extensión natural do modelo ANOVA visto no Capítulo 1. De aquí en diante, seguirase supoñendo que se dispón dunha mostra aleatoria de datos multinivel, Yij, cunha estrutura de dous niveis, onde o segundo nivel está constituído por Jgrupos e o nivel 1confórmano njindividuos de cada grupo j, con j= 1, ... , J. Xa se reparou en que existen ocasións nas que os grupos non teñen un interés intrínseco, senón que constitúen unha mostra (aleatoria) dun conxunto de moitos máis grupos e o que se quere é facer Inferencia sobre todos os grupos. Por exemplo, en canto á base de datos mates7, o que realmente interesa non son eses 7colexios, senón toda a poboación de escolas dos Estados Unidos de América. Nestas circunstancias, o ANOVA con efectos fixos presentado en (1.1) deixa de ser relevante e é necesario extendelo a un modelo que incorpore efectos aleatorios, isto é, un modelo no cal as medias dos grupos deixen de ser constantes para converterse en variables aleatorias que seguen unha distribución normal cuxa media é a media xeral e cuxa varianza determina a capacidade de influenza de cada grupo. 3.1. Análise da varianza con efectos aleatorios Neste contexto xorde o modelo RANOVA ou modelo de análise da varianza con efectos aleatorios (coñecido na literatura anglosaxoa como Random effects ANOVA), que non é máis que o modelo mixto ou multinivel máis sinxelo de todos e que pode escribirse como Yij =µ+uj+εij,con i= 1, ... njej= 1, ... , J;(3.1) onde µé a media global, os uj∈N(0, σ2 u)son independentes e identicamente distribuídos, os εij ∈N(0, σ2 ε)son independentes, e tamén ujeεij son variables aleatorias independentes entre si 19 20 3. RANOVA para cada j= 1, ... , J ei= 1, ... , nj. Con σ2 ueσ2 εdenotamos as varianzas entre grupos e dentro dos grupos, respectivamente. O modelo (3.1) é similar ao modelo ANOVA presentado en (1.1). A maior distinción é que no ANOVA as medias dos grupos µjson constantes, mentres que no RANOVA tense que µj= µ+uj∈N(µ, σ2 u), con j= 1, ... , J, son variables aleatorias. Por ese motivo o modelo (3.1) é chamado modelo de análise da varianza con efectos aleatorios. Precisamente unha das vantaxes dos modelos mixtos é a habilidade de combinar os datos introducindo efectos aleatorios multinivel. No RANOVA pénsase a priori nunha cantidade non fixa e sen límite de grupos, o interesante non son os µ1, ... , µJparticulares do estudo, senón toda a posible poboación de µj, especialmente a media dos µj,µ; e a variabilidade dos µj, medida por σ2 u. Mentres que σ2 ué unha medida directa da variabilidade dos µj, o efecto desta variabilidade sóese medir pola razón V PC =σ2 u σ2 u+σ2 ε =variabilidade entre grupos variabilidade total ,(3.2) denominado coeficiente de partición da varianza (coñecido polas siglas en inglés VPC). Nótese que esta razón toma valores entre 0(cando σ2 u= 0, logo µj=µpara todo j= 1, ... , J e non hai diferenzas entre grupos, co cal (3.1) non é máis que unha regresión linear ordinaria) e 1 (cando σ2 ε= 0 ou Yij =Yjpara todo i, é dicir, non hai diferenzas dentro de cada grupo). Nótese que o denominador de V PC representa a variabilidade da variable resposta Y, pois V ar(Yij) = V ar(uj+εij) = V ar(uj) + V ar(εij)+2Cov(uj, εij) = σ2 u+σ2 ε, xa que V ar(uj) = σ2 u,V ar(εij) = σ2 εe as variables aleatorias ujeεij son independentes entre si. Por conseguinte, o modelo RANOVA verifica a hipótese de homocedasticidade. En vista destas propiedades, a razón VPC mide a proporción da variabilidade total dos Yij que é explicada pola variabilidade dos µj, isto é, a proporción total de varianza atribuíble ás diferenzas entre grupos. Deste xeito, cando o cociente (3.2) toma valores próximos a cero, o efecto das diferenzas entre os grupos na variabilidade total é insignificante, e cando a razón é maior ou igual a 0.5considerarase que unha cantidade considerable da variabilidade total é explicada polas diferenzas entre grupos. Por outra banda, a covarianza entre dous individuos (que se denotarán por iei′) de distintos grupos (que se denotarán por jej′) é nula, isto é: Cov(Yij, Yi′j′) = Cov(uj+εij, uj′+εi′j′) = =Cov(uj, uj′) + Cov(uj, εi′j′) + Cov(εij, uj′) + Cov(εij, εi′j′)=0, por causa de tódalas hipóteses de independencia asumidas previamente; mentres que para dous individuos iei′dun mesmo grupo j, resulta que Cov(Yij, Yi′j) = Cov(uj+εij, uj+εi′j) = =Cov(uj, uj) + Cov(uj, εi′j) + Cov(εij, uj) + Cov(εij, εi′j) = =Cov(uj, uj) = σ2 u. 3.2. Estimación dos parámetros 21 Deste xeito, a correlación entre dous individuos dun mesmo grupo resulta ser p=Cor(Yij, Yi′j) = Cov(Yij, Yi′j) pV ar(Yij)pV ar(Yi′j)=σ2 u σ2 u+σ2 ε , a que se denominará correlación intra-grupos de nivel 2. Ao longo deste traballo, na análise de datos reais denotarémolo por correlación intra-escola1, ao igual que decide facelo Goldstein en [13]. Por conseguinte, en modelos con só dous niveis o VPC coincide co coeficiente de correlación intra-escola p; que proporciona unha medida da homoxeneidade dos inidividuos dentro de cada grupo. En particular, segundo [16], os valores da correlación intra-escola no eido educativo soen estar entre 0.05 e0.20. Finalmente, nótese que isto é certo para modelos de dous niveis e só para estes; é sinxelo darse de conta de que ao engadir outro nivel ao modelo, por exemplo as clases nas que están ubicadas/os as/os estudantes, a correlación pxa non tería a mesma expresión, senón unha notablemente máis complexa. 3.2. Estimación dos parámetros No modelo de análise da varianza con efectos aleatorios, os erros de primeiro e segundo nivel son variables aleatorias N(0, σ2 u)eN(0, σ2 ε), respectivamente; e o seu comportamento queda caracterizado entón polas varianzas σ2 ueσ2 ε. Así, os parámetros a estimar no RANOVA son precisamente µ,σ2 ueσ2 ε. 3.2.1. Estimación da media global µ Considéranse as medias grupais Y•j=µ+uj+ε•j, con j= 1, ... , J. Así, E(Y•j) = E(µ+uj+ε•j) = µ+E(uj) + E(ε•j) = µ, e V ar(Y•j) = V ar(µ+uj+ε•j) = V ar(uj+ε•j) = =V ar(uj) + V ar(ε•j)+2Cov(uj, ε•j) = σ2 u+σ2 ε nj ; onde a última igualdade é certa debido a que as variables aleatorias ujeεij son independentes entre si e mais ao teorema de Fisher2. Deste xeito, a media da mostra de cada grupo é un estimador innesgado da media global; pero non é consistente, no sentido de que a varianza nunca será nula por moito que aumente o número de datos do grupo j. Para solventar o problema e mellorar 1No ámbito da estatística este termo é frecuentemente coñecido por correlación intra-clase (ou ICC, empregando a notación de [9]), pero isto pode resultar confuso ao empregalo no eido educativo. 2Se X= (X1, ... , Xn)é unha mostra aleatoria simple dunha poboación N(µ, σ2), entón X∈N(µ, σ2/n)e nV ar(X)/σ2∈χ2 n−1. A demostración do Teorema de Fisher pódese consultar en [19, p. 66]. 22 3. RANOVA o estimador, denotando V ar(Y•j) = σ2 u+σ2 ε nj= ∆j, basta con considerar unha combinación das medias da mostra, bµ= J X j=1 ∆−1 j PJ h=1 ∆−1 h Y•j;(3.3) onde ∆−1 jé a “precisión” de Y•j. Cando a precisión é igual en todos os grupos, o estimador obtido non é máis que o promedio global, isto é, bµ=∆−1PJ j=1 Y•j J∆−1=1 J J X j=1 Y•j=Y•• se ∆−1 j= ∆−1para todo j= 1, ... , J. Pola contra, se algunha precisión difire das outras; é necesario obter un estimador para V ar(Y•j) = ∆j, e este depende de σ2 ueσ2 ε; cuxas estimacións serán as seguintes en procurar obterse. 3.2.2. Estimación das compoñentes da varianza Analogamente ao feito para o modelo ANOVA en (1.4), considérase a seguinte descomposición da variabilidade da variable de interese Y: J X j=1 nj X i=1 (Yij −Y••)2= J X j=1 nj X i=1 (Yij −Y•j)2+ J X j=1 nj(Y•j−Y••)2;(3.4) e seguindo a mesma notación, (3.4) pode escribirse de xeito equivalente como: V T |{z} variabilidade total =V DG |{z} variabilidade dentro dos grupos +V EG |{z} variabilidade entre grupos . Para estimar a varianza de primeiro nivel σ2 ε, considérase bσ2 ε=V DG N−J, onde segue a denotarse N=PJ j=1 nj. Para a variabilidade entre grupos (ou de segundo nivel), se todos os grupos constan do mesmo número de observacións (datos balanceados), entón podemos considerar S2 u=V EG n(J−1), poñendo n=njpara un j∈ {1, ... , J}calquera. Ocorre que E(bσ2 ε) = σ2 ε, mentres que E(S2 u) = σ2 u+σ2 ε n> σ2 u; en consecuencia, o estimador bσ2 εé innesgado, mentres que non o é S2 u. Para anular o nesgo deste último, basta con considerar bσ2 u=S2 u− bσ2 ε n, establecendo bσ2 u= 0 no caso de que o resultado fose negativo, entendéndoo como que non hai efecto de grupos. A demostración destes feitos atópase deseguido e está baseada nos cálculos que fixo Kutner para o modelo ANOVA con efectos fixos, que poden consultarse en [20, pp. 538–542]. Proposición 3.1. Baixo as hipóteses formuladas para o modelo RANOVA e dados bσ2 ε=V DG N−Je S2 u=V EG n(J−1), verifícase que (i) E(bσ2 ε) = σ2 ε, 3.2. Estimación dos parámetros 23 (ii) E(S2 u) = σ2 u+σ2 ε n. Demostración. (i) Denotando a cuasivarianza das observacións nun grupo jpor S2 j, tense que bσ2 ε=1 N−J J X j=1 nj X i=1 (Yij −Y•j)2=1 N−J J X j=1 (nj−1)Pnj i=1 (Yij −Y•j)2 nj−1= =1 N−J J X j=1 (nj−1)S2 j, e empregando que S2 jé un estimador innesgado da varianza dentro dos grupos (corolario inmediato do Teorema de Fisher), séguese que, efectivamente; E(bσ2 ε) = 1 N−J J X j=1 (nj−1) E(S2 j) = 1 N−J J X j=1 (nj−1)σ2 ε=σ2 ε. (ii) Consideraranse datos balanceados, así pois, nj=npara todo j. Deste xeito, S2 u=V EG n(J−1) =1 J−1 J X j=1 (Y•j−Y••)2. Agora ben,    Y•j=µ+uj+ε•j,con ε•j=1 njPnj i=1 εij,e Y•• =µ+u•+ε••, logo Y•j−Y•• = (uj−u•) + (ε•j−ε••); e elevando ao cadrado e sumando por grupos séguese que J X j=1 (Y•j−Y••)2= J X j=1 (uj−u•)2 | {z } (a) + J X j=1 (ε•j−ε••)2 | {z } (b) + 2 J X j=1 (uj−u•)(ε•j−ε••) | {z } (c) . Debido a que a esperanza é un operador linear, para calcular a esperanza do termo da esquerda pódese calcular a de todos estes sumandos por separado e logo sumalas. (a) Analogamente ao feito en (i), J X j=1 (uj−u•)2= (J−1)PJ j=1 (uj−u•)2 J−1= (J−1)S2 uj; onde S2 ujdenota a cuasivarianza dos uje polo tanto E  J X j=1 (uj−u•)2 = (J−1) E(S2 uj)=(J−1)σ2 u. 24 3. RANOVA (b) Basta darse de conta de que PJ j=1 (ε•j−ε••)2/(J−1) é unha varianza na mostra, e polo tanto un estimador innesgado da varianza da variable ε•j; pero ε•jé precisamente a media de nerros independentes εij; e entón do Teorema de Fisher séguese que V ar(ε•j) = V ar(εij) n=σ2 ε n⇒E  J X j=1 (ε•j−ε••)2 =(J−1)σ2 ε n. (c) Pola independencia entre si das variables aleatorias ujeεij e mais a linearidade da esperanza, pódese escribir E 2 J X j=1 (uj−u•)(ε•j−ε••) = 2 J X j=1 E(uj−u•)E(ε•j−ε••). Agora ben, xa que E(εij) = 0, entón E(ε•j) = 0 eE(ε••) = 0; co cal E(ε•j−ε••) = 0 e a esperanza de (c) é nula. Finalmente, pódese concluír que efectivamente E(S2 u) = 1 J−1E  J X j=1 (Y•j−Y••)2 =σ2 u+σ2 ε n. Para poder establecer o estimador da varianza entre grupos σ2 uinicialmente proposto (S2 u) supúxose que o número de observacións en cada grupo era o mesmo. En caso contrario, a estimación de σ2 urequerirá dunha estimación previa do efecto fixo (media global), que á súa vez precisa das estimacións da varianza. Neste caso, para obter as estimacións dos parámetros do RANOVA é necesario recorrer a procedementos iterativos. Cando se estima un modelo multinivel mediante o método de máxima verosimilitude (ML na literatura inglesa), a dependencia entre os parámetros (efectos fixos e compoñentes da varianza) reflectida previamente para o caso particular do modelo RANOVA leva implícito recorrer a procedementos iterativos; primeiramente estímanse os efectos fixos, por exemplo mediante o algoritmo EM (esperanza-maximización) que considera os efectos aleatorios como faltantes, ou ben empregando o método de mínimos cadrados xeralizados iterativo (IGLS, acrónimo de Iterative Generalised Least Squares), sobre os cales se pode atopar máis información en [28, pp. 41–42] ou en [25, pp. 164–165], respectivamente. Deseguido, empréganse as estimacións dos efectos fixos obtidas nesta primeira iteración para estimar as compoñentes da varianza, e á súa vez estas estimacións das compoñentes da varianza empréganse para obter unhas novas estimacións dos efectos fixos na seguinte iteración. Continuando deste xeito até que non haxa cambios significativos entre iteracións consecutivas obteranse as estimacións procuradas, tanto dos efectos fixos 3.2. Estimación dos parámetros 25 como das compoñentes da varianza. Na Figura 3.1 amósase un bosquexo do que fai cando se emprega o método ML para estimar os parámetros dun modelo mixto. Pola contra, o método de máxima verosimilitude produce estimacións nesgadas dos parámetros aleatorios, posto que estimar os efectos fixos en primeiro lugar implica que toda a variabilidade da mostra debida a estes sexa ignorada, co cal as compoñentes da varianza son subestimadas, incrementándose así tanto os t-valores como os p-valores ou niveis críticos asociados. Isto é realmente preocupante en bases de datos pequenas, en concreto, cun número de grupos reducido, xa que a variabilidade da mostra dos efectos fixos tende a ser máis grande; asunto que se trata máis profundamente en [22], onde tamén se propón unha solución; o método de máxima verosimilitude restrinxida, coñecido polas súas siglas en inglés REML. O método REML, ao contrario que ML, estima os efectos fixos e as compoñentes da varianza por separado. Ambos procedementos son asintoticamente equivalentes e igualmente complexos no relativo ao eido da computación. Para evitar a estimación simultánea que realiza ML, REML non fai máis que restrinxir a 0os efectos fixos nun primeiro momento, permitindo estimar así as compoñentes da varianza separadamente. Finalmente, estímanse os efectos fixos empregando as estimacións das compoñentes da varianza obtidas previamente, mediante unha modificación do método IGLS coñecida por RIGLS (acrónimo de Restricted Iterative Generalised Least Squares) e sen empregar procedementos iterativos. Toda a teoría concernente ao método RIGLS explícaa Goldstein en [12]. Na parte inferior da Figura 3.1 ilústrase o procedemento que se segue cando se emprega o método REML no contexto dos modelos multinivel. Efectos fixos Compoñentes da varianza Efectos fixos Compoñentes da varianza ... · EM ou IGLS aa Compoñentes da varianza Compoñentes da varianza ... Compoñentes da varianza Efectos fixos RIGLS Figura 3.1: Comparación entre os procedementos ML (parte superior) e REML (parte inferior). En definitiva, por unha banda, tal e como afirma [7] convén buscar a sinxeleza dos modelos; mentres que por outra banda, en canto ao procedemento por máxima verosimilitude, estimando os efectos fixos en primeiro lugar, ademais de ignorarse a variabilidade da mostra debida a estes, tampouco se teñen en conta os graos de liberdade que se perden estimándoos. É por iso que ao longo deste documento se decida traballar co método ML para os modelos mixtos máis sinxelos mentres que do Capítulo 4 en diante, ao engadir efectos aleatorios nas pendentes asociadas a covariables relativas ao nivel 1, asemade os graos de liberdade comecen a escasear, se procure empregar máis o método REML. 32 3. RANOVA 3.3. Contraste sobre os efectos grupais Cando o V PC se atopa próximo a cero, así pois, σ2 ué case nulo, parece haber evidencias de que realmente non hai ningunha diferenza verdadeira entre os distintos grupos. Para contrastar esta afirmación pódese empregar un contraste de hipóteses a fin de comprobar se o efecto das diferenzas entre os grupos é nulo ou difire máis do que se podería esperar só por azar. O contraste a estudar é o seguinte:    H0:non hai diferenzas entre os grupos, Ha:hai diferenzas entre os grupos. ⇐⇒    H0:σ2 u= 0, Ha:σ2 u= 0. (3.5) Como uj∈N(0, σ2 u), se a hipótese nula é certa entón necesariamente todos os ujson nulos. Logo, o contraste anterior pode pensarse como un contraste entre dous modelos diferentes, que ademais están aniñados6;   H0:Yij =µ+εij, Ha:Yij =µ+uj+εij. (3.6) en ambos casos con i= 1, ... njej= 1, ... , J. En canto ao modelo baixo a hipótese nula, coñécese habitualmente como modelo dun só nivel para a media, posto que tamén se podería expresar como Yi=µ+εicon i= 1, ... , n; é dicir, sen ter en conta os grupos e reflectindo simplemente a desviación de cada dato con respecto á media global. Tales diferenzas entre o valor da variable resposta Ye a media global veñen dadas por εi, que non é máis que o erro para cada individuo ino modelo; e estes seguen unha distribución normal de media cero e unha varianza σ2, isto é, εi∈N(0, σ2)para cada dato i= 1, ... , n. A varianza σ2simboliza a variabilidade ao redor da media; deste xeito, se fose cero, todos os puntos terían o mesmo valor. Do mesmo xeito, canto máis grande é a varianza, máis grandes son as desviacións sobre a media. Na Figura 3.4 ilústrase a natureza do modelo dun só nivel para a media, en concreto, para o caso particular das sete escolas. É evidente a diferenza imperante entre ámbolos dous modelos do contraste (3.6). Mentres que na Figura 3.2 concernente ao RANOVA os residuos ao nivel alumnado resultaban ser a distancia de cada dato á media do seu respectivo grupo; no modelo baixo a hipótese nula, ao non reparar no efecto das escolas, os residuos (de cor negra) non son máis que as distancias de cada dato á media global µ, tal e como se amosa na Figura 3.4. 6Dous modelos dinse aniñados se o modelo máis complexo se pode construír a partir do máis simple engadindo un ou máis parámetros. Por exemplo, considérense tres modelos M1,M2eM3.M1considera como variable explicativa a X1,M2considera a X1e a X2eM3considera a X1e a X3.M1eM2están aniñados, así como tamén M1eM3; pero non o están M2eM3. 3.3. Contraste sobre os efectos grupais 33 ε17,1 ε20,2 ε42,5 ε24,6 Y17,1 Y20,2 Y42,5 Y24,6 µ 0 5 10 15 20 25 Nota en Matemáticas Escola E1 E2 E3 E4 E5 E6 E7 Figura 3.4: Modelo dun só nivel para a media (global), sen considerar o efecto das escolas, xunto con algúns residuos (de cor negra) e cos seus respectivos datos da cor correspondente á escola. Agora ben, xa que ámbolos dous modelos están aniñados, á hora de construír un estatístico para realizar o contraste (3.5), podemos comparar os dous modelos de (3.6) mediante un test de razón de verosimilitudes, isto é, empregando o estatístico LR =−2 log L1 L2=−2 log L1−(−2 log L2)∈χ2 1, onde L1é a verosimilitude7para o modelo baixo a hipótese nula en (3.6), igualmente L2é a verosimilitude para o modelo de análise da varianza con efectos aleatorios e “log” refírese ao logaritmo natural. O estatístico LR segue unha chi cadrado cun único grao de liberdade, xa que o modelo máis complexo soamente incorpora un parámetro adicional con respecto ao máis sinxelo, que precisamente é a varianza entre grupos σ2 u. A demostración deste feito omítese neste traballo pero pode atoparse en [26, pp. 417–419]. Rexeitar a hipótese nula implica que hai evidencias estatisticamente significativas a favor da existencia de diferenzas entre os distintos grupos, en cuxo caso un modelo multinivel sería máis axeitado que un modelo que non tivese en conta os grupos. Agora ben, se cun determinado conxunto de datos non existen evidencias en contra da hipótese nula e entón non se pode rexeitar, isto non quere dicir que non se deban ter en conta os grupos á hora de axustar un modelo para 7Considérese un vector aleatorio Xcon función de densidade fθ, sendo θ∈Rqun vector de parámetros descoñecido. Entón, dada X1, ... , Xnunha mostra aleatoria simple de X, pódese estimar θempregando a función de verosimilitude; que vén dada por L(θ) = Qn i=1 f(xi(θ)) e de onde se obtén o estimador facendo b θ= m´axθ∈RqL(θ). 34 3. RANOVA eses datos; é posible que as diferenzas entre os diferentes grupos se revelen soamente logo de engadir máis variables explicativas ao modelo e que permanezan agochadas ao considerar tan só as diferenzas entre os grupos relativas á media. Exemplo 3.3. Realizarase o contraste sobre os efectos das distintas escolas para ver se efectivamente a escola á que asista un/unha estudante ten influenza sobre a nota que acade en Matemáticas. Para axustar o modelo dun só nivel para a media das notas en Matemáticas baixo aH0en (3.6) abonda con considerar un modelo de regresión para a media da escola dos xa vistos na Introdución mediante a función lm de : vcmates <- lm(NotaMates ~1,data = mates7) Tras o axuste deste modelo, obtense que a estimación da media global µnon é outra que a media dos valores almacenados na variable NotaMates, isto é, bµ= 16.823. O valor do estatístico para o test de razón de verosimilitudes e o p-valor asociado ao contraste calcúlanse deseguido: LR <- -2*logLik(vcmates)[1]-(-2*logLik(ranovamates))[1] 1-pchisq(LR, df =1) Na Figura 3.5 represéntase a función de densidade χ2 1que segue o estatístico LR, cuxo valor vén dado pola liña vermella e, como se pode apreciar, deixa unha probabilidade moi baixa á súa dereita (1.96 ×10−9). Así pois, existen evidencias estatisticamente significativas a favor da hipótese alternativa de (3.5) de que hai diferenzas entre as distintas escolas en canto á nota acadada polas/os súas/seus estudantes na materia de Matemáticas. χ1 2 LRobs 1.96e−09 0.00 0.05 0.10 0.15 0 10 20 30 40 Figura 3.5: Función de densidade dunha chi-cadrado cun grao de liberdade en laranxa xunto co valor do estatístico do contraste sobre os efectos das escolas LR (de cor vermella) e mais o p-valor asociado a tal contraste. Capítulo 4 Modelos mixtos con covariables relativas ao primeiro nivel De xeito análogo ao que se fixo no Capítulo 1 introducindo os modelos ANOVA eANCOVA, considerarase agora unha extensión do modelo de análise da varianza con efectos aleatorios ou RANOVA engadindo unha ou máis variables continuas que estén relacionadas coa variable resposta Y. Consideraranse neste capítulo modelos que inclúan soamente unha variable explicativa relativa ao primeiro nivel, o cal será suficiente para explicar con detalle a natureza destes. Os datos desta covariable denotaranse por Xij, con j= 1, ... , J ei= 1, ... , nj; e o efecto desta pode ser ou ben fixo (modelo con intercepto aleatorio e pendente fixa) ou ben aleatorio (modelo con intercepto e pendente aleatorias). A continuación describiranse este tipo de modelos. 4.1. Modelo con intercepto aleatorio Como se viu anteriormente, un modelo de análise da varianza con efectos aleatorios pode escribirse como Yij =µ+uj+εij,con j= 1, ... , J ei= 1, ... , nj. Pois ben, para ter en conta a información sobre os individuos que proporciona unha certa variable Xe considerando fixo o efecto desta, basta con engadila ao modelo RANOVA multiplicada por un coeficiente fixo, do xeito: Yij =β0j+β1Xij +εij,con j= 1, ... , J ei= 1, ... , nj;(4.1) onde β0j=µ+u0j∈N(µ, σ2 u0)é un intercepto aleatorio que varía entre os distintos grupos; mentres que a pendente β1é a mesma para todos os grupos. No eido dos modelos multinivel, a este modelo coñéceselle por modelo con intercepto aleatorio (e pendente fixa). De maneira 35 36 4. Modelos mixtos con covariables relativas ao primeiro nivel análoga ao presentado no Capítulo 3, µrepresenta a media global, os u0j∈N(0, σ2 u0)son independentes e identicamente distribuídos, os εij ∈N(0, σ2 ε)son independentes, e tamén u0j eεij son variables aleatorias independentes entre si para cada j= 1, ... , J ei= 1, ... , nj. Destas propiedades séguese ademais que os interceptos aleatorios β0json independentes entre si e independentes dos erros de primeiro nivel. Como interpretación xeométrica, a representación gráfica do axuste dun modelo multinivel con intercepto aleatorio sería bastante similar á mostrada na Figura 3.2 elaborada para o RANOVA, só que agora tanto a recta de regresión axustada para todos os datos sen ter en conta os grupos (a media global) como as rectas de regresión axustadas para os distintos grupos terán todas unha pendente β1; como consecuencia de considerar o efecto fixo da covariable X. É dicir, as liñas da forma b Yij =bµ+bu0j+b β1Xij, con j= 1, ... , J;conformarán polo tanto un conxunto de rectas paralelas. Na Figura 4.1 ilústranse estes feitos, amosando a recta de regresión xeral xunto coas rectas estimadas para os colexios E3, E4, E5 e E7. Recórdese que o modelo de análise da covarianza con efectos fixos ou ANCOVA, explicado na Sección 1.2, vén dado por Yij =µ+τj+γXij +εij,con j= 1, ... , J ei= 1, ... , nj; e, polo tanto, fixado un grupo, non é máis ca unha recta de regresión con intercepto (µ+τj) e pendente γ(a mesma para todos os grupos). Este modelo é bastante similar ao modelo con intercepto aleatorio que se está a considerar pero con interceptos fixos en lugar de aleatorios. Por tanto, aínda que se construíu o modelo con intercepto aleatorio como unha extensión natural do RANOVA, tamén se podería elaborar a partir do modelo de análise da covarianza con efectos fixos ou ANCOVA. Deste xeito, o modelo con intercepto aleatorio non é máis ca un modelo ANCOVA (sen interacción) con efectos aleatorios, é dicir, un RANCOVA1sen interacción. O modelo con intercepto aleatorio e pendente fixa verifica a hipótese de homocedasticidade, xa que V ar(Yij|X=Xij) = V ar(µ+u0j+β1Xij +εij|X=Xij) = V ar(u0j+εij) = =V ar(u0j) + V ar(εij)+2Cov(u0j, εij) = σ2 u0+σ2 ε, por ser as variables aleatorias u0jeεij independentes entre si. Ademais, a covarianza entre 1Non é común o emprego desta denominación na literatura sobre os modelos multinivel, só se emprega neste caso particular para clarificar que o modelo con intercepto aleatorio non é máis ca un ANCOVA con efectos aleatorios (Random effects ANCOVA). 4.1. Modelo con intercepto aleatorio 37 observacións de distintos grupos sempre é nula2: Cov(Yij, Yi′j′|Xij, Xi′j′) = Cov(µ+u0j+β1Xij +εij, µ +u0j′+β1Xi′j′+εi′j′|Xij, Xi′j′) = =Cov(u0j+εij, u0j′+εi′j′) = Cov(u0j, u0j′) + Cov(u0j, εi′j′)+ +Cov(εij, u0j′) + Cov(εij, εi′j′)=0, por ser as variables aleatorias u0jeεij independentes separadamente. Por outra banda, entre observacións (iei′) dun mesmo grupo jresulta ser Cov(Yij, Yi′j|Xij, Xi′j) = Cov(µ+u0j+β1Xij +εij, µ +u0j+β1Xi′j+εi′j|Xij, Xi′j) = =Cov(u0j+εij, u0j+εi′j) = Cov(u0j, u0j) + Cov(u0j, εi′j)+ +Cov(εij, u0j) + Cov(εij, εi′j) = V ar(u0j) = σ2 u0. Exemplo 4.1. Nas seguintes liñas, extenderase o modelo RANOVA construído previamente mediante a introdución dunha variable explicativa concernente ao nivel 1do alumnado, o status socio-económico (que se vén denotando por StSE). Xa se viu na Sección 1.2 relativa á construción do modelo ANCOVA que o efecto do StSE nas notas acadadas por un/unha alumna/o era significativo. Na Figura 4.1 represéntase o axuste do modelo obtido coa seguinte sintaxe en : matesmm1 <- lmer(NotaMates ~StSE +(1|Escola), REML =FALSE) A pendente b β1= 1.397 é efectivamente a mesma para todas as escolas, mentres que o intercepto é diferente para cada colexio; o que ocasiona rectas paralelas. A media global estímase por bµ= 16.156, que non é máis que o intercepto da recta negra asociada á escola “media” (bu0j= 0), entendéndoo como a nota estimada cando StSE = 0. Así, para obter o intercepto asociado á escola E3, por exemplo, non hai máis que sumarlle a bµa predición do efecto aleatorio para tal colexio bu3obtido mediante a función ranef de ; resultando 16.156 + 3.037 = 19.193, valor que predí a recta descontinua azul cando StSE = 0. Por outra banda, as estimacións da varianza dos coeficientes aleatorios resultan ser bσ2 u0= 4.845 ebσ2 ε= 24.168. Reparando nas disimilitudes deste modelo e do RANOVA advírtese que a adición do StSE reduciu a varianza a nivel de alumnado e mais a varianza total, o cal era esperado porque StSE é unha variable concernente ao nivel das/os estudantes. A varianza a nivel escolar, pola contra, incrementouse lixeiramente. Isto é debido a que o status socio-económico non se distribúe de forma moi regular entre as escolas. Aínda que en todos os centros hai alumnas e alumnos máis desfavorecidas/os e máis beneficiadas/os, non todos os colexios se ubican en lugares 2Denótase Cov(Yij, Yi′j′|Xij , Xi′j′)por comodidade o que realmente son covarianzas condicionais: Cov(Yij, Yi′j′|X=Xij , X =Xi′j′), e o mesmo acontece con todos os termos sucesivos. 38 4. Modelos mixtos con covariables relativas ao primeiro nivel ε4,3 ε32,4 ε4,5 ε51,7 Y4,3 Y32,4 Y4,5 Y51,7 β1 u3 u4 u5 u7 0 5 10 15 20 25 −1 0 1 Status socio−económico Nota en Matemáticas Escola E1 E2 E3 E4 E5 E6 E7 Figura 4.1: Ilustración do modelo mixto con intercepto aleatorio e pendente fixa para a nota de Matemáticas nas 7escolas con participación unánime, xunto coa recta de regresión para a escola media, E3, E4, E5, E7 e mais algúns residuos concernentes a ambos niveis (alumnado e escolas). coas mesmas características sociais ou económicas nin tampouco son frecuentados polo mesmo tipo de alumnado. Pode atoparse máis información arredor deste fenómeno que ocorre sobre as varianzas ao engadir variables explicativas de nivel 1aos modelos nos apuntamentos do curso LEMMA impartido polo Centro de Modelos Multinivel da Universidade de Bristol [30]. A Figura 4.1 é un claro exemplo deste fenómeno, xa que os puntos de determinadas cores tenden a situarse máis cara a parte esquerda do gráfico mentres que outros se espallan máis cara a dereita, até os valores máis altos do StSE. Tal é o caso, por exemplo, das escolas E4 e E5. O alumnado do centro educativo E4 (de cor rosa) ten, en promedio, un status socio-económico claramente menor que o da escola E5 (de cor verde). De feito, a media do StSE para o colexio E4 resulta ser 0.367, mentres que para o centro E5 é 0.759. Outra fonte de variación entre escolas pode verse tamén na Figura 4.1, posto que as rectas correspontes a centros educativos con maior intercepto (liñas da parte superior), tenden a ter máis alumnas/os con StSE maiores que 0.5; mentres que as escolas con interceptos máis pequenos tenden a ter máis alumnas/os cun menor status socio-económico. Logo de ter en conta os efectos do StSE, a proporción total de varianza atribuíble ás diferenzas entre as distintas escolas sería: V PC =bσ2 u0 bσ2 u0+bσ2 ε =4.845 4.845 + 24.168 = 0.167 = 16.7 %, fronte ao 15.84 % sen considerar o efecto do StSE. Tal incremento é debido ao acrecentamento 4.1. Modelo con intercepto aleatorio 39 da varianza a nivel dos centros educativos. Finalmente, nótese que no caso de ser a covariable Xunha variable dicotómica, µsería a media total de Ypara os individuos con x= 0,µ+u0jsería a media para os individuos con x= 0 no grupo j, e a “pendente” β1sería a diferenza na media que posúen os datos con x= 1 relativa á dos datos con x= 0 (en calquer grupo). Ilústrase a continuación esta situación. Exemplo 4.2. Co propósito de ver que ocorre ao considerar unha variable explicativa dicotómica e de mellorar o modelo construído no Exemplo 4.1, engádeselle a este o efecto fixo asociado á variable SocialMin, que lémbrese indicaba se a/o estudante é membro dun grupo racial minoritario ou non. Codificarase o valor “Non” por 0e o “Si” por 1(cando a/o menor pertence a un grupo racial minoritario). Así, o modelo a considerar é o seguinte, NotaMatesij =µ+β1StSEij +β2SocialMinij | {z } efectos fixos +uj |{z} efectos aleatorios +εij. A continuación, amósase un resumo do modelo axustado. matesmm2 <- lmer(NotaMates ~StSE +SocialMin +(1|Escola), REML =FALSE) summary(matesmm2) ## Linear mixed model fit by maximum likelihood ['lmerMod'] ## Formula: NotaMates ~ StSE + SocialMin + (1 | Escola) ## ## AIC BIC logLik deviance df.resid ## 1861.4 1880.0 -925.7 1851.4 302 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -3.8713 -0.6260 0.1161 0.7342 2.1362 ## ## Random effects: ## Groups Name Variance Std.Dev. ## Escola (Intercept) 4.35 2.086 ## Residual 23.17 4.814 ## Number of obs: 307, groups: Escola, 7 ## ## Fixed effects: ## Estimate Std. Error t value ## (Intercept) 16.9501 0.8817 19.223 40 4. Modelos mixtos con covariables relativas ao primeiro nivel ## StSE 0.9937 0.4509 2.204 ## SocialMinSi -2.7527 0.7455 -3.692 ## ## Correlation of Fixed Effects: ## (Intr) StSE ## StSE -0.240 ## SocialMinSi -0.244 0.241 Para un/unha mesma/o estudante, acrecentar o StSE exactamente nun punto suporía un incremento na nota de b β1= 0.99 puntos; independentemente de se a/o alumna/o procede dun grupo racial minoritario ou non. De xeito semellante, para estudantes co mesmo StSE; o simple feito de pertencer a un grupo racial minoritario (SocialMin = 1) automaticamente augura unha mingua na nota acadada de |b β2|= 2.75 puntos. Por outra banda, a interpretación do intercepto fixo bµ= 16.95 debería ser a do prognóstico para un/unha alumno/a calquera cun status socioeconómico igual a 0, non procedente dun grupo racial minoritario e sen considerar o efecto da escola á que acode. Ademais, todos estes coeficientes resultan ser significativos segundo o criterio de [27], xa que os t-valores resultan todos maiores que 1.96 en valor absoluto. Na Figura 4.2 represéntanse de cor negra as rectas de regresión xerais para a escola media (bu0j= 0), unha relativa aos datos sobre alumnado minoritario e outra diferente para o non minoritario. Posto que b β2=−2.753 <0, a liña relativa ás/aos menores procedentes de grupos raciais minoritarios é a que se atopa por debaixo; e o mesmo ocorre para cada par de liñas de regresión axustadas relativo a un determinado centro, en particular, para os pares de liñas asociadas ás escolas E4 e E5; que son as dúas que a modo de exemplo se representan na gráfica. 4.2. Modelo con intercepto e pendente aleatorios Extenderase agora o modelo con intercepto aleatorio visto previamente mediante a incorporación de efectos aleatorios na pendente asociada á variable explicativa X, o que consistirá en formular un modelo linear en cada grupo j= 1, ... , J e ademais pode representarse como Yij =β0j+β1jXij +εij,con j= 1, ... , J ei= 1, ... , nj;(4.2) onde β0j=γ00 +u0jeβ1j=γ10 +u1json variables aleatorias independentes separadamente (e independentes dos erros), con distribución normal bivariante; β0j β1j!∈N γ00 γ10!,Σu= σ2 u0σu01 σu01 σ2 u1!!, 4.2. Modelo con intercepto e pendente aleatorios 41 ε33,4 ε26,5 Y33,4 Y26,5 β1 β1 β2 β2 β2 u4 u5 u4 u5 0 5 10 15 20 25 −1 0 1 Status socio−económico Nota en Matemáticas Escola E1 E2 E3 E4 E5 E6 E7 Racial minori− tario Non Si Figura 4.2: Ilustración do modelo mixto con intercepto aleatorio para a nota de Matemáticas nas 7escolas considerando como variables explicativas StSE eSocialMin, xunto coas rectas de regresión para a escola media, E4, E5 e mais algúns residuos concernentes a ambos niveis. onde γ00 eγ10 son o intercepto medio e a pendente media, respectivamente; σ2 u0é a varianza do intercepto, σ2 u1é a varianza da pendente e σu01 é a covarianza entre intercepto e pendente. Seguiranse a supor as mesmas hipóteses sobre a independencia que ata agora; as variables aleatorias u0j,u1jeεij son todas independentes e os efectos aleatorios do intercepto e da pendente son independentes dos erros de primeiro nivel. A diferenza reside en que agora hai variables non independentes entre si, posto que dentro de cada grupo pode existir correlación entre os efectos aleatorios do intercepto e da pendente. Analogamente ao símil do modelo con intercepto aleatorio e un suposto RANCOVA sen interacción, poderíase aludir a este modelo como un RANCOVA con interacción, agora con efectos aleatorios tanto no intercepto como na pendente. Precisamente por variar o intercepto e a pendente entre os distintos grupos é polo que (4.2) se denomina modelo con intercepto e pendente aleatorios. Ao igual que acontecía no Exemplo 4.1 coa adición dos efectos fixos dunha variable concernente ao nivel 1, a introdución de efectos aleatorios nunha variable a nivel de individuo reducirá a varianza nese nivel (σ2 ε, varianza dentro dos grupos); mentres que a varianza entre grupos (σ2 u0 eσ2 u1) pode manterse e mesmo aumentar. Considerando que as fontes de variabilidade que interveñen no modelo se poden estruturar en 48 5. Modelos mixtos con covariables relativas ao segundo nivel A variable explicativa de nivel 2pode engadirse a un modelo multinivel exactamente do mesmo xeito que unha variable explicativa de nivel 1. Así, se Wjé a variable contextual e Xij é a variable asociada ao nivel 1; o modelo con intercepto aleatorio (4.1) convértese en Yij =β0j+β1Xij +β2Wj+εij,con j= 1, ... , J ei= 1, ... , nj;(5.1) onde de novo β0j=γ00 +u0j∈N(γ00, σ2 u0)eεij ∈N(0, σ2 ε), e ademais estas variables aleatorias son independentes tanto por separado como entre si mesmas, para cada j= 1, ... , J ei= 1, ... , nj. 5.1. Variables contextuais composicionais No caso particular de que a variable contextual sexa a media dunha variable Xdo primeiro nivel, entón a variable composicional non é outra que Wj=X•jpara cada j= 1, ... , J. Así, (5.1) reescríbese do seguinte xeito, Yij =β0j+β1Xij +β2X•j+εij,con j= 1, ... , J ei= 1, ... , nj;(5.2) onde β1é o efecto de dentro dos grupos de X,β1+β2é o efecto entre grupos e β2é o efecto contextual, é dicir, o efecto das medias grupais de X(X•j) sobre Yque non se manifestaba considerando só os efectos máis elementais. Co obxectivo de que o efecto entre grupos β1+β2 sexa un coeficiente máis do modelo, a prol de que mediante a implementación informática se estime directamente, o modelo (5.2) pode expresarse equivalentemente como Yij =β∗ 0j+β∗ 1(Xij −X•j) + β∗ 2X•j+εij = =γ00 +γ10(Xij −X•j) + δX•j | {z } efectos fixos +u0j |{z} efectos aleatorios +εij;(5.3) para cada j= 1, ... , J ei= 1, ... , nj; sen máis que facer β∗ 0j=β0j=γ00 +u0j,β∗ 1=β1=γ10 eβ∗ 2=β1+β2. O efecto contextual é frecuente denotalo por δ=β1+β2. Isto tamén se pode pensar considerando efectos aleatorios na pendente, sen máis que ter en conta o termo β∗ 1j=β1j=γ10 +u1jen lugar de β∗ 1. Realmente, o que se fai é considerar como variables explicativas as diferenzas entre as observacións e as medias dos seus grupos (Xij −X•j), no que na literatura concernente aos modelos multinivel se soe denominar un modelo de Cronbach, posto que a idea de centrar a variable de primeiro nivel antes de introducila no modelo co propósito de distinguir o efecto entre grupos e o efecto dentro dos grupos foi inicialmente proposta polo experto en Psicoloxía da Educación homónimo ao modelo; segundo afirman Hox et al. en [18]. Xa que tanto o modelo presentado en (5.2) como o modelo de Cronbach (5.3) son equivalentes, ao longo deste TFG farase uso do primeiro, que é o que permite unha sintaxe máis sinxela e eficiente da función lmer de . En consecuencia, a estimación do efecto entre grupos (b β1+b β2) haberá que calculala por separado. 5.1. Variables contextuais composicionais 49 Exemplo 5.1. Sobre a mesma base de datos coa que se viña traballando (mates7), no Exemplo 4.3 considerábase o modelo que explicaba a nota en Matemáticas en función do status socioeconómico (variable que se vén denotando por StSE) con efectos aleatorios, e mais do efecto fixo concernente á pertenza a un grupo racial minoritario (variable categórica dada por SocialMin). Agora engadiráselle ao mesmo a variable composicional MStSE, que non é máis que a media dos status socio-económico das/os estudantes para cada escola e que polo tanto verifica MStSEj= StSE•jpara cada j= 1, ... , J. Xa que o seu valor é o mesmo para todas/os as/os menores dun mesmo centro educativo, a variable MStSE é unha variable contextual, e por medir o status socioeconómico promedio dos colexios séguese que é ademais unha variable composicional2. O modelo construído virá dado pola expresión NotaMatesij =γ00 +γ10StSEij +β2SocialMinij +β3MStSEj | {z } efectos fixos +u0j+u1jStSEij | {z } efectos aleatorios +εij;(5.4) de novo para todo j= 1, ... , J ei= 1, ... , nj. Así, a única diferenza coa expresión do previo modelo con intercepto e pendente aleatorios (4.3) estudado no Exemplo 4.3 é a adición do termo β3MStSEj, co cal a diferenza coa Figura 4.3 do modelo (4.3) consistirá puramente nos interceptos das distintas rectas. Transcribindo esta fórmula na sintaxe requerida pola función lmer de , obtense o seguinte axuste do modelo (5.4): matesmm4 <- lmer(NotaMates ~StSE +SocialMin +MStSE +(1+StSE |Escola), REML =TRUE) summary(matesmm4) ## Linear mixed model fit by REML ['lmerMod'] ## Formula: NotaMates ~ StSE + SocialMin + MStSE + (1 + StSE | Escola) ## ## REML criterion at convergence: 1843.9 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -3.8761 -0.6031 0.1107 0.7389 2.1732 ## ## Random effects: ## Groups Name Variance Std.Dev. Corr 2Algúns autores non concordan en que as variables composicionais sexan un caso particular das contextuais, e tratan estas dúas situacións conxuntamente sen facer ningún tipo de distinción (como [28]). Neste traballo, pola contra, deféndese esta clasificación; ao igual que deciden facelo moitos outros autores salientábeis na literatura dos modelos multinivel (como [18], [10] ou [21]). Á hora de analizar modelos mixtos non é importante que clasificación se esté a empregar, é soamente unha cuestión conceptual; esclarece a que nivel pertence unha certa variable. 50 5. Modelos mixtos con covariables relativas ao segundo nivel ## Escola (Intercept) 6.69176 2.5868 ## StSE 0.05674 0.2382 -1.00 ## Residual 23.30565 4.8276 ## Number of obs: 307, groups: Escola, 7 ## ## Fixed effects: ## Estimate Std. Error t value ## (Intercept) 17.0893 1.6662 10.256 ## StSE 1.0275 0.4670 2.200 ## SocialMinSi -2.7128 0.7537 -3.599 ## MStSE -0.4411 3.3888 -0.130 ## ## Correlation of Fixed Effects: ## (Intr) StSE SclMnS ## StSE -0.124 ## SocialMinSi -0.102 0.241 ## MStSE -0.773 -0.150 -0.037 ## optimizer (nloptwrap) convergence code: 0 (OK) ## boundary (singular) fit: see help('isSingular') O coeficiente fixo asociado á variable MStSE estímase por b β3=−0.441, co cal un neno ou nena, procedente dun grupo racial minoritario ou non e cun status socio-económico fixado de antemán; ao incrementarse nun punto o promedio dos status socio-económico das/os demais alumnas/os da súa escola, espérase que experimente unha diminución na súa nota de Matemáticas de case medio punto. En definitiva, se un membro do alumnado se queda por detrás dos seus compañeiros e compañeiras; ben socialmente ou ben economicamente, non só ficará inamovible no eido da Educación mentres o resto perfecciona as súas notas, senón que incluso sufrirá unha mingua nos seus resultados; en particular, na materia de Matemáticas. Na Figura 5.1 ilústrase a natureza do modelo (5.4). Como xa se adiantaba, é moi semellante á Figura 4.3; a diferencia reside enteiramente nos interceptos asociados ás distintas escolas. Mentres que no Exemplo 4.3 as estimacións dos efectos aleatorios do intercepto bu0jcon j= 1, ... , J se interpretaban como a desviación do intercepto para a escola j-ésima con respecto ao intercepto medio (o da recta de cor negra, entendéndoo como o intercepto habitual cando StSE = 0); neste contexto é preciso ter en conta a estimación do efecto fixo asociado á variable contextual composicional MStSE á hora de estudar os interceptos dos distintos colexios, posto que para obter o intercepto dunha escola jtamén haberá que sumarlle agora a bu0jo termo b β3MStSEj=b β3StSE•j. Na Figura 5.1 tamén se representa (de cor rosa e punteada) a modo de exemplo a liña constituída 5.1. Variables contextuais composicionais 51 polas predicións das notas en Matemáticas para as/os alumnas/os da escola E4 procedentes dun grupo racial minoritario, é dicir, con SocialMin = 1. Xuntando o dito anteriormente co que se fai no Capítulo 4, o intercepto desta recta (entendéndoo de novo cando StSE = 0) calcúlase sumando ao intercepto da escola media o termo bu0j+b β3StSE•j+b β2, polo que o intercepto concernente á liña dunha escola calquera será |b β2|= 2.71 unidades menor polo simple feito de considerar só o alumnado procedente dun grupo racial minoritario. ε7,3 ε26,4 Y7,3 Y26,4 u 03 +β 3 ×StSE • 3 u 04 +β 3 ×StSE • 4 γ10 γ10 +u14 γ10 +u14 γ 10 +u 13 β2 0 5 10 15 20 25 −1 0 1 Status socio−económico Nota en Matemáticas Escola E1 E2 E3 E4 E5 E6 E7 Racial minori− tario Non Si Figura 5.1: Ilustración do modelo mixto con intercepto e pendente aleatorios para a nota de Matemáticas nas 7escolas considerando pendente fixa asociada á variable categórica SocialMin; así como o efecto da variable composicional MStSE, xunto coa recta de regresión para a escola “media” cando SocialMin = 0 e para as escolas E3 e E4. Por outra banda, empregando o criterio xa explicado de [27]; como o t-valor asociado á estimación do coeficiente β3é de -0.13 (moito menor que 1.96 en valor absoluto), séguese que tal coeficiente non resulta significativo. Máis formalmente, empregando un test de razón de verosimilitudes onde a hipótese nula é o modelo (4.3) e a alternativa é o actual (5.4); obtense un nivel crítico de 0.04, maior que o nivel de significación α= 0.001 que se viña empregando ata o momento. Ademais, segundo [27, p. 295] ou [9, p. 175] hai moitos indicios a favor de que este p-valor é demasiado pequeno, é dicir, de que se subestimou e realmente é maior de 0.043. En consecuencia, como o nivel crítico é maior que o nivel de significación considerado, non tería senso rexeitar a hipótese nula, isto é, o promedio dos status socio-económico do alumnado dunha 3Para obter p-valores asociados aos tests de razón de verosimilitudes máis precisos pode empregarse o denominado procedemento bootstrap paramétrico, sobre o cal se pode atopar máis información en [27, pp. 294–301]. 52 5. Modelos mixtos con covariables relativas ao segundo nivel escola non inflúe de xeito significativo na nota dun/dunha estudante calquera. Conseguintemente, a estimación das varianzas tamén aumenta con respecto ás obtidas no modelo (4.3) axustado no Exemplo 4.3 (excepto a concernente á pendente aleatoria asociada á variable StSE,bσu1= 0.057). 5.2. Variables contextuais globais O seguinte paso natural no incremento da complexidade do modelo (5.1) é permitir que a pendente asociada a Xvaríe entre os diferentes grupos relativos ao segundo nivel; co cal o modelo con intercepto e pendente aleatorios tórnase agora en Yij =β0j+β1jXij +β2Wj+εij,con j= 1, ... , J ei= 1, ... , nj;(5.5) onde se seguirán a supor as mesmas hipóteses sobre a independencia que até agora, é dicir, de novo εij ∈N(0, σ2 ε)son independentes e ademais β0j=γ00 +u0j∈N(γ00, σ2 u0)e β1j=γ10 +u1j∈N(γ10, σ2 u1), con j= 1, ... , J, son variables aleatorias independentes separadamente (e independentes dos erros). Así, dentro de cada grupo pode existir correlación entre os efectos aleatorios do intercepto e da pendente, sendo σu01 a súa covarianza. De aquí en diante, o coeficiente fixo asociado á variable global Wescribirase tamén β2=γ01. Nótese que o modelo (5.5) pode escribirse como: Yij =γ00 +γ10Xij +γ01Wj | {z } efectos fixos +u0j+u1jXij | {z } efectos aleatorios +εij. Exemplo 5.2. Partindo do modelo (4.3) axustado na Sección 4.2, xa que o modelo axustado previamente e que lémbrese incluía como covariable a variable contextual composicional MStSE non resultou ser significativo; engadiránselle agora variables contextuais de tipo global a prol de construír un novo modelo que axuste as notas en Matemáticas de xeito máis eficiente. Na base de datos mates7 disponse de 5variables globais; aínda que as variables Sector (carácter público ou privado da escola) e Particip (proporción de estudantes da escola que participan no estudo) non resultan útiles, posto que as 7escolas que se están a ter en consideración son todas católicas e tiveron unha participación unánime no estudo [4]. Pódense empregar para axustar o novo modelo, polo tanto, as variables Tamaño,AmbDiscrim eMaioriaMin; que lémbrese reflectían o número de estudantes, unha medida do ambiente discriminatorio e se a cantidade das/os matriculadas/os membros dun grupo racial minoritario supera o 40 % do total da escola, respectivamente. Nótese que a variable MaioriaMin, tal e como se describe, ben podería tratarse dunha variable composicional do mesmo estilo que a do Exemplo 5.1. Isto é debido a que MaioriaMin podería construírse para as sete escolas da base de datos mates7 simplemente calculando proporcións sobre a variable de primeiro nivel SocialMin (que indicaba se un/unha estudante era membro dun grupo racial minoritario). Pola contra, isto só pode facerse no caso particular das escolas E1, 5.2. Variables contextuais globais 53 E2, ... , E7 debido a que tiveron unha participación unánime no estudo presentado en [4]; mais con todos os outros centros educativos presentes na base de datos orixinal isto non é posible, posto que non se ten información acerca de todas/os as/os estudantes; indicio de que tal variable se extraeu nun primeiro momento dunha fonte externa ao estudo académico descrito en [4]. Polo tanto, neste traballo considerarase que MaioriaMin é unha variable contextual global. De novo saliéntase que á hora de axustar un modelo multinivel, a clasificación que se empregue non ten importancia algunha, polo que se pola contra MaioriaMin se considerase como unha variable composicional; as conclusións que se acadarían non variarían no máis mínimo coas obtidas a continuación. Deseguido resumiranse as conclusións obtidas grazas ao entorno estatístico . O código empregado pode atoparse no Anexo A.8.2. Por unha banda, considerando o efecto fixo do Tamaño ao engadilo ao modelo (4.3) obtense que a estimación do coeficiente asociado a este resulta ser case nulo. Ademais, tras a realización dun test de razón de verosimilitudes compróbase que non resulta ser significativo. Conclúese entón que o tamaño dunha escola non inflúe de xeito significativo na nota acadada por un/unha estudante na materia de Matemáticas. Unha explicación a dito resultado podería ser a dada por [14] en relación ao tamaño das clases. Ao contrario do que se soe pensar, que as clases sexan máis numerosas non sempre implica un menor rendemento académico. Isto ocorre tan só cando o número de estudantes con malas notas nesa mesma clase supera un certo umbral; e no caso particular dos 7centros de mates7, as porcentaxes de alumnas/os con notas menores que 5puntos (por exemplo) son, respectivamente: ## Escola ## E1 E2 E3 E4 E5 E6 E7 ## 0.00 4.76 0.00 8.82 0.00 0.00 6.78 Ningunha porcentaxe supera o umbral do 10 % (por exemplo), isto podería ser a xustificación de que o tamaño das escolas non resulte significativo. Por outra banda, engadindo agora o efecto da variable global AmbDiscrim ao modelo (4.3), tal e como era de esperar a estimación do coeficiente resulta negativa; isto é, o aumento do ambiente discriminatorio nunha certa escola ten efectos perniciosos sobre as notas das/os súas/seus estudantes. Pola contra, novamente este coeficiente tampouco resulta significativo. A variable que resta por engadir ao modelo (4.3) é MaioriaMin, variable dicotómica que se codifica por 1para as escolas con máis do 40 % do alumnado membro dun grupo racial minoritario, e por 0no caso contrario. No Exemplo 4.2 viuse que o simple feito de pertencer a un grupo racial minoritario auguraba unha mingua importante na nota de Matemáticas. Pola contra, ao axustar un novo modelo tendo en conta a variable MaioriaMin, para as escolas con máis do 40 % do alumnado procedente dun grupo racial minoritario augúrase un incremento do intercepto 54 5. Modelos mixtos con covariables relativas ao segundo nivel de bγ01 = 1.175 puntos na nota de Matemáticas (entendéndoo cando StSE = 0). Isto non goza de senso ningún e pode deberse á estrutura dos datos; posto que só a escola E2 (a de menor tamaño) verifica que MaioriaMin = 1, e estes datos poderían resultar insuficientes para extraer conclusións axeitadas. De feito, engadir a variable MaioriaMin resulta máis significativo que a anexión das dúas anteriores; aínda que non o suficiente, pois o test de razón de verosimilitudes entre o modelo considerando MaioriaMin e o previo (4.3) fornece un p-valor de 0.094; maior que o nivel de significación α= 0.001 considerado neste traballo; e ademais lémbrese que segundo [27] ou [9], o p-valor de 0.094 foi subestimado e o real aínda podería ser aínda máis grande. 5.3. Interacción entre niveis Afondando aínda máis no terreo dos modelos multinivel, poderíase permitir que o efecto dunha covariable dependa do valor doutra variable explicativa. No caso particular de que unha variable estea asociada ao nivel 1e a outra ao nivel 2isto denomínase interacción entre niveis. Deste xeito, partindo do previo modelo (5.5) con intercepto e pendente aleatorios no que se introduce información sobre algunha característica do grupo a través dunha variable de segundo nivel; e considerando interacción entre niveis, estase a construír un modelo cuxas formulacións xerárquica e combinada son, para todo i= 1, ... , nje para todo j= 1, ... , J: (i) Formulación xerárquica: Nivel 1 : Yij =β0j+β1jXij +εij. Nivel 2 :    β0j=γ00 +γ01Wj+u0j,con u0j∈N(0, σ2 u0), β1j=γ10 +γ11Wj+u1j,con u1j∈N(0, σ2 u1). (ii) Formulación conxunta: Yij =γ00 +γ10Xij +γ01Wj+γ11WjXij | {z } efectos fixos +u0j+u1jXij | {z } efectos aleatorios +εij. De novo se supón que os erros son independentes dos efectos aleatorios do intercepto e pendente, pero estes últimos poden estar correlados. Séguense impoñendo tamén as hipóteses de normalidade sobre os erros de ambos niveis. Nótese que a parte aleatoria coincide coa do modelo con intercepto e pendente aleatorios (4.2), sen considerar a variable de segundo nivel W. A principal característica deste modelo e asemade a diferenza co modelo sen interacción reside na inserción do termo WjXij, que se pode entender como a interacción entre variables de ambos niveis; sendo γ11 o coeficiente de tal interacción. Finalmente, no caso de ser a covariable Wunha variable dicotómica, a interpretación dos coeficientes é semellante á feita no Exemplo 4.2 para SocialMin. Capítulo 6 Conclusións Ao longo deste traballo viuse que os modelos de regresión clásicos, como o modelo de regresión linear ou os modelos de análise da varianza e covarianza non son suficientes cando se está a traballar con datos aniñados xerarquicamente, os cales son moi frecuentes no eido da Educación, a Medicina ou as Ciencias Medioambientais. Cando os individuos forman grupos, é obvio pensar que os individuos clasificados nun mesmo grupo tenderán a ter un comportamento máis semellante que uns individuos calesquera de grupos diferentes e polo tanto con menos información en común. Ademais, habitualmente o interesante non son os grupos presentes no conxunto de datos, senón que se pretende aplicar técnicas da Inferencia Estatística sobre unha poboación máis grande de grupos. É aquí onde xorden os denominados modelos mixtos, modelos multinivel ou modelos de efectos aleatorios. En particular, neste traballo tratáronse os modelos mixtos lineares con resposta continua e púxose de manifesto a súa utilidade para estudar bases de datos cunha estrutura xerárquica de dous niveis, onde os individuos se atopan no primeiro nivel e están aniñados en grupos no segundo nivel, mediante a incorporación de efectos aleatorios. Cabe salientar que os modelos mixtos lineares vistos neste traballo poderían ampliarse a formas máis complexas, como os modelos mixtos lineares xeralizados (denotados habitualmente como GLMM ), que xorden cando a variable resposta Yposee certas restricións, como por exemplo, que se trate dunha variable de reconto. Estes modelos son a continuación natural dos modelos mixtos lineares, e toda a información concernente aos modelos GLMM poden consultarse na extensa obra de Demidenko [6]. Outra posible extensión dos modelos mixtos lineares sería asumir que en cada grupo non se ten unha recta, senón un modelo polinómico, por exemplo. Deste xeito, daríase máis flexibilidade para modelar o comportamento da variable resposta en cada grupo de interese. Ao longo deste traballo ilustrouse a utilidade dos modelos mixtos lineares con resposta continua no eido da Educación mediante a base de datos mates7 creada para esta ocasión e mais 55 56 6. Conclusións a ferramenta estatística . Lémbrese que o código desenvolto ao longo deste traballo está dispoñible no repositorio Modelos_Mixtos_con_R de GitHub e é de carácter enteiramente reproducible. Isto serviu para estudar que características e situacións inflúen na nota de Matemáticas dun/dunha alumno/a calquera, asemade se construíu un modelo mixto para tratar de predecir as notas das/os mesmas/os. Obtívose que a escola á que acode un/unha alumno/a inflúe enormemente na nota acadada na materia de Matemáticas por este/a, así como a situación socioeconómica da súa familia. Ademais, estudouse como inflúe a integridade social das/os estudantes na nota dende un punto de vista racial; obtendo que as/os alumnas/os procedentes dun grupo racial minoritario acadan en promedio unha nota significativamente menor que as/os outras/os estudantes en xeral. Por outra banda, tamén se tocaron aspectos concernentes á política no eido da Educación. Ao contrario do que se soe crer, tanto o tamaño das escolas como o número de alumnas/os por clase non inflúen na nota das/os estudantes. Isto só ocorre cando a porcentaxe de menores de baixo rendemento na mesma clase ou na mesma escola supera certo umbral. É por isto que os diferentes gobernos do Mundo non deberían invertir os seus bens en diminuír o tamaño dos centros educativos ou o ratio de alumnas/os por clase, senón en mellorar a calidade da Educación que se ofrece nas/os mesmas/os. Anexo A Código de R Neste Anexo reproducirase parte do código de empregado ao longo de todo o traballo. A súa integridade pode consultarse no repositorio Modelos_Mixtos_con_R de GitHub, tanto o relativo á construción dos diversos modelos que se estudaron como o concernente a todas as figuras ilustradas, todas programadas empregando o paquete ggplot2. A sintaxe empregada á hora de redactar o código é a defendida por Hadley Wickham en [31]. Tamén merecen especial recoñecemento [32] e [5]; posto que moitas das seguintes liñas de código puideron elaborarse grazas ás ideas proporcionadas nestas obras. Por último, o código disponse seguindo o formato que emprega o paquete knitr, o cal tamén se empregou para xerar este documento e sobre o cal se pode atopar máis información na obra orixinal do autor do paquete Yihui [33]. Todos estes libros son de código aberto e atópanse dispoñibles en GitHub. A.1. Creación da base de datos mates pkgs <- c('lme4','nlme','lmtest','ggplot2','car','dplyr','RColorBrewer', 'viridis','gtools','patchwork','knitr','ggthemes','lattice','latex2exp') install.packages(setdiff(pkgs, installed.packages()[,"Package"]), dependencies =TRUE) library(lme4); library(nlme); library(lmtest); library(ggplot2); library(car) library(dplyr); library(RColorBrewer); library(viridis); library(gtools) library(patchwork); library(knitr); library(ggthemes); library(lattice) library(latex2exp) mates <- merge(MathAchieve, MathAchSchool, by ="School") 57 64 A. Código de R coefs_anc1 <- ancova1$coefficients #estimacións ANCOVA sen interacción dark_2 <- brewer.pal(n=7,name ="Dark2")#paleta de cores rang <- data.frame(i=numeric(7), f=numeric(7)) for (k in 1:length(levels(Escola))){ rang[k ,] <- range(StSE[Escola == levels(Escola)[k]]) }#Función para os rangos das rectas recta <- function(x=0,grupo =1){ if (!grupo %in% 1:7){ stop("Tal grupo non existe.") } if (grupo == 1){ as.numeric(coefs_anc1[1]) +as.numeric(coefs_anc1[2]) *x } else{ as.numeric(coefs_anc1[1]) +as.numeric(coefs_anc1[grupo +1]) + as.numeric(coefs_ancova1[2]) *x } }#Función para as ecuacións das rectas ancova_si <- ggplot(mates7, aes(x= StSE, y= NotaMates, color = Escola)) + geom_point() + labs(x="Status socio-económico",y="Notas en Matemáticas", color ="Escola")+ coord_cartesian(ylim =c(-.5,25) ) + scale_color_brewer(palette ="Dark2")+ geom_segment(x= rang[1,1], xend = rang[1,2], y=recta(rang[1,1]), yend = recta(rang[1,2]), col = dark_2[1], lwd =1.25)+ geom_segment(x= rang[2,1], xend = rang[2,2], y=recta(rang[2,1], 2), yend = recta(rang[2,2], 2), col = dark_2[2], lwd =1.25)+ geom_segment(x= rang[3,1], xend = rang[3,2], y=recta(rang[3,1], 3), yend = recta(rang[3,2], 3), col = dark_2[3], lwd =1.25)+ geom_segment(x= rang[4,1], xend = rang[4,2], y=recta(rang[4,1], 4), yend = recta(rang[4,2], 4), col = dark_2[4], lwd =1.25)+ geom_segment(x= rang[5,1], xend = rang[5,2], y=recta(rang[5,1], 5), yend = recta(rang[5,2], 5), col = dark_2[5], lwd =1.25)+ geom_segment(x= rang[6,1], xend = rang[6,2], y=recta(rang[6,1], 6), yend = A.6. RANOVA 65 recta(rang[6,2], 6), col = dark_2[6], lwd =1.25)+ geom_segment(x= rang[7,1], xend = rang[7,2], y=recta(rang[7,1], 7), yend = recta(rang[7,2], 7), col = dark_2[7], lwd =1.25)+ theme( legend.position ="none", axis.title.x =element_text(size =16), axis.text.x =element_text(size =14), axis.title.y =element_text(size =16), axis.text.y =element_text(size =14)) ancova_ci <- ggplot(mates7, aes(x= StSE, y= NotaMates, color = Escola)) + geom_point() + labs(x="Status socio-económico",color ="Escola")+ coord_cartesian(ylim =c(-.5,25) ) + scale_color_brewer(palette ="Dark2")+ geom_smooth(method ="lm",se =FALSE,lwd =1.25)+ theme( legend.title =element_text(size =16), legend.text =element_text(size =14), axis.title.x =element_text(size =16), axis.text.x =element_text(size =14), axis.title.y =element_blank(), axis.text.y =element_text(size =14)) ancova_si +ancova_ci #sintaxe válida grazas ao paquete de R patchwork A.6. RANOVA ranovamates <- lmer(NotaMates ~(1|Escola), data = mates7, REML =FALSE) summary(ranovamates) sigma2_eps_ran =summary(ranovamates)$sigma^2 sigma2_u_ran =as.data.frame(summary(ranovamates)$varcor)[1,4] VT_ran =sigma2_u_ran +sigma2_eps_ran VPC_ran =sigma2_u_ran /VT_ran 66 A. Código de R A.6.1. VARCOMP e contraste sobre os efectos das escolas vcmates <- lm(NotaMates ~1,data = mates7) summary(vcmates) #Contraste sobre o efecto da escola mu_t <- vcmates$coef[[1]] #estimación media global (a da mostra) LR <- -2*logLik(vcmates)[1]-(-2*logLik(ranovamates))[1] 1-pchisq(LR, df =1)#p-valor A.6.2. Figura 3.5 x<- seq(0,45,length =1000) y<- dchisq(x, df =1) df <- data.frame(x,y) ggplot(df, aes(x=x,y= y)) + geom_line(color ="orange",size =1.25)+ xlab(" ")+ coord_cartesian(ylim =c(0,.15), clip ="off")+ geom_hline(yintercept =0,col="green2")+ geom_vline(xintercept =0,col="green2")+ geom_segment(x= LR, y=-0.004,xend = LR, yend =-.0115,col="tomato3", linetype ="solid",size =1.5)+ annotate("text",x=6,y=.1,parse =TRUE,label =expression(chi[1]^2), col ="orange",size =17)+ annotate("text",x= LR, y=-0.017,parse =TRUE,label =expression(LR[obs]), col ="tomato3",size =10)+ annotate("text",x=LR+6,y=.007,parse =TRUE,label = expression(1.96e-09), col ="lightsalmon",size =8)+ theme_classic() + theme( axis.title.x =element_text(size =10), axis.text.x =element_text(size =17), axis.title.y =element_blank(), axis.text.y =element_text(size =17)) A.6.3. Conxunto de datos u0df A.6. RANOVA 67 u0 <- ranef(ranovamates, postVar =TRUE) u0se <- sqrt(attr(u0[[1]], "postVar")[1, , ]) escolaidentif <- rownames(u0[[1]]) u0df <- cbind(escolaidentif, u0[[1]], u0se) colnames(u0df) <- c("escolaidentif","u0","u0se") u0df <- u0df[order(u0df$u0), ] u0df <- cbind(u0df, c(1:dim(u0df)[1])) colnames(u0df)[4]<- "u0pos" for (k in 1:nrow(u0df)){ numero <- as.numeric(strsplit(u0df$escolaidentif[k], split ="")[[1]][2]) u0df$escolaidentif[k] <- numero } u0df <- u0df[order(u0df$escolaidentif), ] #reordenación estloc <- summary(ranovamates)$coef[1]+u0df$u0 u0df$estloc <- estloc u0df$escolaidentif <- paste("E",1:7,sep ="") naive_res <- mu_local -summary(ranovamates)$coef[1] u0df$naive_res <- naive_res shrink <- numeric(7) for (i in 1:length(levels(mates7$Escola))){ shrink[i] <- u0[[1]][i, 1]/naive_res[i] } u0df$shrinkage <- shrink kable(u0df, format ="pipe",digits =3,row.names =FALSE) A.6.4. Figura 3.3: Gráfico de eiruga ggplot(u0df, aes(x= u0pos, y= u0)) + geom_segment(x= u0df$u0pos, y= u0df$u0 -1.96 *u0df$u0se, xend = u0df$u0pos, yend = u0df$u0 +1.96 *u0df$u0se, size =1,col ="darkcyan")+ geom_point(size =3,col ="tomato3")+ xlab("")+ ylab("Residuos a nivel de Escola")+ scale_x_continuous(breaks =1:7)+ geom_hline(yintercept =0,size =1)+ annotate("text",x= u0df$u0pos +0.3,y= u0df$u0, parse =TRUE,label = 68 A. Código de R paste("E",1:7,sep =""), col ="tomato3",size =7)+ theme( axis.title.x =element_text(size =19), axis.text.x =element_text(size =17), axis.title.y =element_text(size =19), axis.text.y =element_text(size =17)) A.6.5. Figura 3.2 O código de relativo á Figura 3.4 concernente ao modelo dun só nivel para a media é análogo ao disposto deseguido. mates7$Estudante <- 1:nrow(mates7) attach(mates7) eps_ran <- residuals(ranovamates) ggplot(mates7, aes(x= Estudante, y= NotaMates, color = Escola)) + geom_point(size =1.25)+ coord_cartesian(xlim =c(0,350), ylim =c(1,25)) + xlab("Indentificador de estudante")+ ylab("Nota en Matemáticas")+ scale_color_brewer(palette ="Dark2")+ #Media global geom_segment(x=-7,y=summary(ranovamates)$coef[1], xend =317, yend =summary(ranovamates)$coef[1], size =1.35,col =1)+ #Medias locais geom_segment(x=-7,y= estloc[3], xend =294,yend = estloc[3], lwd =1, linetype ="twodash",color = dark_2[3]) + geom_segment(x=-7,y= estloc[4], xend =294,yend = estloc[4], lwd =1, linetype ="twodash",color = dark_2[4]) + geom_segment(x=0,y= estloc[5], xend =294,yend = estloc[5], lwd =1, linetype ="twodash",color = dark_2[5]) + geom_segment(x=0,y= estloc[7], xend =294,yend = estloc[7], lwd =1, linetype ="twodash",color = dark_2[7]) + #Frechas geom_segment(x= Estudante[72], y= NotaMates[72], xend = Estudante[72], yend = NotaMates[72]-eps_ran[72], color =1,arrow =arrow(), A.6. RANOVA 69 size =1.15)+ geom_segment(x= Estudante[140], y= NotaMates[140], xend = Estudante[140], yend = NotaMates[140]-eps_ran[140], color =1,arrow =arrow(), size =1.15)+ geom_segment(x= Estudante[185], y= NotaMates[185], xend = Estudante[185], yend = NotaMates[185]-eps_ran[185], color =1,arrow =arrow(), size =1.15)+ geom_segment(x= Estudante[264], y= NotaMates[264], xend = Estudante[264], yend = NotaMates[264]-eps_ran[264], color =1,arrow =arrow(), size =1.15)+ #Erros annotate("text",x= Estudante[72]+20,y= (NotaMates[72]+estloc[3])/2, parse =TRUE,label =expression(widehat(epsilon)[7][","][3]), col =1,size =7)+ annotate("text",x= Estudante[140]-23,y= (NotaMates[140]+estloc[4])/2, parse =TRUE,label =expression(widehat(epsilon)[26][","][4]), col =1,size =7)+ annotate("text",x= Estudante[185]+24.5,y= (NotaMates[185]+estloc[5])/2 -.1,parse =TRUE,label =expression(widehat(epsilon)[37][","][5]), col =1,size =7)+ annotate("text",x= Estudante[264]-23,y= (NotaMates[264]+estloc[7])/2 -.5,parse =TRUE,label =expression(widehat(epsilon)[16][","][7]), col =1,size =7)+ #Observacións geom_point(aes(x= Estudante[72], y= NotaMates[72]), col = dark_2[3], size =5)+ geom_point(aes(x= Estudante[140], y= NotaMates[140]), col = dark_2[4], size =5)+ geom_point(aes(x= Estudante[185], y= NotaMates[185]), col = dark_2[5], size =5)+ geom_point(aes(x= Estudante[264], y= NotaMates[264]), col = dark_2[7], size =5)+ annotate("text",x= Estudante[72]-18,y= NotaMates[72]-.55,parse = TRUE,label =expression(Y[7][","][3]), col = dark_2[3], size =6)+ annotate("text",x= Estudante[140]-12,y= NotaMates[140]-.65,parse = TRUE,label =expression(Y[26][","][4]), col = dark_2[4], size =6)+ annotate("text",x= Estudante[185]-12,y= NotaMates[185]-.65,parse = 70 A. Código de R TRUE,label =expression(Y[37][","][5]), col = dark_2[5], size =6)+ annotate("text",x= Estudante[264]-12,y= NotaMates[264]-.65,parse = TRUE,label =expression(Y[16][","][7]), col = dark_2[7], size =6)+ #Media global símbolo annotate("text",x=329,y=summary(ranovamates)$coef[1], parse =TRUE, label =expression(widehat(mu)), size =7)+ #Medias locais símbolos annotate("text",x=330,y= estloc[3]+.25,parse =TRUE,label = expression(paste(widehat(mu) +widehat(u)[3], "=",widehat(mu)[3])), size =6,col = dark_2[3]) + annotate("text",x=330,y= estloc[4], parse =TRUE,label = expression(paste(widehat(mu) +widehat(u)[4], "=",widehat(mu)[4])), size =6,col = dark_2[4]) + annotate("text",x=330,y= estloc[5]-.1,parse =TRUE,label = expression(paste(widehat(mu) +widehat(u)[5], "=",widehat(mu)[5])), size =6,col = dark_2[5]) + annotate("text",x=330,y= estloc[7], parse =TRUE,label = expression(paste(widehat(mu) +widehat(u)[7], "=",widehat(mu)[7])), size =6,col = dark_2[7]) + #Residuos a nivel de escola u^gorro_{0j} geom_segment(x=-7,y=summary(ranovamates)$coef[1], xend =-7,yend = estloc[3], color = dark_2[3], arrow =arrow(length =unit(0.2,"inches")), linetype ="solid",size =.9)+ geom_segment(x=-7,y=summary(ranovamates)$coef[1], xend =-7,yend = estloc[4], color = dark_2[4], arrow =arrow(length =unit(0.2,"inches")), linetype ="solid",size =.9)+ geom_segment(x=20,y=summary(ranovamates)$coef[1], xend =20,yend = estloc[5], color = dark_2[5], arrow =arrow(length =unit(0.15,"inches")), linetype ="solid",size =.75)+ geom_segment(x=20,y=summary(ranovamates)$coef[1], xend =20,yend = estloc[7], color = dark_2[7], arrow =arrow(length =unit(0.15,"inches")), linetype ="solid",size =.75)+ #Símbolos residuos a nivel de escola u^gorro_{0j} annotate("text",x=5,y= (summary(ranovamates)$coef[1]+estloc[3]) /2 -.3,parse =TRUE,label =expression(widehat(u)[3]), col = dark_2[3], size =6.5)+ annotate("text",x=5,y= (summary(ranovamates)$coef[1]+estloc[4]) /2, A.7. Modelos mixtos con covariables relativas ao nivel 171 parse =TRUE,label =expression(widehat(u)[4]), col = dark_2[4], size =6.5)+ annotate("text",x=36,y= (summary(ranovamates)$coef[1]+estloc[5]) /2, parse =TRUE,label =expression(widehat(u)[5]), col = dark_2[5], size =6.5)+ annotate("text",x=36,y= (summary(ranovamates)$coef[1]+estloc[7]) /2, parse =TRUE,label =expression(widehat(u)[7]), col = dark_2[7], size =6.5)+ theme( legend.title =element_text(size =18), legend.text =element_text(size =16), axis.title.x =element_blank(), axis.text.x =element_blank(), axis.title.y =element_text(size =18), axis.text.y =element_text(size =16)) A.7. Modelos mixtos con covariables relativas ao nivel 1 A.7.1. Modelo con intercepto aleatorio e pendente fixa matesmm1 <- lmer(NotaMates ~StSE +(1|Escola), data = mates7, REML =FALSE) summary(matesmm1) coefsmm1 <- summary(matesmm1)$coef sigma2_eps_mm1 <- summary(matesmm1)$sigma^2 sigma2_u_mm1 <- as.data.frame(summary(matesmm1)$varcor)[1,4] VT_mm1 <- sigma2_u_mm1 +sigma2_eps_mm1 VPC_mm1 <- sigma2_u_mm1 /VT_mm1 #Cálculo dos efectos aleatorios e dos erros (nivel 1) u0_mm1 <- ranef(matesmm1, postVar =TRUE) uj_mm1 <- u0_mm1[[1]][,1] u0se_mm1 <- sqrt(attr(u0_mm1[[1]], "postVar")[1, , ]) eps_mm1 <- residuals(matesmm1) A.7.2. Modelo con intercepto aleatorio e pendente fixa con SocialMin 72 A. Código de R matesmm2 <- lmer(NotaMates ~StSE +SocialMin +(1|Escola), data = mates7, REML =FALSE) summary(matesmm2) coefsmm2 <- summary(matesmm2)$coef sigma2_eps_mm2 <- summary(matesmm2)$sigma^2 sigma2_u0_mm2 <- as.data.frame(summary(matesmm2)$varcor)[1,4] VT_mm2 <- sigma2_u0_mm2 +sigma2_eps_mm2 VPC_mm2 <- (sigma2_u0_mm2) /VT_mm2 #Cálculo dos efectos aleatorios e dos erros (nivel 1) u0_mm2 <- ranef(matesmm2, postVar =TRUE) uj_mm2 <- u0_mm2[[1]][,1] u0se_mm2 <- sqrt(attr(u0_mm2[[1]], "postVar")[1, , ]) eps_mm2 <- residuals(matesmm2) A.7.3. Modelo con intercepto e pendente aleatorios e mais variable categórica matesmm3 <- lmer(NotaMates ~StSE +SocialMin +(1+StSE |Escola), REML =TRUE) summary(matesmm3) coefsmm3 <- summary(matesmm3)$coef sigma2_eps_mm3 <- summary(matesmm3)$sigma^2 sigma2_u0_mm3 <- as.data.frame(summary(matesmm3)$varcor)[1,4] sigma2_u1_mm3 <- as.data.frame(summary(matesmm3)$varcor)[2,4] sigma_u01_mm3 <- as.data.frame(summary(matesmm3)$varcor)[3,4] VT_mm3 <- sigma2_u0_mm3 +sigma2_u1_mm3 +sigma2_eps_mm3 VPC_mm3 <- (sigma2_u0_mm3 +sigma2_u1_mm3) /VT_mm3 #Cálculo dos efectos aleatorios e dos erros (nivel 1) u0_mm3 <- ranef(matesmm3, postVar =TRUE) u0j_mm3 <- u0_mm3[[1]][,1] u1j_mm3 <- u0_mm3[[1]][,2] eps_mm3 <- residuals(matesmm3) A.7.4. Figura 4.3 O código de relativo á Figura 4.1 concernente ao modelo mixto con intercepto aleatorio e pendente fixa matesmm1 (construído en A.7.1), e mais o código relativo á Figura 4.2 concernente A.7. Modelos mixtos con covariables relativas ao nivel 173 ao modelo mixto matesmm2 (construído en A.7.2) son análogos ao código disposto deseguido. rectamm3 <- function(x=0,escola =NULL,sm =0){ if (!sm %in% 0:1){ stop("Só hai dúas posibilidades para sm, non (0) ou si (1).") } if (!is.null(escola)){ if (!escola %in% 1:7){ stop("Tal escola non existe, só do 1 ao 7.") } } if (is.null(escola)){ if (sm == 0){ coefsmm3[1]+coefsmm3[2]*x } else{ coefsmm3[1]+coefsmm3[3]+coefsmm3[2]*x } } else{ if (sm == 0){ coefsmm3[1]+u0j_mm3[escola] +(coefsmm3[2]+u1j_mm3[escola]) *x } else{ coefsmm3[1]+coefsmm3[3]+u0j_mm3[escola] + (coefsmm3[2]+u1j_mm3[escola]) *x } } } ggplot(mates7, aes(x= StSE, y= NotaMates, color = Escola, shape = SocialMin)) + geom_point(size =1.25)+ coord_cartesian(xlim =c(range(StSE)[1], range(StSE)[2]), ylim =c(-.5,25))+ xlab("Status socio-económico")+ ylab("Nota en Matemáticas")+ labs(shape ="Racial\nminori-\ntario")+ scale_color_brewer(palette ="Dark2")+ 80 BIBLIOGRAFÍA [13] Goldstein, H.: Multilevel statistical models. John Wiley & Sons (2011). [14] Goldstein, H., Sc, B.: Models for reality: New approaches to the understanding of educational processes. University of London, Institute of Education (1998). [15] Grilli, L., Rampichini, C.: A handful of critical choices in multilevel modelling. BEIO, Boletín de Estadística e Investigación Operativa 34(1), 7–24 (2018). [16] Hedges, L.V., Hedberg, E.C.: Intraclass correlation values for planning group-randomized trials in education. Educational Evaluation and Policy Analysis 29(1), 60–87 (2007). [17] Holm, S.: A simple sequentially rejective multiple test procedure. Scandinavian journal of statistics pp. 65–70 (1979). [18] Hox, J.J., Moerbeek, M., Van de Schoot, R.: Multilevel analysis: Techniques and applications. Routledge (2017). [19] Ibarrola, R.V., Pérez, A.G.: Principios de inferencia estadística. UNED, Universidad Nacional de Educación a Distancia (2006). [20] Kutner, M.H., Nachtsheim, C.J., Neter, J., Li, W., et al.: Applied linear statistical models. McGraw-Hill New York (2005). [21] Leyland, A.H., Groenewegen, P.P.: Multilevel modelling for public health and health services research: health in context. Springer Nature (2020). [22] McNeish, D.: Small sample methods for multilevel modeling: A colloquial elucidation of REML and the Kenward-Roger correction. Multivariate Behavioral Research 52(5), 661– 670 (2017). [23] Petrov, V.V., Ernesto, M.P.: Teoría de la probabilidad. 519.2 PET. Dirac (2008). [24] Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., Team, R.C.: Linear and nonlinear mixed effects models. R package version 3(57), 1–89 (2007). URL https://cran.r-project.org/ web/packages/nlme/nlme.pdf. [25] Rabe-Hesketh, S., Skrondal, A.: Multilevel and longitudinal modeling using Stata. STATA press (2008). [26] Rao, C.R.: Linear statistical inference and its applications, vol. 2. Wiley New York (1973). [27] Roback, P., Legler, J.: Beyond multiple linear regression: applied generalized linear models and multilevel models in R. Chapman and Hall/CRC (2021). URL https://bookdown. org/roback/bookdown-BeyondMLR/. BIBLIOGRAFÍA 81 [28] Scott, M.A., Simonoff, J.S., Marx, B.D.: The SAGE handbook of multilevel modeling. Sage (2013). [29] Snijders, T.A., Bosker, R.J.: Multilevel analysis: An introduction to basic and advanced multilevel modeling. sage (2011). [30] Steele, F.: Module 5: Introduction to Multilevel modelling concepts. LEMMA (Learning Environment for Multilevel Methodology and Applications), Centre for Multilevel Modelling, University of Bristol (2008). [31] Wickham, H.: Advanced R. CRC press (2019). URL https://adv-r.hadley.nz/. [32] Wickham, H., Grolemund, G.: R for data science: import, tidy, transform, visualize, and model data. .Reilly Media, Inc. (2016). URL https://r4ds.had.co.nz/. [33] Xie, Y.: Dynamic Documents with R and knitr. Chapman and Hall/CRC (2017). URL https://yihui.org/knitr/. [34] Zuur, A.F., Ieno, E.N., Walker, N.J., Saveliev, A.A., Smith, G.M., et al.: Mixed effects models and extensions in ecology with R, vol. 574. Springer (2009).