4.2.1 Especificação dos pesos espaciais

Apesar de algumas análises exploratórias no espaço não exigirem a matriz W de vizinhança (matriz de pesos espaciais), optou-se inicialmente por especificá-la. Embora se trate de um processo ad hoc, foram seguidos os procedimentos recomendados pela literatura (^{Almeida, 2012}; ^{Anselin; Rey, 2014}; ^{Golgher, 2015}), como a busca pelo tipo de matriz que maximiza a autocorrelação espacial (evidência de interação) presente nos dados (^{Vieira, 2009}). Para isso, para cada uma das duas variáveis dependentes (MAT e POR), os modelos descritos na Tabela 6 e na Tabela 8 foram estimados por Mínimos Quadrados Ordinários (OLS). A partir dos resíduos obtidos, calcularam-se e avaliaram-se os índices de Moran globais para diversas matrizes de pesos espaciais.

Iniciou-se pelas matrizes contíguas mais simples, Rainha e Torre. Os resultados foram similares e indicaram evidências de autocorrelação espacial positiva para as duas variáveis dependentes (Tabela 5). Essas duas matrizes geraram número irrelevantes de ilhas e um número parcimonioso de vizinhos. Com base na média de vizinhos desses dois tipos de matrizes espaciais (5,27 e 5,14, respectivamente), foram simulados três outros tipos de matrizes, com base nos quatro, cinco e seis vizinhos mais próximos (K-nearest neighbors [KNN]). Nessa abordagem, não se geram ilhas, mas o número de vizinhos é fixo. Nessas três matrizes, os valores do índice de Moran foram similares, indicando a estabilidade dos resultados. Adicionalmente, foram simulados três cenários extremos, com 10, 15 e 20 vizinhos mais próximos, e os valores do índice de Moran caíram ligeiramente.

Para complementar a busca simulamos três outros tipos de matrizes baseados em vizinhos situados a até 282, 80 e 40 km de distância. O primeiro tipo, KM-282, representa a distância mínima necessária para evitar a formação de ilhas. Esse valor foi sugerido pelo próprio software GeoDa, com base nas variáveis de latitude e longitude, e está alinhado com outras evidências empíricas sobre a mesma variável dependente (^{Fujita, Bagolin e Fochezatto, 2021}). Os dois outros tipos, KM-80 e KM-40, foram considerados com base na inspeção do correlograma espacial. Na primeira matriz (KM-282), os valores do índice de Moran diminuíram substancialmente. Nas duas últimas (KM-80 e KM-40), apesar de isso não ocorrer, os números médio e máximo de vizinhos aumentaram significativamente. No caso da matriz KM-40, o número de ilhas tornou-se inaceitável.

Com base nessas análises, optou-se pela matriz W mais parcimoniosa: Rainha. ^{Golgher (2015)} sugere que, em casos de dúvida, deve-se escolher uma matriz simples, como a de contiguidade. O autor também relata que, em trabalhos empíricos em estudos regionais, o número de conexões (vizinhos) varia entre 3,8 e 6,0 (^{Golgher, 2015}). Há algumas evidências empíricas favoráveis às estimações com matrizes menos conectadas, como as de contiguidade ou aquelas com os cinco ou sete vizinhos mais próximos (^{Golgher, 2015}). Em trabalho semelhante, ^{Cavalcanti et al. (2020)} consideraram uma matriz com os cinco vizinhos mais próximos, valor que é praticamente o mesmo da média de vizinhos (5,27) da matriz W considerada neste estudo.

4.2.2 Mapas descritivos das variáveis dependentes

Dado que o objetivo do trabalho é realizar uma análise espacial, focada principalmente na variação da proficiência em Matemática e Língua Portuguesa entre os ciclos de avaliação do Saeb de 2013 e 2017, as análises descritivas se concentram em ilustrações cartográficas das variáveis dependentes. A Figura 1 e a Figura 2 apresentam a distribuição das variações da proficiência em Matemática e em Língua Portuguesa, respectivamente. Ambas estão divididas em seis classes uniformes, representando os valores observados da variável (a) e seu desvio padrão (b).

Fonte: Elaborada pelos autores.

Figura 1 Mapa quantílico (a) e de desvio padrão (b) da variação da proficiência em Matemática. 

Fonte: Elaborada pelos autores.

Figura 2 Mapa quantílico (a) e de desvio padrão (b) da variação da proficiência em Língua Portuguesa. 

