Apresentação dos dados

Este estudo tem como objetivo modelar os dados de diplomados por curso, na UFERSA. Considerando que os números de diplomados por curso representam contagens, um método apropriado para descrever este conjunto de dados é a regressão Poisson, um dos casos de Modelos Lineares Generalizados, em que a variável resposta segue distribuição Poisson.

Os dados estão disponíveis no Relatório de Gestão 2017 e no pacote dadosUfersa. Para baixar o pacote diretamente do repositório do github, instale o pacote devtools e utilize os comandos abaixo:

install.packages("devtools")
devtools::install_github("Kassio-Ferreira/dadosUfersa")

O pacote dadosUfersa contém uma séries de conjuntos de dados retirados do Relatório de Gestão 2017, aqui utilizaremos os dados sobre o número de alunos diplomados. A tabela foi construída com base no item 7 do anexo Indicadores do TCU do Relatório de Gestão 2017 da UFERSA. Os cursos de engenharia não possuem concorrência disponível, pois os alunos ingressam no curso de Ciência e Tecnologia para depois, por meio de seleção interna, iniciarem os estudos em alguma área da engenharia.

library(dadosUfersa)
library(dplyr)
library(ggplot2)
library(summarytools)
library(MASS)
library(ggfortify)

data(diplomados)

diplomados %>% knitr::kable(., caption = "Estrutura dos dados.")
Estrutura dos dados.
curso campus duracao ingressos enade area turno candidatos vagas concorrencia diplomados
agronomia mossoro 5 123 4.00 agrarias integral 5268 120 43.90000 43
direito mossoro 5 98 5.00 humanas noturno 7304 80 91.30000 50
engenharia agricola e ambiental mossoro 5 37 3.00 agrarias integral 635 50 12.70000 1
engenharia de pesca mossoro 5 51 NA agrarias integral 1546 30 51.53333 10
engenharia florestal mossoro 5 50 4.00 agrarias integral 1954 50 39.08000 9
medicina veterinaria mossoro 5 57 4.00 biologicas integral 2690 50 53.80000 31
zootecnia mossoro 5 52 3.00 agrarias integral 3210 50 64.20000 10
biotecnologia mossoro 4 51 NA biologicas integral 7552 50 151.04000 24
administracao mossoro 4 121 3.00 humanas noturno 6793 100 67.93000 31
ciencia da computacao mossoro 4 54 4.00 tecnologicas integral 5638 50 112.76000 16
ciencias contabeis mossoro 4 84 4.00 humanas noturno 1433 200 7.16500 34
computacao e informatica angicos 4 54 4.00 tecnologicas noturno 1694 50 33.88000 5
ecologia mossoro 4 50 NA biologicas integral 2155 50 43.10000 4
educacao no campo mossoro 4 59 NA interdisciplinar integral 120 60 2.00000 0
sistemas de informacao angicos 4 51 3.00 tecnologicas noturno 1493 50 29.86000 8
ciencia e tecnologia angicos 3 200 3.00 tecnologicas integral 2719 200 13.59500 82
ciencia e tecnologia caraubas 3 200 3.33 tecnologicas integral 2927 200 14.63500 48
ciencia e tecnologia mossoro 3 425 3.80 tecnologicas integral 5608 400 14.02000 213
ciencia e tecnologia pau dos ferros 3 200 2.50 tecnologicas integral 3282 200 16.41000 67
ciencia e tecnologia angicos 3 100 3.00 tecnologicas noturno 1766 100 17.66000 17
ciencia e tecnologia caraubas 3 101 3.33 tecnologicas noturno 1539 100 15.39000 13
ciencia e tecnologia mossoro 3 227 3.80 tecnologicas noturno 3753 80 46.91250 64
ciencia e tecnologia pau dos ferros 3 102 2.50 tecnologicas noturno 2289 100 22.89000 27
engenharia da computacao pau dos ferros 2 5 2.00 engenharias integral NA NA NA 3
engenharia civil angicos 2 37 3.00 engenharias integral NA NA NA 35
engenharia civil caraubas 2 29 3.00 engenharias integral NA NA NA 31
engenharia civil mossoro 2 60 4.00 engenharias integral NA NA NA 52
engenharia civil pau dos ferros 2 30 3.00 engenharias integral NA NA NA 37
engenharia eletrica mossoro 2 11 4.00 engenharias integral NA NA NA 4
engenharia eletrica caraubas 2 16 4.00 engenharias integral NA NA NA 13
engenharia de energia mossoro 2 11 2.00 engenharias integral NA NA NA 8
engenharia mecanica caraubas 2 11 3.00 engenharias integral NA NA NA 9
engenharia mecanica mossoro 2 57 3.00 engenharias integral NA NA NA 36
engenharia de petroleo mossoro 2 12 NA engenharias noturno NA NA NA 6
engenharia de producao angicos 2 11 3.00 engenharias integral NA NA NA 13
engenharia de producao mossoro 2 30 4.00 engenharias noturno NA NA NA 16
engenharia quimica mossoro 2 34 4.00 engenharias integral NA NA NA 41

