## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## 
## Anexando pacote: 'modelr'
## 
## 
## O seguinte objeto é mascarado por 'package:broom':
## 
##     bootstrap

Dados da CAPES sobre avaliação da pós-graduação

A CAPES é um órgão do MEC que tem a atribuição de acompanhar a pós-graduação na universidade brasileira. Uma das formas que ela encontrou de fazer isso e pela qual ela é bastante criticada é através de uma avaliação quantitativa a cada x anos (era 3, mudou para 4).

Usaremos dados da penúltima avaliação da CAPES:

cacc_tudo = read_projectdata()

glimpse(cacc_tudo)
## Rows: 73
## Columns: 31
## $ Instituição                  <chr> "UNIVERSIDADE FEDERAL DO AMAZONAS", "UNIV…
## $ Programa                     <chr> "INFORMÁTICA (12001015012P2)", "CIÊNCIA D…
## $ Nível                        <int> 5, 4, 3, 3, 3, 5, 4, 3, 3, 3, 5, 3, 3, 3,…
## $ Sigla                        <chr> "UFAM", "UFPA", "UFMA", "UEMA", "FUFPI", …
## $ `Tem doutorado`              <chr> "Sim", "Sim", "Não", "Não", "Não", "Sim",…
## $ `Docentes colaboradores`     <dbl> 0.25, 5.50, 3.00, 6.25, 1.75, 2.00, 1.00,…
## $ `Docentes permanentes`       <dbl> 24.75, 14.00, 10.00, 14.00, 9.50, 20.75, …
## $ `Docentes visitantes`        <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.75, 0.50,…
## $ `Resumos em conf`            <int> 20, 23, 15, 5, 4, 10, 6, 136, 0, 24, 27, …
## $ `Resumos expandidos em conf` <int> 25, 24, 7, 10, 1, 68, 9, 13, 4, 6, 16, 5,…
## $ `Artigos em conf`            <int> 390, 284, 115, 73, 150, 269, 179, 0, 120,…
## $ Dissertacoes                 <int> 108, 77, 50, 25, 31, 75, 60, 129, 45, 3, …
## $ Teses                        <int> 14, 0, 0, 0, 0, 24, 5, 0, 0, 0, 29, 0, 0,…
## $ periodicos_A1                <int> 15, 19, 5, 1, 7, 21, 21, 0, 3, 8, 44, 0, …
## $ periodicos_A2                <int> 19, 21, 11, 1, 4, 32, 13, 0, 9, 2, 23, 2,…
## $ periodicos_B1                <int> 19, 38, 7, 3, 6, 26, 16, 2, 6, 4, 32, 4, …
## $ periodicos_B2                <int> 1, 12, 2, 6, 0, 0, 11, 0, 0, 2, 1, 0, 0, …
## $ periodicos_B3                <int> 3, 16, 2, 2, 3, 16, 15, 0, 4, 6, 9, 0, 2,…
## $ periodicos_B4                <int> 0, 4, 0, 3, 3, 0, 1, 3, 1, 6, 0, 0, 4, 5,…
## $ periodicos_B5                <int> 10, 16, 8, 4, 12, 4, 16, 2, 6, 2, 11, 0, …
## $ periodicos_C                 <int> 9, 34, 12, 5, 2, 3, 11, 9, 5, 10, 16, 1, …
## $ periodicos_NA                <int> 7, 15, 8, 11, 12, 6, 19, 31, 7, 14, 19, 0…
## $ per_comaluno_A1              <int> 4, 1, 0, 0, 1, 7, 5, 0, 1, 0, 10, 0, 0, 2…
## $ per_comaluno_A2              <int> 5, 5, 5, 0, 2, 15, 3, 0, 3, 0, 3, 0, 0, 1…
## $ per_comaluno_B1              <int> 4, 2, 5, 2, 2, 14, 6, 0, 2, 0, 17, 0, 1, …
## $ per_comaluno_B2              <int> 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,…
## $ per_comaluno_B3              <int> 2, 2, 0, 1, 0, 7, 9, 0, 2, 0, 4, 0, 0, 1,…
## $ per_comaluno_B4              <int> 0, 0, 0, 0, 2, 0, 1, 0, 1, 3, 0, 0, 2, 0,…
## $ per_comaluno_B5              <int> 5, 0, 4, 0, 8, 3, 6, 0, 4, 0, 4, 0, 2, 5,…
## $ per_comaluno_C               <int> 6, 5, 3, 1, 2, 3, 7, 1, 2, 4, 8, 0, 11, 3…
## $ per_comaluno_NA              <int> 6, 14, 2, 2, 9, 3, 6, 4, 5, 1, 10, 0, 17,…