Sem a intenção de identificar padrões nos mapas apresentados na Figura 1 e na Figura 2, alguns fatos estilizados se destacam e até fazem sentido à luz da percepção empírica (^{Cavalcanti, André e Araújo, 2020}; ^{Fujita, Bagolin e Fochezatto, 2021}). Entre eles, estão:

4.2.3 Mapas Lisa das variáveis dependentes

Para inferências mais contundentes sobre padrões ou clusters espaciais, deve-se avaliar o indicador de correlação espacial de Moran local, que é materializado nos mapas Lisa. A Figura 3 e a Figura 4 apresentam o mapa Lisa dos regimes espaciais (a) e das indicações de significâncias (p-valor < 0,01) para os indicadores de Moran local (b) relacionados às variações da proficiência em Matemática (Figura 3) e em Língua Portuguesa (Figura 4). Os valores em vermelho no mapa Lisa (a) são os hot spots (high-high), onde altas variações da proficiência em Matemática ou em Língua Portuguesa são acompanhadas por altas variações nos vizinhos. Já os valores em azul representam os cold spots (low-low), onde baixas variações da proficiência em Matemática ou em Língua Portuguesa são acompanhadas por baixas variações dos vizinhos. Geralmente, os regimes low-high ou high-low são considerados outliers espaciais, conforme descrito por ^{Anselin e Rey (2014)}.

Fonte: Elaborada pelos autores.

Figura 3 Mapa Lisa (a) e de significância (p < 0,01) (b) da variação em Matemática. 

Fonte: Elaborada pelos autores.

Figura 4 Mapa Lisa (a) e de significância (p < 0,01) (b) da variação em Português. 

De forma geral, percebem-se mais pontos de cold spots do que de hot spots. Isso é mais evidente para a variação da proficiência em Matemática (n = 211 vs. n = 125) do que para a da proficiência em Língua Portuguesa (n = 173 vs. n = 130). A maior concentração de pontos low-low, no caso da proficiência em Língua Portuguesa (Figura 4), parece ocorrer no Noroeste de Minas Gerais e na Bahia. Quanto à variação da proficiência em Matemática, além dos clusters espaciais nessas mesmas localidades, também parece haver outro na Região Norte (Figura 3). Em relação aos hot spots, percebe-se que há uma concentração de pontos high-high no Ceará, tanto para a variação da proficiência em Matemática como para a da proficiência em Língua Portuguesa. No caso da variação da proficiência em Matemática, também existe um cluster espacial em Rondônia (Figura 3).

4.2.4 Modelos espaciais

Cabe ressaltar que, para a estimação dos modelos, utilizaram-se as variações das variáveis entre 2017 e 2013. Esse procedimento, que é recomendado quando se têm dados longitudinais com apenas dois pontos no tempo, visa eliminar as influências de fatores não observados que não mudam com o tempo, para melhor controle das variáveis omitidas. Com isso, utiliza-se, de fato, um tipo de estimador de efeito fixo, especificamente o método de primeiras diferenças para dados em painel (^{Almeida, 2012}). O resultado, por se ter apenas dois pontos no tempo, gera dados cross-section e, assim, podem-se aplicar os procedimentos usuais (clássicos) para a modelagem de dados espaciais (^{Almeida, 2012}; ^{Golgher, 2015}).

Os componentes espaciais do modelo podem aparecer basicamente em três formas: 1) na forma de defasagem espacial na variável dependente (ρWy); 2) na forma de defasagem nas variáveis explicativas (θWx); 3) ou como defasagem no termo de erro (λWu). Tais componentes podem aparecer de forma isolada ou em conjunto (^{Vieira, 2009}). Foram considerados para as análises espaciais apenas modelos que levassem em conta a variável dependente espacialmente defasada (ρWy) e/ou o erro espacialmente defasado (λWu). Isso porque estima-se, e as evidências sugerem, que as variáveis independentes consideradas – principalmente aquelas relacionadas com a percepção docente sobre a escola e as condições extraescolares e do aluno – não exercem influência espacial (spillover) na proficiência em Língua Portuguesa ou em Matemática do município.