Os cursos de Engenharia não possuem dados sobre a concorrência, uma vez que os estudantes ingressam, via SiSU, no curso de Ciência e Tecnologia e passam por seleção interna para ocupar as vagas nas Engenharias. Por este motivo, os cursos de Engenharia foram excluídos do modelo. Por outro lado, os cursos de Ciência e Tecnologia não possuem conceito no Enade. O conceito para estes cursos foi computado como e média dos conceitos das Engenharias no mesmo campus. São removidos os cursos com dados faltantes e as Engenharias, resultando num dataframe com 19 cursos. As variáveis apresentadas são:

diplomados = diplomados %>% filter(complete.cases(.))
diplomados %>% knitr::kable(., caption = "Dados após remoção dos itens faltantes.")
Dados após remoção dos itens faltantes.
curso campus duracao ingressos enade area turno candidatos vagas concorrencia diplomados
agronomia mossoro 5 123 4.00 agrarias integral 5268 120 43.9000 43
direito mossoro 5 98 5.00 humanas noturno 7304 80 91.3000 50
engenharia agricola e ambiental mossoro 5 37 3.00 agrarias integral 635 50 12.7000 1
engenharia florestal mossoro 5 50 4.00 agrarias integral 1954 50 39.0800 9
medicina veterinaria mossoro 5 57 4.00 biologicas integral 2690 50 53.8000 31
zootecnia mossoro 5 52 3.00 agrarias integral 3210 50 64.2000 10
administracao mossoro 4 121 3.00 humanas noturno 6793 100 67.9300 31
ciencia da computacao mossoro 4 54 4.00 tecnologicas integral 5638 50 112.7600 16
ciencias contabeis mossoro 4 84 4.00 humanas noturno 1433 200 7.1650 34
computacao e informatica angicos 4 54 4.00 tecnologicas noturno 1694 50 33.8800 5
sistemas de informacao angicos 4 51 3.00 tecnologicas noturno 1493 50 29.8600 8
ciencia e tecnologia angicos 3 200 3.00 tecnologicas integral 2719 200 13.5950 82
ciencia e tecnologia caraubas 3 200 3.33 tecnologicas integral 2927 200 14.6350 48
ciencia e tecnologia mossoro 3 425 3.80 tecnologicas integral 5608 400 14.0200 213
ciencia e tecnologia pau dos ferros 3 200 2.50 tecnologicas integral 3282 200 16.4100 67
ciencia e tecnologia angicos 3 100 3.00 tecnologicas noturno 1766 100 17.6600 17
ciencia e tecnologia caraubas 3 101 3.33 tecnologicas noturno 1539 100 15.3900 13
ciencia e tecnologia mossoro 3 227 3.80 tecnologicas noturno 3753 80 46.9125 64
ciencia e tecnologia pau dos ferros 3 102 2.50 tecnologicas noturno 2289 100 22.8900 27

A variáveis categóricas devem ser declaradas como fatores:

diplomados$campus <- as.factor(diplomados$campus)
diplomados$duracao <- as.factor(diplomados$duracao)
diplomados$area <- as.factor(diplomados$area)
diplomados$turno <- as.factor(diplomados$turno)

Também é possível definir os níveis das variáveis categóricas que serão utilizados como referência. Queremos que a referência da variável “campus” seja Mossoró e da variável “area” seja Tecnológicas:

# colocando mossoró como nível de referencia da variável "campus"
diplomados <- within(diplomados, campus <- relevel(campus, ref = "mossoro"))

# colocando tecnologicas como nível de referencia da variável "area"
diplomados <- within(diplomados, area <- relevel(area, ref = "tecnologicas"))

É possível avaliar as observações obtidas por campus e área no boxplot abaixo:

df = diplomados[,c("area", "campus","diplomados")]
p <- ggplot(data = df, aes(x=area, y=diplomados)) + geom_boxplot(aes(fill=campus), alpha=0.3) +
  geom_jitter() +  facet_wrap( ~ area, scales="free")

p

Segue tabela com a análise descritiva dos dados:

print(dfSummary(diplomados, graph.magnif = 0.75), method = 'render')