Produção e produtividade de artigos

Uma das maneiras de avaliar a produção dos docentes que a CAPES utiliza é quantificando a produção de artigos pelos docentes. Os artigos são categorizados em extratos ordenados (A1 é o mais alto), e separados entre artigos em conferências e periódicos. Usaremos para esse lab a produção em periódicos avaliados com A1, A2 e B1.

cacc = cacc_tudo %>%
  transmute(
    docentes = `Docentes permanentes`,
    producao = (periodicos_A1 + periodicos_A2 + periodicos_B1),
    produtividade = producao / docentes,
    mestrados = Dissertacoes,
    doutorados = Teses,
    tem_doutorado = tolower(`Tem doutorado`) == "sim",
    mestrados_pprof = mestrados / docentes,
    doutorados_pprof = doutorados / docentes,
    nivel = Nível
  )

cacc_md = cacc %>% 
  filter(tem_doutorado)

EDA

skimr::skim(cacc)
Data summary
Name cacc
Number of rows 73
Number of columns 9
_______________________
Column type frequency:
logical 1
numeric 8
________________________
Group variables None

Variable type: logical

skim_variable n_missing complete_rate mean count
tem_doutorado 0 1 0.47 FAL: 39, TRU: 34

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
docentes 0 1 20.63 12.27 8.25 11.25 16.75 25.75 67.25 ▇▃▁▁▁
producao 0 1 58.03 65.44 0.00 18.00 42.00 67.00 355.00 ▇▂▁▁▁
produtividade 0 1 2.36 1.37 0.00 1.40 2.27 3.20 5.66 ▆▇▇▅▂
mestrados 0 1 75.79 63.23 0.00 39.00 58.00 103.00 433.00 ▇▃▁▁▁
doutorados 0 1 14.96 30.98 0.00 0.00 0.00 14.00 152.00 ▇▁▁▁▁
mestrados_pprof 0 1 3.66 1.81 0.00 2.57 3.58 4.88 8.19 ▂▇▇▃▂
doutorados_pprof 0 1 0.43 0.73 0.00 0.00 0.00 0.57 2.69 ▇▁▁▁▁
nivel 0 1 3.84 1.17 3.00 3.00 3.00 4.00 7.00 ▇▅▁▁▁ 
cacc %>% 
  ggplot(aes(x = docentes)) + 
  geom_histogram(bins = 15, fill = paleta[1])

cacc %>% 
  ggplot(aes(x = producao)) + 
  geom_histogram(bins = 15, fill = paleta[2])

cacc %>% 
  ggplot(aes(x = produtividade)) + 
  geom_histogram(bins = 15, fill = paleta[3])

Como produtividade é uma razão (producao/docentes), vale checar se ela tá muito enviesada ou se tem outlier chamativo antes de sair usando ela como variável resposta

# skewness calculada na mão, sem depender de pacote externo
skewness_manual <- function(x) {
    xbar <- mean(x)
    m3 <- mean((x - xbar)^3)
    dp_pop <- sqrt(mean((x - xbar)^2))
    m3 / dp_pop^3
}

skewness_manual(cacc$produtividade)
## [1] 0.4331583
shapiro.test(cacc$produtividade)
## 
##  Shapiro-Wilk normality test
## 
## data:  cacc$produtividade
## W = 0.97371, p-value = 0.1305

A skewness deu 0,43 (levinha assimetria à direita, nada absurdo) e o teste de Shapiro-Wilk não rejeitou normalidade (p = 0,13, bem acima de 0,05). Então produtividade não está gritantemente fora do normal, dá pra seguir usando ela direto num modelo linear sem se preocupar com transformação.