A partir dessa decisão, temos menor risco ao excluir modelos espaciais mais complexos para estimar (SDM, SLX etc.) e ficamos restritos aos modelos SAR, SEM e SAC. Assim, seguindo o procedimento clássico (^{Almeida, 2012}) de especificação de modelos espaciais ou, em outros termos, a estratégica específica geral (^{Golgher, 2015}), começando pelo modelo OLS básico, estimamos os outros três modelos espaciais (SAR, SEM e SAC), um para cada variável dependente (MAT e POR). A partir dos testes diagnósticos dos resíduos do modelo OLS e de medidas de ajuste recomendadas, tomamos as decisões sobre os estimadores utilizados e os modelos finais para interpretação.

De forma geral, algumas evidências foram comuns e contundentes em todas as análises (modelos) realizadas: 1) rejeita-se veemente a hipótese de normalidade dos resíduos, conforme apontado pelo teste Jarque-Bera; 2) rejeita-se a hipótese de homogeneidade dos resíduos, conforme evidenciado pelo teste Koenker-Basset; 3) há evidências de algum tipo de dependência espacial, como apontado pelos outros testes diagnósticos (I de Moran dos resíduos e testes LM). O primeiro e o segundo pontos levam à decisão de considerar os métodos 2SLS/GMM, em detrimento do método ML, para estimar os parâmetros dos modelos. Isso se dá pelo fato de o método ML requerer normalidade e pela evidência de heterogeneidade, que só poderia ser tratada nos métodos 2SLS/GMM, a partir das correções dos erros padrão HAC e do procedimento de Kelejian e Prucha. Optou-se por seguir o caminho sugerido por ^{Anselin (1988)}: tratar os problemas gerados pela heterogeneidade espacial com o emprego de instrumentos fornecidos pela economia clássica, como o uso de erros padrão robustos (^{Vieira, 2009}), em detrimento de um modelo de heterogeneidade espacial. O terceiro ponto indica a pertinência dos modelos espaciais em relação aos dados, em detrimento de um modelo OLS básico.

4.2.5 Proficiência em Matemática

Partindo do modelo específico, estimado por OLS, e diagnosticando seus resíduos, podem-se constatar, conforme a Tabela 6, alguns direcionamentos para ajustar o modelo espacial da variação da proficiência em Matemática: 1) o teste difuso I de Moran dos resíduos indicou algum tipo de dependência espacial (I = 0,22; Z = 23,51; p-valor < 0,001); 2) os testes focados LM direcionaram para um SAR ou SAC.

Essa última conclusão surgiu porque todos os testes LM, inclusive os robustos, foram altamente significativos. Os testes LMρ e LMρ robusto foram, no entanto, mais significativos (maior valor da estatística χ²), direcionando assim para um modelo SAR. Como, entretanto, também os testes LMλ e LMλ robusto – bem como o teste LM (Sarma), que testa, conjuntamente, ρ = 0 e λ = 0 – foram significativos, isso poderia direcionar para um modelo de ordem superior do tipo SAC. Esse último deve ser considerado apenas se houver ganhos efetivos em eliminar a dependência espacial que ficou nos resíduos (^{Anselin e Rey, 2014}).

Os testes Jarque-Bera e Koenker-Basset rejeitaram as hipóteses de normalidade e homogeneidade dos resíduos, respectivamente. De acordo com os valores da estatística χ², a falta de normalidade e a presença de heterogeneidade parecem ser exacerbadas, corroborando o uso de estimadores que não levem em conta a normalidade dos dados e cujos erros padrão possam ser corrigidos para a presença de heterogeneidade. Assim, os modelos espaciais estimados por 2SLS são apesentados na Tabela 7, inclusive o modelo SEM para comparabilidade.

Os erros padrão de todos os modelos, para as devidas inferências, foram ajustados em relação à heterogeneidade, seja pelo método HAC, seja pelo método de Kelejian e Prucha. A única medida de ajuste existente quando se usam estimadores 2SLS/GMM, o pseudo R², ratifica as suspeitas iniciais: os modelos SAR e SAC são os que melhor se ajustam aos dados.

Vale lembrar que, para as medidas AFD, ICG e IED, não foram consideradas as variáveis dummies AFD1, ICG1 e IED1, a fim de evitar colinearidade perfeita. Assim, essas variáveis servem como categorias base. Desse modo, tomando o modelo SAR como o modelo final, preferido em detrimento do modelo SAC, observa-se que, em municípios com maior proporção de “docentes com licenciatura em área diferente daquela em que leciona, ou com bacharelado nas disciplinas da base curricular comum e complementação pedagógica concluída em área diferente daquela em que leciona” (AFD3), em relação a “docentes com formação superior de licenciatura na mesma disciplina que lecionam, ou bacharelado na mesma disciplina com curso de complementação pedagógica concluído” (AFD1), a variação da proficiência em Matemática foi menor.