Data Frame Summary

diplomados

N: 19
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 curso [character] 1. ciencia e tecnologia 2. administracao 3. agronomia 4. ciencia da computacao 5. ciencias contabeis 6. computacao e informatica 7. direito 8. engenharia agricola e amb 9. engenharia florestal 10. medicina veterinaria [ 2 others ] 8 (42.1%) 1 (5.3%) 1 (5.3%) 1 (5.3%) 1 (5.3%) 1 (5.3%) 1 (5.3%) 1 (5.3%) 1 (5.3%) 1 (5.3%) 2 (10.5%) 19 (100%) 0 (0%)
2 campus [factor] 1. mossoro 2. angicos 3. caraubas 4. pau dos ferros 11 (57.9%) 4 (21.1%) 2 (10.5%) 2 (10.5%) 19 (100%) 0 (0%)
3 duracao [factor] 1. 3 2. 4 3. 5 8 (42.1%) 5 (26.3%) 6 (31.6%) 19 (100%) 0 (0%)
4 ingressos [numeric] mean (sd) : 122.95 (94.35) min < med < max : 37 < 100 < 425 IQR (CV) : 107.5 (0.77) 16 distinct values 19 (100%) 0 (0%)
5 enade [numeric] mean (sd) : 3.49 (0.64) min < med < max : 2.5 < 3.33 < 5 IQR (CV) : 1 (0.18) 2.50 : 2 (10.5%) 3.00 : 6 (31.6%) 3.33 : 2 (10.5%) 3.80 : 2 (10.5%) 4.00 : 6 (31.6%) 5.00 : 1 (5.3%) 19 (100%) 0 (0%)
6 area [factor] 1. tecnologicas 2. agrarias 3. biologicas 4. humanas 11 (57.9%) 4 (21.1%) 1 (5.3%) 3 (15.8%) 19 (100%) 0 (0%)
7 turno [factor] 1. integral 2. noturno 10 (52.6%) 9 (47.4%) 19 (100%) 0 (0%)
8 candidatos [numeric] mean (sd) : 3262.89 (1953) min < med < max : 635 < 2719 < 7304 IQR (CV) : 2780.5 (0.6) 19 distinct values 19 (100%) 0 (0%)
9 vagas [numeric] mean (sd) : 117.37 (89.12) min < med < max : 50 < 100 < 400 IQR (CV) : 110 (0.76) 50 : 7 (36.8%) 80 : 2 (10.5%) 100 : 4 (21.1%) 120 : 1 (5.3%) 200 : 4 (21.1%) 400 : 1 (5.3%) 19 (100%) 0 (0%)
10 concorrencia [numeric] mean (sd) : 37.79 (29.24) min < med < max : 7.17 < 29.86 < 112.76 IQR (CV) : 35.34 (0.77) 19 distinct values 19 (100%) 0 (0%)
11 diplomados [numeric] mean (sd) : 40.47 (47.73) min < med < max : 1 < 31 < 213 IQR (CV) : 37.5 (1.18) 18 distinct values 19 (100%) 0 (0%)

Generated by summarytools 0.8.8 (R version 3.5.1)
2018-12-10

Regressão Poisson

Neste caso, temos que a variável resposta \(Y = diplomados\) é uma contagem (inteiro positivo). Para explicar \(Y\) em termos de um conjunto de preditores, podemos utilizar um modelo de regressão para contagem. Os casos especiais de Modelos Lineares Generalizados para esta situação são os modelos Poisson e Binomial Negativo.

A distribuição de Poisson é utilizada para modelar contagens de eventos em um certo período de tempo. Sua função massa de probabilidade é dada por

\[\begin{align} f(x;\mu) &= \frac{e^{-\mu}\mu^x}{x!} \end{align}\]

com \(x=1,2,\ldots\). A média e a variância são dadas por \(\mathbb{E}(x) = var(x) = \mu\). Para modelar a variável \(Y\) em termos de um vetor de preditores \(x_i\), é necessária uma função de ligação entre \(\mathbb{E}(Y_i) = \mu_i\) e \(x_i\). Utiliza-se, portanto, a função de ligação log:

\[\begin{align} \log{\mu_i} &= \eta_i = x_i^\top\boldsymbol{\beta} \end{align}\]

A função abaixo fornece uma tabela com os intervalos de confiança para as rates ratio dos modelos.