Também testei se tem outlier isolado usando os dois testes que vimos em aula, teste de Grubbs e teste de Dixon Q, implementados na mão (não achei pacote pronto no ambiente que eu tenho acesso aqui):

# Teste de Grubbs (checa se o ponto mais distante da média é um outlier)
grubbs_test <- function(x, alpha = 0.05) {
  n <- length(x)
  xbar <- mean(x)
  s <- sd(x)
  devs <- abs(x - xbar)
  G <- max(devs) / s
  suspeito <- x[which.max(devs)]
  t_crit <- qt(1 - alpha / (2 * n), n - 2)
  G_crit <- ((n - 1) / sqrt(n)) * sqrt(t_crit^2 / (n - 2 + t_crit^2))
  tibble(G = G, G_critico = G_crit, valor_suspeito = suspeito, e_outlier = G > G_crit)
}

grubbs_test(cacc$produtividade)
# Teste de Dixon Q (checa o ponto mais extremo de cada ponta, olhando o "gap" em relação ao range)
dixon_q_test <- function(x) {
  xs <- sort(x)
  n <- length(xs)
  amplitude <- xs[n] - xs[1]
  tibble(
    ponta = c("menor valor", "maior valor"),
    valor_suspeito = c(xs[1], xs[n]),
    Q = c((xs[2] - xs[1]) / amplitude, (xs[n] - xs[n - 1]) / amplitude)
  )
}

dixon_q_test(cacc$produtividade)

Pelo Grubbs, o ponto mais distante da média (produtividade = 5,66) tem G = 2,40, abaixo do G crítico de 3,27, então não é considerado outlier estatisticamente significativo. O Dixon Q deu valores bem baixos nas duas pontas (0,02 e 0,016), bem longe de indicar outlier também — mas vale um adendo: o teste de Dixon Q foi criado pra amostras pequenas (a tabela clássica de valores críticos só cobre até uns 30 pontos), e aqui produtividade tem N = 73. Então trato o resultado do Dixon Q aqui como só mais um indício descritivo, e confio mais no Grubbs (que tem fórmula válida pra qualquer N) como teste formal. De qualquer forma, os dois convergem pra mesma conclusão: não parece ter outlier isolado distorcendo a variável antes de eu montar o modelo.

Se quisermos modelar o efeito do tamanho do programa em termos de docentes (permanentes) na quantidade de artigos publicados, podemos usar regressão.

Importante: sempre queremos ver os dados antes de fazermos qualquer modelo ou sumário:

cacc %>% 
  ggplot(aes(x = docentes, y = producao)) + 
  geom_point()

Parece que existe uma relação. Vamos criar um modelo então:

modelo1 = lm(producao ~ docentes, data = cacc)

tidy(modelo1, conf.int = TRUE, conf.level = 0.95)
glance(modelo1)

Para visualizar o modelo:

cacc_augmented = cacc %>% 
  add_predictions(modelo1) 

cacc_augmented %>% 
  ggplot(aes(x = docentes)) + 
  geom_line(aes(y = pred), colour = "brown") + 
  geom_point(aes(y = producao)) + 
  labs(y = "Produção do programa")

Se considerarmos que temos apenas uma amostra de todos os programas de pós em CC no Brasil, o que podemos inferir a partir desse modelo sobre a relação entre número de docentes permanentes e produção de artigos em programas de pós?

Normalmente reportaríamos o resultado da seguinte maneira, substituindo VarIndepX e todos os x’s e y’s pelos nomes e valores de fato:

Regressão múltipla foi utilizada para analisar se VarIndep1 e VarIndep2 tem uma associação significativa com VarDep. Os resultados da regressão indicam que um modelo com os 2 preditores no formato VarDep = XXX.VarIndep1 + YYY.VarIndep2 explicam XX,XX% da variância da variável de resposta (R2 = XX,XX). VarIndep1, medida como/em [unidade ou o que é o 0 e o que é 1] tem uma relação significativa com o erro (b = [yy,yy; zz,zz], IC com 95%), assim como VarIndep2 medida como [unidade ou o que é o 0 e o que é 1] (b = [yy,yy; zz,zz], IC com 95%). O aumento de 1 unidade de VarIndep1 produz uma mudança de xxx em VarDep, enquanto um aumento…