No mesmo sentido, em municípios com maior proporção de “docentes que não possuem curso superior completo” (AFD5), em comparação com “docentes com formação superior de licenciatura na mesma disciplina que lecionam, ou bacharelado na mesma disciplina com curso de complementação pedagógica concluído” (AFD1), a variação da proficiência em Matemática também foi menor. Isso ocorreu tanto para AFD3 quanto para AFD5 no modelo SAR, devido aos coeficientes significativos e negativos. Esses achados estão em consonância com ^{Fujita et al. (2021)} e ^{Vernier (2016)}, que relataram que uma maior proporção de docentes pós-graduados nos municípios se relacionou positivamente com a performance dos alunos na escola. Vale ressaltar que o valor dos coeficientes em modelos SAR não é interpretado em termos de variações marginais (derivadas parciais), como nos modelos OLS ou SEM, devido aos efeitos de spillovers de realimentação dos vizinhos, materializados no coeficiente ρ. Por essa razão, optou-se por realizar esse tipo de análise somente na seção seguinte.

4.2.6 Proficiência em Língua Portuguesa

No caso da variação da proficiência em Língua Portuguesa, o modelo preferido também foi do tipo SAR, pelos mesmos motivos e argumentos discutidos na seção anterior: 1) os testes Jarque-Bera e Koenker-Basset rejeitaram as hipóteses de normalidade e homogeneidade dos resíduos (Tabela 7), respectivamente, e todos os modelos espaciais foram estimados por 2SLS, conforme Tabela 8; 2) os erros padrão de todos os modelos foram ajustados quanto à heterogeneidade; 3) o teste I de Moran dos resíduos, os testes focados LM e a medida de ajuste pseudo R² direcionaram para os modelos SAR e SAC, sendo avaliado o primeiro (modelo final), porque, apesar da significância de λ no modelo SAC, o erro espacialmente defasado não melhorou a explicação da variação da proficiência em Língua Portuguesa. Por esse motivo, escolheu-se o modelo mais parcimonioso.

A interpretação do resultado significativo da variável AFD3 é feita da mesma forma que na seção anterior. No caso da variável AFD4, que também é significativa, observa-se que, em municípios com maior proporção de “docentes com formação superior não considerada nas categorias anteriores” (AFD4), em relação a “docentes com formação superior de licenciatura na mesma disciplina que lecionam, ou bacharelado na mesma disciplina com curso de complementação pedagógica concluído” (AFD1), a variação da proficiência em Língua Portuguesa foi menor. Parece-nos que esse achado não tem muita relevância prática, devido à definição da variável AFD4.

A significância, no entanto, da variável Aluno no modelo SAR (Tabela 9) – e, diga-se de passagem, em todos os outros modelos, inclusive nos modelos SEM e SAC de proficiência em Matemática (Tabela 7) –, com significância marginal (p-valor < 0,15) no modelo SAR (Tabela 7), traz implicações práticas interessantes. Em todos os modelos, o coeficiente da variável Aluno foi negativo e teve valores similares, indicando estabilidade dos resultados. Assim, conclui-se inicialmente que, em municípios onde a percepção docente sobre a indisciplina dos alunos (Q81) e o alto nível de faltas (Q82) é maior, a proficiência em Língua Portuguesa é menor. Isso também ocorre no caso da proficiência em Matemática, mas essa avaliação é marginalmente significativa. Esses achados sobre a percepção docente são consistentes com os de ^{Vidal e Vieira (2017)}, que argumentam que os fatores externos à escola são mais relevantes do que os internos para as deficiências de aprendizagem.

Devido aos spillovers de realimentação, o valor de - 0,382 do parâmetro da variável Aluno no modelo SAR deve ser avaliado em termos de efeito direto (ED), efeito indireto (EI) ou efeito total (ET), como propõe ^{Golgher (2015)}. Assim, obtiveram-se os efeitos ED = - 1,630, EI = - 1,426 e ET = - 3,056. Ao se considerar apenas o ET, observa-se que, em municípios com desvio padrão acima da média em relação aos outros municípios no escore do componente Aluno – indicando aumento da percepção docente sobre a indisciplina e o nível de faltas dos alunos –, a proficiência em Língua Portuguesa diminuiu aproximadamente três pontos.