glm.RR <- function(GLM.RESULT, digits = 2) {

    COEF      <- stats::coef(GLM.RESULT)
    CONFINT   <- stats::confint(GLM.RESULT)
    TABLE     <- cbind(coef=COEF, CONFINT)
    TABLE.EXP <- round(exp(TABLE), digits)

    colnames(TABLE.EXP)[1] <- "RR"

    TABLE.EXP
}
log.fit <- glm(diplomados ~ 
                 campus + 
                 enade + 
                 area +
                 concorrencia +
                 turno +
                 offset(log(ingressos)), family = poisson(link=log), data = diplomados)
log.fit %>% summary.glm()
## 
## Call:
## glm(formula = diplomados ~ campus + enade + area + concorrencia + 
##     turno + offset(log(ingressos)), family = poisson(link = log), 
##     data = diplomados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0542  -0.6387  -0.0932   0.3692   1.7363  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.922331   0.389697  -4.933 8.10e-07 ***
## campusangicos        -0.101135   0.132717  -0.762  0.44604    
## campuscaraubas       -0.588527   0.149108  -3.947 7.91e-05 ***
## campuspau dos ferros  0.130244   0.178032   0.732  0.46443    
## enade                 0.336429   0.103919   3.237  0.00121 ** 
## areaagrarias         -0.616973   0.144071  -4.282 1.85e-05 ***
## areabiologicas        0.132266   0.196085   0.675  0.49997    
## areahumanas           0.358555   0.140686   2.549  0.01081 *  
## concorrencia         -0.003062   0.001736  -1.764  0.07779 .  
## turnonoturno         -0.582647   0.103673  -5.620 1.91e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 132.652  on 18  degrees of freedom
## Residual deviance:  29.265  on  9  degrees of freedom
## AIC: 144.24
## 
## Number of Fisher Scoring iterations: 4

Os resultados mostram que foram signicativos, ao nível de 5%, o campus de Caraúbas, a áreas de ciências agrárias (em relação à categoria de referência) e o turno noturno como fatores que reduzem o número de diplomados. Os fatores estatisticamente significativos que aumentam o número de formandos são a área de humanas (em relação à categoria de referência) e a nota do enade. A análise dos resíduos após o ajuste da regressão Poisson mostram que há 3 outliers, isto é, resíduos de Pearson fora do intervalo \([-1.96, 1.96]\): as observações 1, 3 e 10. O QQ-Plot oferece indícios de que os resíduos são Normais.

log.fit %>% autoplot(., which = 1:4, label.size = 3, data = diplomados,
                     smooth.colour = 'black', smooth.linetype = 'dashed',
                     colour = 'campus')

O modelo Poisson assume que a variância e a média são iguais, mas nem sempre a variância observada nos dados corrobora com esta hipótese. Quando a variância observada é maior do que a variância assumida pelo modelo, a regressão Poisson torna-se inadequada. Este fenômeno é conhecido como overdispersion. Uma possível solução para o problema é o modelo quasipoisson:

Modelo Quasipoisson:

O modelo quasipoisson introduz um parâmetro de dispersão \(\phi\) no modelo Poisson, de tal forma que a variância agora é dada por

\[var(Y_i|\eta_i) = \phi \mu_i\]

Se \(\phi \geq 1\), então a variância cresce mais rapidamente do que a média. Um estimador de \(\phi\) é dado por

\[\hat{\phi} = \frac{1}{n-k}\sum\frac{(Y_i - \hat{\mu}_i)^2}{\hat{\mu}_i}\] em que \(n\) é o tamanho da amostra, \(k\) é o número de parâmetros estimados e \(\hat{\mu_i} = g^-1(\hat{\eta}_i)\), é a esperança ajustada de \(Y_i\), em que \(g^-1(.)\) é a inversa da função de ligação. Ao ajustar o modelo quasipoisson, notamos que as estimativas dos parâmetros serão as mesmas obtidos no modelo de Poisson, no entanto, os erros padrão dos estimadores serão maiores, resultando em intervalos de confiança com maior amplitude.

## quasi-Poisson to allow the scale parameter to change from 1
model.1q <- glm(diplomados ~ 
                  campus + 
                  enade + 
                  area +
                  concorrencia +
                  turno +
                  offset(log(ingressos)), family = quasipoisson(link=log), data = diplomados)