Produza aqui a sua versão desse texto, portanto:

Regressão linear simples foi utilizada pra analisar se o número de docentes permanentes tem uma associação significativa com a produção de artigos (periódicos A1, A2 e B1) de um programa de pós em CC. Os resultados da regressão indicam que um modelo com o preditor no formato producao = -41,3 + 4,81 x docentes explica 81,5% da variância da produção (R² = 0,815). Docentes, medido como a quantidade de docentes permanentes do programa, tem uma relação significativa com a produção (b = [4,27; 5,36], IC com 95%). O aumento de 1 docente permanente produz uma mudança média de 4,81 artigos na produção do programa, mantendo tudo mais constante.

Dito isso, o que significa a relação que você encontrou na prática para entendermos os programas de pós graduação no Brasil? E algum palpite de por que a relação que encontramos é forte?

Basicamente, quanto mais docente permanente um programa tem, mais artigo ele produz no total, o que não chega a ser surpreendente: é quase um efeito mecânico de soma, mais gente pesquisando tende a gerar mais artigo no total. É por isso que essa métrica (produção total) não é muito boa pra comparar programas de tamanhos diferentes, um programa de 30 docentes vai produzir mais artigo que um de 5 só pelo tamanho, não necessariamente porque é “melhor” ou mais eficiente. Pra comparar de fato a eficiência dos programas entre si, faz mais sentido olhar produtividade (produção dividida pelo número de docentes), que é o que fazemos mais pra frente no lab.

Mais fatores

modelo2 = lm(producao ~ docentes + mestrados_pprof + doutorados_pprof + tem_doutorado, 
             data = cacc_md)

tidy(modelo2, conf.int = TRUE, conf.level = 0.95)
glance(modelo2)

Reparei que tem_doutorado sai com todos os valores NA nesse modelo. Isso não é erro de conta, é porque cacc_md já é o subconjunto filtrado só com os programas que têm doutorado (filter(tem_doutorado)), então dentro desse subconjunto a variável tem_doutorado é sempre TRUE, não varia. Sem variação não dá pra estimar coeficiente nenhum pra ela. Então nesse modelo só valem de fato os coeficientes de docentes, mestrados_pprof e doutorados_pprof.

E se considerarmos também o número de alunos?

modelo2 = lm(producao ~ docentes + mestrados + doutorados, data = cacc)

tidy(modelo2, conf.int = TRUE, conf.level = 0.95)
glance(modelo2)

Visualizar o modelo com muitas variáveis independentes fica mais difícil

para_plotar_modelo = cacc %>% 
  data_grid(producao = seq_range(producao, 10), # Crie um vetor de 10 valores no range
            docentes = seq_range(docentes, 4),  
            # mestrados = seq_range(mestrados, 3),
            mestrados = median(mestrados),
            doutorados = seq_range(doutorados, 3)) %>% 
  add_predictions(modelo2)

glimpse(para_plotar_modelo)
## Rows: 120
## Columns: 5
## $ producao   <dbl> 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.000…
## $ docentes   <dbl> 8.25000, 8.25000, 8.25000, 27.91667, 27.91667, 27.91667, 47…
## $ mestrados  <int> 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58,…
## $ doutorados <dbl> 0, 76, 152, 0, 76, 152, 0, 76, 152, 0, 76, 152, 0, 76, 152,…
## $ pred       <dbl> 3.199123, 79.257725, 155.316327, 72.026777, 148.085378, 224…
para_plotar_modelo %>% 
  ggplot(aes(x = docentes, y = pred)) + 
  geom_line(aes(group = doutorados, colour = doutorados)) + 
  geom_point(data = cacc, aes(y = producao, colour = doutorados))

Considerando agora esses três fatores, o que podemos dizer sobre como cada um deles se relaciona com a produção de um programa de pós em CC? E sobre o modelo? Ele explica mais que o modelo 1?

EXPLICAÇÃO: Regressão múltipla foi utilizada pra analisar se docentes permanentes, dissertações de mestrado e teses de doutorado têm associação significativa com a produção de artigos. Os resultados indicam que um modelo com os três preditores no formato producao = -14,4 + 3,50 x docentes - 0,195 x mestrados + 1,00 x doutorados explica 87,1% da variância da produção (R² = 0,871, R² ajustado = 0,865), mais que o modelo 1, que só com docentes já explicava 81,5%. Docentes (b = [2,58; 4,42], IC 95%) e doutorados (b = [0,636; 1,37], IC 95%) têm relação positiva e significativa com a produção. Já mestrados tem uma relação negativa e significativa (b = [-0,358; -0,032], IC 95%). Ou seja, mantendo os outros fatores constantes: cada docente permanente a mais está associado a 3,50 artigos a mais; cada tese de doutorado a mais está associado a 1,00 artigo a mais, mas cada dissertação de mestrado a mais está associada a uma queda média de 0,195 artigo.

Esse coeficiente negativo de mestrados chama minha atenção, dando a entender que isso reflita um trade-off de uso do tempo do docente: orientar muita dissertação de mestrado consome tempo e esforço que poderiam estar sendo usados pra escrever artigo ou orientar doutorado (que tende a gerar publicação de mais impacto). Não necessariamente quer dizer que “mestrado atrapalha o programa”, mas sim que programas mais focados/pesados em mestrado (proporcionalmente) tendem a produzir menos artigo do que programas mais focados em doutorado, quando o número de docentes já está controlado no modelo.

Agora produtividade

Diferente de medirmos produção (total produzido), é medirmos produtividade (produzido / utilizado). Abaixo focaremos nessa análise. Para isso crie um modelo que investiga como um conjunto de fatores que você julga que são relevantes se relacionam com a produtividade dos programas. Crie um modelo que avalie como pelo menos 3 fatores se relacionam com a produtividade de um programa. Pode reutilizar fatores que já definimos e analizamos para produção. Mas cuidado para não incluir fatores que sejam função linear de outros já incluídos (ex: incluir A, B e um tercero C=A+B)

Produza abaixo o modelo e um texto que comente (i) o modelo, tal como os que fizemos antes, e (ii) as implicações - o que aprendemos sobre como funcionam programas de pós no brasil?.

Escolha dos fatores

Os três fatores que escolhi pra explicar produtividade (produção de artigos por docente permanente) foi:

  • tem_doutorado: o programa tem doutorado ou não (0/1)
  • mestrados_pprof: dissertações de mestrado por docente permanente
  • doutorados_pprof: teses de doutorado por docente permanente

Eu não incluí docentes como preditor aqui, de propósito, como produtividade = producao / docentes, usar docentes puro como preditor de uma variável que já divide por docentes pode gerar uma correlação meio artificial (correlação espúria de razão), então preferi ficar só com fatores que não têm essa relação direta de construção com a variável resposta.

Também segui a dica do enunciado: doutorados_pprof é 0 pra quase metade dos programas (os que não têm doutorado), e isso não é dado faltante, é zero de verdade (não existe doutorando pra contar). Por isso incluí tem_doutorado no modelo junto com doutorados_pprof, assim o modelo consegue separar “o efeito de ter doutorado, ponto” do “efeito de quanto o programa produz de doutorado proporcionalmente”, em vez de confundir os dois.

cacc %>% 
  select(produtividade, tem_doutorado, mestrados_pprof, doutorados_pprof) %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 73
Number of columns 4
_______________________
Column type frequency:
logical 1
numeric 3
________________________
Group variables None

Variable type: logical

skim_variable n_missing complete_rate mean count
tem_doutorado 0 1 0.47 FAL: 39, TRU: 34

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
produtividade 0 1 2.36 1.37 0 1.40 2.27 3.20 5.66 ▆▇▇▅▂
mestrados_pprof 0 1 3.66 1.81 0 2.57 3.58 4.88 8.19 ▂▇▇▃▂
doutorados_pprof 0 1 0.43 0.73 0 0.00 0.00 0.57 2.69 ▇▁▁▁▁ 
cacc %>% 
  ggplot(aes(x = tem_doutorado, y = produtividade, fill = tem_doutorado)) + 
  geom_boxplot(show.legend = FALSE) + 
  labs(x = "Tem doutorado", y = "Produtividade (artigos por docente)")