model.1q %>% summary.glm()
## 
## Call:
## glm(formula = diplomados ~ campus + enade + area + concorrencia + 
##     turno + offset(log(ingressos)), family = quasipoisson(link = log), 
##     data = diplomados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0542  -0.6387  -0.0932   0.3692   1.7363  
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)   
## (Intercept)          -1.922331   0.639056  -3.008  0.01476 * 
## campusangicos        -0.101135   0.217640  -0.465  0.65319   
## campuscaraubas       -0.588527   0.244519  -2.407  0.03945 * 
## campuspau dos ferros  0.130244   0.291952   0.446  0.66605   
## enade                 0.336429   0.170415   1.974  0.07981 . 
## areaagrarias         -0.616973   0.236260  -2.611  0.02821 * 
## areabiologicas        0.132266   0.321555   0.411  0.69045   
## areahumanas           0.358555   0.230708   1.554  0.15457   
## concorrencia         -0.003062   0.002847  -1.075  0.31014   
## turnonoturno         -0.582647   0.170012  -3.427  0.00754 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 2.689209)
## 
##     Null deviance: 132.652  on 18  degrees of freedom
## Residual deviance:  29.265  on  9  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
#model.1q %>% autoplot(., which = 1:4, label.size = 3, data = diplomados,
#                      smooth.colour = 'black', smooth.linetype = 'dashed',
#                      colour = 'campus')

Analisando os resultados, vemos que ao nível de \(\alpha = 0.05\), foram significantes o campus Caraúbas, a área de ciências agrárias e o turno noturno como fatores que reduzem o total de diplomados. Ao nível de 10%, a nota do curso no enade também foi significativa, com notas maiores contribuindo para mais diplomados.

Note que, a análise dos resíduos é a mesma para os dois modelos, uma vez que as estimativas dos parâmetros são iguais.

data.frame(curso = diplomados$curso,
           campus = diplomados$campus,
           turno = diplomados$turno,
           n_diplomados = diplomados$diplomados,
           poisson = round(fitted(log.fit),3),
           quasipoisson = round(fitted(model.1q),3)) %>% knitr::kable(.,caption = "Resultados")
Resultados
curso campus turno n_diplomados poisson quasipoisson
agronomia mossoro integral 43 32.596 32.596
direito mossoro noturno 50 46.580 46.580
engenharia agricola e ambiental mossoro integral 1 7.706 7.706
engenharia florestal mossoro integral 9 13.447 13.447
medicina veterinaria mossoro integral 31 31.000 31.000
zootecnia mossoro integral 10 9.250 9.250
administracao mossoro noturno 31 31.522 31.522
ciencia da computacao mossoro integral 16 21.480 21.480
ciencias contabeis mossoro noturno 34 36.899 36.899
computacao e informatica angicos noturno 5 13.803 13.803
sistemas de informacao angicos noturno 8 9.427 9.427
ciencia e tecnologia angicos integral 82 69.583 69.583
ciencia e tecnologia caraubas integral 48 47.606 47.606
ciencia e tecnologia mossoro integral 213 213.847 213.847
ciencia e tecnologia pau dos ferros integral 67 73.483 73.483
ciencia e tecnologia angicos noturno 17 19.188 19.188
ciencia e tecnologia caraubas noturno 13 13.394 13.394
ciencia e tecnologia mossoro noturno 64 57.672 57.672
ciencia e tecnologia pau dos ferros noturno 27 20.517 20.517

Abaixo temos os intervalos de confiança fornecidos para as rate ratios ao ajustar os modelos Poisson e Quasipoisson. Note que os intervalos de confiança da Quasipoisson são maiores, dado o tratamento que o método dá ao problema de overdispersion:

knitr::kable(as.data.frame(glm.RR(log.fit, 3)), 
             caption = "Intervalo de confiança para RR - Poisson")
Intervalo de confiança para RR - Poisson
RR 2.5 % 97.5 %
(Intercept) 0.146 0.068 0.312
campusangicos 0.904 0.695 1.170
campuscaraubas 0.555 0.411 0.738
campuspau dos ferros 1.139 0.803 1.614
enade 1.400 1.145 1.721
areaagrarias 0.540 0.404 0.711
areabiologicas 1.141 0.763 1.650
areahumanas 1.431 1.085 1.885
concorrencia 0.997 0.993 1.000
turnonoturno 0.558 0.454 0.682 
knitr::kable(as.data.frame(glm.RR(model.1q, 3)), 
             caption = "Intervalo de confiança para RR - Quasipoisson")
Intervalo de confiança para RR - Quasipoisson
RR 2.5 % 97.5 %
(Intercept) 0.146 0.041 0.505
campusangicos 0.904 0.585 1.376
campuscaraubas 0.555 0.336 0.881
campuspau dos ferros 1.139 0.641 2.018
enade 1.400 1.008 1.971
areaagrarias 0.540 0.332 0.842
areabiologicas 1.141 0.577 2.060
areahumanas 1.431 0.908 2.249
concorrencia 0.997 0.991 1.002
turnonoturno 0.558 0.397 0.774