cacc %>% 
  ggplot(aes(x = doutorados_pprof, y = produtividade, colour = tem_doutorado)) + 
  geom_point() + 
  labs(x = "Doutorados por docente", y = "Produtividade (artigos por docente)", colour = "Tem doutorado")

Dá pra ver nos gráficos que quem tem doutorado já parte de uma produtividade mais alta, e que mais teses de doutorado por docente parece vir junto com mais produtividade também.

O modelo

modelo_produtividade = lm(produtividade ~ tem_doutorado + mestrados_pprof + doutorados_pprof, 
                           data = cacc)

tidy(modelo_produtividade, conf.int = TRUE, conf.level = 0.95)
glance(modelo_produtividade)
para_plotar_produtividade = cacc %>% 
  data_grid(tem_doutorado, 
            mestrados_pprof = median(mestrados_pprof),
            doutorados_pprof = seq_range(doutorados_pprof, 10)) %>% 
  add_predictions(modelo_produtividade)

para_plotar_produtividade %>% 
  ggplot(aes(x = doutorados_pprof, y = pred, colour = tem_doutorado)) + 
  geom_line() + 
  geom_point(data = cacc, aes(y = produtividade, colour = tem_doutorado)) + 
  labs(x = "Doutorados por docente", y = "Produtividade prevista", colour = "Tem doutorado")

Antes de confiar no modelo, vale checar se os resíduos dele se comportam mais ou menos como o esperado (aproximadamente normais), do jeito que vimos no material sobre regressão:

modelo_produtividade_aug = augment(modelo_produtividade)

modelo_produtividade_aug %>% 
  ggplot(aes(sample = .std.resid)) + 
  stat_qq(color = paleta[4]) + 
  stat_qq_line(color = paleta[5], linetype = "dashed") + 
  labs(x = "Quantis teóricos (normal)", y = "Resíduos padronizados",
       title = "QQ-plot dos resíduos - modelo_produtividade")

shapiro.test(residuals(modelo_produtividade))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_produtividade)
## W = 0.9666, p-value = 0.04996

Os pontos seguem bem de perto a linha nas duas pontas, com um leve desvio nas caudas. O teste de Shapiro-Wilk deu W = 0,967, p = 0,050 — bem em cima da fronteira do 0,05, então a normalidade dos resíduos é um pouco questionável, mas não claramente violada. É uma limitação honesta do modelo, não motivo pra descartá-lo.

(i) Sobre o modelo: Regressão múltipla foi utilizada pra analisar se ter doutorado, a quantidade de dissertações de mestrado por docente e a quantidade de teses de doutorado por docente têm associação significativa com a produtividade de um programa. Os resultados indicam que um modelo com os três preditores no formato produtividade = 1,47 + 1,10 x tem_doutorado + 0,018 x mestrados_pprof + 0,730 x doutorados_pprof explica 51,3% da variância da produtividade (R² = 0,513, R² ajustado = 0,492). tem_doutorado, medido como 1 = programa tem doutorado e 0 = não tem, tem relação significativa e positiva com a produtividade (b = [0,508; 1,69], IC 95%), assim como doutorados_pprof, medido como número de teses de doutorado por docente permanente (b = [0,323; 1,14], IC 95%). Já mestrados_pprof, medido como número de dissertações de mestrado por docente permanente, não teve relação significativa com a produtividade (b = [-0,109; 0,145], IC 95%, esse intervalo cruza o zero). Ter doutorado está associado a um aumento médio de 1,10 artigo por docente na produtividade do programa (mantendo os outros fatores constantes); cada tese de doutorado a mais por docente está associada a um aumento médio de 0,73 artigo por docente.

Como teste de robustez, testei também incluir o conceito CAPES (nivel) como quarto fator:

modelo_produtividade_nivel = lm(produtividade ~ tem_doutorado + mestrados_pprof + doutorados_pprof + nivel, 
                                 data = cacc)

tidy(modelo_produtividade_nivel, conf.int = TRUE, conf.level = 0.95)
glance(modelo_produtividade_nivel)

Pra verificar isso de um jeito mais rigoroso que só olhar de longe, segui o que o slide de multicolinearidade ensina: olhar a correlação entre as preditoras e ver se o p-valor de algum coeficiente fica esquisito (alto) quando incluímos uma preditora correlacionada:

cor(cacc$doutorados_pprof, cacc$nivel)
## [1] 0.8320627
tidy(modelo_produtividade, conf.int = TRUE) %>% select(term, estimate, p.value)
tidy(modelo_produtividade_nivel, conf.int = TRUE) %>% select(term, estimate, p.value)

A correlação entre doutorados_pprof e nivel é de 0,83, bem alta. E dá pra ver exatamente o “indício” que o slide descreve: o p-valor de doutorados_pprof pula de 0,0006 (bem significativo, no modelo de 3 fatores) pra 0,52 (nada significativo) quando nivel entra no modelo — mesmo o R² geral do modelo tendo melhorado. Isso é sinal de colinearidade: as duas variáveis carregam bastante informação repetida (faz sentido, já que conceito CAPES alto tende a vir de programa com bastante produção de doutorado), então o modelo passa a “brigar” pra decidir qual das duas leva o crédito pelo efeito, e isso infla o erro padrão e derruba a significância de doutorados_pprof. Preferi manter o modelo de 3 fatores como principal por ser mais simples e não ter esse problema, mas achei válido mostrar que incluir o conceito CAPES melhora o ajuste, com essa ressalva.

(ii) Implicações: O que dá pra aprender daqui sobre os programas de pós em CC no Brasil é que ter doutorado parece funcionar como uma espécie de selo de programa mais maduro/consolidado, mesmo controlando pela intensidade de orientação de doutorado (doutorados_pprof), só o fato de ter doutorado já eleva a produtividade média do programa. Isso sugere que não é só “ter aluno de doutorado que produz mais artigo”, tem algo estrutural também, tipo programas com doutorado tendendo a atrair docente mais experiente, ter mais infraestrutura de pesquisa, financiamento melhor, etc, coisas que não estão diretamente nos dados mas que provavelmente colaboram.

Já o mestrado conta uma história diferente, mestrados_pprof não teve relação significativa com produtividade, reforçando o que já tinha aparecido no modelo de produção (coeficiente negativo de mestrados lá atrás). Formar muito mestrando por docente não parece se traduzir em mais artigo publicado, então talvez o esforço de orientação de mestrado compita pelo mesmo tempo do docente que poderia render artigo, sem necessariamente resultar num ganho de produtividade equivalente ao que a orientação de doutorado traz.

De forma geral, dá pra dizer que a política de avaliação da CAPES parece favorecer bastante os programas com doutorado consolidado, e que produtividade em pesquisa (pelo menos medida como artigo por docente) parece estar mais ligada à pós-graduação em nível de doutorado do que de mestrado.

Resumo: comparando todos os modelos

Por fim, um resumão de todos os modelos desse relatório=

bind_rows(
  glance(modelo1) %>% mutate(modelo = "producao ~ docentes"),
  glance(modelo2) %>% mutate(modelo = "producao ~ docentes + mestrados + doutorados"),
  glance(modelo_produtividade) %>% mutate(modelo = "produtividade ~ tem_doutorado + mestrados_pprof + doutorados_pprof"),
  glance(modelo_produtividade_nivel) %>% mutate(modelo = "produtividade ~ tem_doutorado + mestrados_pprof + doutorados_pprof + nivel")
) %>% 
  select(modelo, r.squared, adj.r.squared) %>% 
  mutate(across(where(is.numeric), ~round(., 3)))

Os modelos de producao (bruta) têm R² bem mais alto que os de produtividade, o que faz sentido, produção total é dominada pelo tamanho do programa (quantos docentes ele tem), então é “mais fácil” de prever. Produtividade, por outro lado, já tirou o efeito de tamanho de cima e sobrou uma variável mais difícil de explicar, então R² menor aqui não é sinal de modelo ruim, é só um problema mais difícil. Entre os dois modelos de produtividade, o que inclui nivel tem R² ajustado um pouco maior (0,54 vs 0,49), mas como vimos, isso vem com o custo da colinearidade com doutorados_pprof, então a escolha de qual usar depende do que se valoriza mais: um pouco mais de poder explicativo (com colinearidade) ou um modelo mais simples e estável (os 3 fatores originais).