Relatório.Rmd

---
title: "Análise de evasão dos cursos de Psicologia no Brasil"
author: "Alisson Rosa Pereira, Caroline Cogo Carneosso, Vítor Bernardo Silveira Pereira"
date: "08/07/2021"
header-includes:
   - \usepackage[brazil]{babel}
geometry: margin=2cm
output:
  bookdown::pdf_document2:
editor_options:
  chunk_output_type: console
---

```{r setup,,include=F,}

library('tidyverse')
library('janitor')
library('kableExtra')
library('cowplot')
library('tidymodels')
library('epiR')
library('DescTools')
knitr::opts_chunk$set(echo=FALSE,message=F,warning=F,fig.pos = 'H',fig.align = 'center',fig.width=7.8, fig.height=4.85)
options(digits=3)
options(OutDec=",")
theme_set(theme_minimal())
scale_fill_discrete = \(...) scale_fill_brewer(... , palette="Set2")



```


```{r,message=F,echo=F,warning=F}
df=read_csv("https://raw.githubusercontent.com/AlissonRP/Categorical--Data/master/ALUNO_PSICO2019.csv")
df = df %>% 
  filter(NU_ANO_INGRESSO==2019,TP_SITUACAO!="7") %>%
      mutate(NU_IDADE=as.integer(NU_IDADE), CO_IES= as.integer(CO_IES),
           Idade=cut(NU_IDADE,breaks=c(15,20,25,30,65,90),
                     include.lowest = TRUE)) %>% 
      mutate(NU_IDADE=as.character(NU_IDADE)) %>% 
      mutate_if(is.double,as.factor) %>% 
      mutate(NU_IDADE=as.double(NU_IDADE),id=as.double(NU_IDADE)) %>%
      select(Turno=TP_TURNO,Sexo=TP_SEXO,Nacionalidade=TP_NACIONALIDADE,
             Situação=TP_SITUACAO,`Apoio social`=IN_APOIO_SOCIAL,
             id, Idade, CO_IES, CO_UF_NASCIMENTO) %>% 
      mutate(Sexo=fct_recode(Sexo,"Feminino"="1","Masculino"="2"),
             Turno=fct_recode(Turno,"Matutino"="1","Vespertino"="2",
                              "Noturno"="3","Integral"="4"),
            Nacionalidade=fct_recode(Nacionalidade,"Brasileira"="1",                              "Exterior/Naturalizado"="2","Estrangeira"="3"),
            Situação=fct_collapse(Situação,Evadido=c("5","4"),
                                  Retido=c("2","3","6")),
            `Apoio social`=fct_recode(`Apoio social`,"Não"= "0","Sim"="1"))

IES <- read_delim('IES.CSV',delim=",")
IES = IES  %>% 
  rename(Codigo = CO_IES)
df <- df %>% 
  inner_join(IES, by = c('CO_IES' = 'Codigo')) %>%
  mutate(Mobilidade = ifelse(CO_UF_NASCIMENTO == CO_UF, 
                               "Mesmo estado", "Estado diferentes")) %>%
  select(-X1, -CO_UF_NASCIMENTO, -CO_IES, -CO_REGIAO, -CO_UF)
```



```{r,tabelas de conting,include=F}
#Função para gerar as tabelas de frequência para uma variavel
freq_table=function(df,v,tit){
  df %>% 
  count({{v}}) %>% 
  mutate(prop=prop.table(n) %>% round(3)) %>%
  rename(`Frequência Absoluta`=n,`Frequência Relativa`=prop) %>% 
  adorn_totals("row") %>% 
  kable(caption=tit,align = "c") %>% 
  kable_classic(latex_options = "HOLD_position") %>% 
  footnote(general = "Elaborado pelos autores ",
         general_title = "Fonte:",
         footnote_as_chunk = T, title_format = c("italic"))
}


freq_table2=function(df,v1,v2,tit,marg){
  if(missing(marg)){
    df %>% 
      group_by({{v1}},{{v2}}) %>%  
      summarise(n=n())   %>%  
      spread({{v2}}, n) %>% 
      adorn_totals("row") %>%  
      adorn_totals("col") %>% 
      ungroup() %>% 
      kable(caption=tit,align = "c") %>%
      footnote(general = "Elaborado pelos autores ",
         general_title = "Fonte:",
         footnote_as_chunk = T, title_format = c("italic")) %>% 
      kable_minimal(latex_options = "HOLD_position") %>% 
      add_header_above(c(" ","Situação"=3),align ="c")
 
  } else {
    df %>% 
     group_by({{v1}},{{v2}}) %>%  
     summarise(n=n())   %>%  
     spread({{v2}}, n) %>% 
     adorn_percentages() %>% 
     ungroup() %>% 
     mutate(Retido=Retido %>% round(3),Evadido=Evadido%>% round(3)) %>% 
     kable(caption=tit,align = "c") %>% 
    footnote(general = "Elaborado pelos autores ",
         general_title = "Fonte:",
         footnote_as_chunk = T, title_format = c("italic")) %>% 
     kable_minimal(latex_options = "HOLD_position") %>% 
     add_header_above(c(" ","Situação"=2),align ="c")}
  }

conf_table=function(fit,p,tit){
  fit %>% 
   p(p) %>% 
 select(Predito=.pred_class,Observado=sit) %>%
    group_by(Predito,Observado) %>% 
     summarise(n=n())   %>%  
      spread(Observado, n) %>% 
      adorn_totals("col") %>%  
  kable(caption=tit,align = "c") %>%
      footnote(general = "Elaborado pelos autores ",
         general_title = "Fonte:",
         footnote_as_chunk = T, title_format = c("italic")) %>% 
      kable_minimal(latex_options = "HOLD_position") %>% 
    add_header_above(c(" ","Observado"=3),align ="c")
 
}

```
```{r, associacao,include=F}
Associar <-function(df,v1,tipo){
  ## tipo = o -> v1 é ordinal
  ## tipo = nd -> v1 é nominal dicotômica
  ## tipo = np -> v1 é nominal politômica
  df <- drop_na(df, {{v1}})
  if(tipo == "nd"){
    Qp <- df %>%
    select({{v1}},Situação) %>%
    table() %>%
    stats::chisq.test(correct = F)
    Asso <- df %>%
      select({{v1}},Situação) %>%
      table() %>%
      epi.2by2(method = "cross.sectional")
    lista <- list('Est' = Qp$statistic, 'p' = Qp$p.value,
                  'RPrevalencia' = Asso$massoc.summary$est[1], 
                  'RChance' = Asso$massoc.summary$est[2], 
                  'ICPrevalencia' = c(Asso$massoc.summary$lower[1],
                                      Asso$massoc.summary$upper[1]),
                  'ICChance' = c(Asso$massoc.summary$lower[2],
                                      Asso$massoc.summary$upper[2]))
    return(lista)
  }
  
  if(tipo == "np"){
    Qp <- df %>%
    select({{v1}},Situação) %>%
    table() %>%
    stats::chisq.test(correct = F)
    lista <- list('Est' = Qp$statistic, 'p' = Qp$p.value)
    return(lista)
  }
  
  if(tipo == "o"){
    y <- df %>% 
    mutate(Situação=fct_recode(Situação,"0"="Retido","1"="Evadido")) %>%
      select(Situação)
    y1 <- ifelse(y == '1',1,0)
    x <- df %>%
      mutate(Idade=fct_recode(Idade,'1'="[15,20]", '2' = "(20,25]", '3' = "(25,30]", 
                              '4' = "(30,65]", '5' = "(65,90]")) %>%
      select(Idade)
    x1 <- ifelse(x == '1',1, ifelse(x == '2', 2, ifelse(x == '3', 3, 
                 ifelse(x == '4', 4,ifelse(x == '5',5,6)))))
    rac<-cor(x1,y1)
    n1<- length(y)
    qcs<-(n1-1)*rac^2
    p<-1-pchisq(qcs,1)  
    lista <- list('Est' = qcs, 'p' = p, 'cor' = rac)
    return(lista)
  }
    
}

#Enviar para a função o banco (df), v1 é o X, vetor são os vetores de variáveis que queremos testar como interferente
Inter <- function(df,v1,vetor){
  #verificando se v1 tem NA
  df <- drop_na(df, {{v1}})
  #Criando um data frame
  a <- NULL
  a <- as.data.frame(a)
  for(i in vetor){
    #testando se as variaveis do vetor tem MA
    df1 <- subset(df, !is.na(eval(parse(text=i))))
    #descobrindo o número de fatores da variavel do vetor
    d1 <- df1 %>% 
      group_by(eval(parse(text=i))) %>%
      summarise(no_rows = length(eval(parse(text=i)))) %>%
      dim()%>%
      .[1]
    #nome da variavel
    n <- df1 %>% 
          select({{v1}}) %>%
          colnames()
    #Fazendo o teste de BreslowDayTest
    BDT <- table(eval(parse(text=paste('df1$',n, sep = ''))), df1$Situação, 
          eval(parse(text=paste('df1$',i, sep = '')))) %>%
          array(dim = c(2,2,d1)) %>%
          BreslowDayTest()
    j = 1
    #Setando os valores no DF
    a[j,2] = BDT$p.value
    a[j,1] = i
    a[j,3] = BDT$statistic
    #Resetando as variaveis
    BDT <- NULL
    n <- NULL
    d1 <- NULL
    df1 <- NULL
    j = j + 1
  }
  colnames(a) <- c('Var. Inter.', 'P', 'Est')
  #retornando o dataframe
  return(a)
}
```

```{r teste_assoc, include = F}
#th sit vs mob 
#nominal dicotomica
a<- Associar(df,Mobilidade,'nd')

#th sit vs idade 
#politomica ordinal
b<-Associar(df,Idade,'o')

#th sit vs turno
#nominal politomica
c<-Associar(df,Turno,'np')

#th sit vs sexo
d<-Associar(df,Sexo,'nd')

#th sit vs Nacionalidade
e<-Associar(df,Nacionalidade,'np')

#th sit vs Apoio Social
f<-Associar(df,`Apoio social`,'nd')


```

\section{Introdução}

O Instituto Nacional de Estudos e Pesquisas Educacionais Anísio Teixeira (INEP) disponibiliza o Censo [^1] da educação superior, o qual fornece informações detalhadas sobre os cursos superiores no Brasil. O objetivo do trabalho é estudar como as variáveis selecionadas influenciam na váriavel de desfecho. Trata-se de um estudo observacional transversal, uma vez que os ingressantes são estudados em um momento específico do tempo, neste caso no ano de 2019.
Neste relatório faremos uma breve análise sobre o curso de Psicologia das Instituições de Ensino Superior (IES), estudando a permanência dos estudantes que ingressaram no ano de 2019, utiliza-se como variáveis de estudo:

Como variável de interesse tem-se a **Situação** do estudante, após o primeiro ano matriculado no ensino superior, ela é classificada em Retido e Evadido.

**Turno**: Uma variável categórica nominal que indica o tipo de turno que o estudante está vinculado, sendo dividido em 4 categorias.

**Sexo**: Uma variável categórica binária  que indica se o estudante é do sexo masculino ou feminino. 
 
**Nacionalidade**: Uma variável categórica nominal  que indica se o aluno é Brasileiro, Exterior/Naturalizado ou Estrangeiro. 

**Apoio Social** Uma variável categórica dicotômica que informa se o estudante possui ou não apoio social.

**Idade** Uma variável que foi discretizada em intervalos, isto é tornando-se uma variável categórica ordinal

E por último a variável **Mobilidade** que foi criada pela junção do banco de dados dos alunos com o banco das Instituições de Ensino Superior, sendo esta uma variável categórica dicotômica para verificar se os alunos são do mesmo estado da universidade ou não.

Adota-se como convenção a partir de agora que todas as observações referem-se ao ano 2019. Todas as tabelas e gráficos foram elaboradas pelos autores com base nas `r nrow(df)` observações do banco de dados, foram selecionadas `r ncol(df)-1` variáveis .


```{r ead graphs,warning=FALSE,message=FALSE}


#função para gerar os graficos de barras
g_t=function(df,v){
  df %>% 
  ggplot(aes({{v}})) + 
  geom_bar(aes(fill={{v}}))+
  geom_text(stat='count',aes(label=..count..),vjust=-0.5)+
  theme_minimal()+
    labs(y="Quantidade",x=df %>%
         select({{v}}) %>% 
         names(),caption = "Fonte: Elaborado pelos autores")+
  scale_fill_brewer(palette="Set2")
    
}

#Gráfico de barras duas variáveis
g_t2=function(df,v1,v2){
df %>% 
  ggplot(aes({{v1}},fill={{v2}})) + 
  geom_bar(position = "dodge")+
  scale_fill_brewer(palette="Set2")+
  theme_minimal()+
  labs(y="Quantidade",x=df %>%
         select({{v1}}) %>% 
         names(),caption = "Fonte: Elaborado pelos autores")
}

```

\section{Estudo Descritivo}
A partir dos dados, foram desenvolvidas tabelas e gráficos de frequência, para melhor compreensão e interpretação das variáveis a serem estudadas:

**Situação**: 
```{r freqsit}
freq_table(df,Situação,"Frequência absoluta e relativa para a variável Situação.")
  
```
```{r barsit, fig.cap = "Gráfico de barras com os valores absolutos da variável Situação."}
g_t(df,Situação)+
  theme(legend.position="none")
```
A partir da tabela \@ref(tab:freqsit) e do gráfico \@ref(fig:barsit), verifica-se que cerca de $82,26\%$ foram retidos, ou seja continuaram no curso de Psicologia, já $17,74\%$ foram evadidos do curso.

**Turno**: 
```{r freqturno}
freq_table(df,Turno,"Frequência absoluta e relativa para a variável Turno.")
  
```
```{r barturno, fig.cap = "Gráfico de barras com os valores absolutos da variável Turno"}
g_t(df,Turno)+
  theme(legend.position="none")
```
Observa-se na  tabela \@ref(tab:freqturno) e no gráfico \@ref(fig:barturno) que aproximadamente $60\%$ dos estudantes são do turno noturno, apenas $2\%$ são do turno vespertino.

**Sexo**: 
```{r freqsexo}
freq_table(df,Sexo,"Frequência absoluta e relativa para a variável Sexo.")  
  
```

```{r  barsex,fig.cap = "Gráfico de barras com os valores absolutos da variável Sexo."}
g_t(df,Sexo)+
  theme(legend.position="none")
```
Como é observado na tabela \@ref(tab:freqsexo) e no gráfico \@ref(fig:barsex), mais de $75\%$ dos alunos ingressantes são do sexo feminino.

**Nacionalidade**: 
```{r freqnac}
freq_table(df,Nacionalidade, "Frequência absoluta e relativa para a variável Nacionalidade.") 
  
```
Ao analisar a tabela \@ref(tab:freqnac) é possivel constatar que a maioria dos ingressantes é brasileira, cerca de $99,7\%$.

**Apoio Social** 

```{r freqapoio}
freq_table(df,`Apoio social`,"Frequência absoluta e relativa para a variável Apoio Social.")
```
Percebe-se pela tabela \@ref(tab:freqapoio) que aproximadamente $97\%$ dos estudantes não recebem apoio social.
 
A tabela \@ref(tab:freqidade) e o gráfico \@ref(fig:baridade) representam a  **idade** discretizada, ou seja a quantidade de observações nos respectivos intervalos para cada idade.
```{r freqidade}
freq_table(df,Idade, "Frequência absoluta e relativa para a variável Idade dividida em classes.")
```

```{r  baridade,fig.cap = "Gráfico de barras com os valores absolutos da variável Idade dividida em classes."}
g_t(df,Idade)+
  theme(legend.position="none")
```
Como existe uma grande concentração de ingressantes na faixa etária entre 15 a 20 anos, cerca de $41\%$, é construído o gráfico abaixo, o qual se refere a quantidade de alunos com a respectiva idade.

```{r freq1520,fig.cap = "Valores Absolutos da variável Idade na faixa etaria 15 a 20 anos"}
df %>% 
  filter(id %in% c(15:20)) %>% 
  select(id) %>%
  mutate(id=as.factor(id)) %>% 
  g_t(id) +
  geom_curve(aes(x = 2, y = 4000, xend = 1, yend = 700),
  arrow = arrow(
  length = unit(0.03, "npc"),type="open"))+
  scale_x_discrete(breaks=seq(15,22,1))+
  theme(legend.position="none")+
  labs(x="Idade")+
  draw_image(
  "https://images.emojiterra.com/twitter/v13.0/512px/1f633.png",
  x = 1.62, y = 3750, width = 0.6, height = 1200)

```

**Mobilidade** 
```{r freqmob}
df %>% 
  subset(!is.na(Mobilidade)) %>%
freq_table(Mobilidade,"Frequência absoluta e relativa para a variável Mobilidade.")
```

```{r  barmob,fig.cap = "Gráfico de barras com os valores absolutos da variável Mobilidade."}
df %>% 
  subset(!is.na(Mobilidade)) %>%
g_t(Mobilidade)+
  theme(legend.position="none")
```
Nota-se que a partir da tabela \@ref(tab:freqmob) e da figura \@ref(fig:barmob), cerca de $82\%$ dos alunos ingressaram na universidade do seu estado de origem.

\section{Estudo de Associação}
A partir do estudo de frequência realizado pelo Estudo Descritivo, iremos construir tabelas, gráficos e testes para verificar a associação entre a variável de desfecho e as variáveis de interesse, com as seguintes hipóteses: 
$$H_0: \mbox{Ausência de associação entre as variáveis}$$ versus a hipótese alternativa $$H_1: \mbox{Presença de associação}$$.

\subsection{Tabelas de Contigência e Testes}
Nesta seção, iremos verificar a relação entre as variáveis com a visualização das tabelas de contingência, gráfico e a análise dos testes qui-quadrado.


**Turno vs Situação** 
```{r freqts,message=F,warning=F}
freq_table2(df,Turno,Situação,"Tabela de contingência com os valores absolutos da variável Turno, considerando a variável de desfecho.")
```
```{r grapht,fig.cap = "Gráfico de barras para o turno dos ingressantes conforme a variável de desfecho."}
g_t2(df,Turno,Situação)
```
A  tabela \@ref(tab:freqts) e o gráfico \@ref(fig:grapht) expõem  que há prevalência dos alunos no curso de Psicologia em todos os turnos, sendo o noturno o turno com maior evasão em valores absolutos.

Para o teste de associação, temos que a estatística de teste é igual a `r round(as.numeric(c$Est),3)`, com p-valor menor que 0.01, ou seja, rejeita-se H0 com todos os níveis de significância, assim há evidência de que há associação entre as variáveis.

**Sexo vs Situação**
```{r freqss,message=F,warning=F}
freq_table2(df,Sexo,Situação,"Tabela de contingência com os valores absolutos da variável Sexo, considerando a variável de desfecho.")
```

```{r graphsex,fig.cap = "Gráfico de barras para o Sexo dos ingressantes de acordo com a variável de desfecho."}
g_t2(df,Sexo,Situação)
```
A  tabela \@ref(tab:freqss) e o gráfico \@ref(fig:graphsex) evidenciam em valores absolutos que a maior parte dos ingressantes permanecem no curso em ambos os sexos.

Para o teste de associação, temos que a estatística de teste é igual a `r round((as.numeric(d$Est)),3)`, com p-valor menor que 0.01 , sendo assim, rejeitamos H0 a todos os níveis de significância, logo há evidência de que há associação entre as variáveis. Sendo uma tabela 2x2 podemos calcular a razão de prevalências que é: `r round(as.numeric(d$RPrevalencia),3)`, com o intervalo de confiança  [`r round(as.numeric(d$ICPrevalencia)[1],3)`;`r round(as.numeric(d$ICPrevalencia)[2],3)`], então podemos ver que 1 não está contido no intervalo, reforçando esta associação. O risco de retençao para alunos do sexo feminino é de 1,02 maior que dos alunos do sexo masculino.

**Nacionalidade vs Situação**
```{r freqsnac}
freq_table2(df,Nacionalidade,Situação,"Tabela de contingência com os valores absolutos da variável Nacionalidade, considerando a variável de desfecho.")
```

```{r graphnac,fig.cap = "Gráfico de barras para o Nacionalidade dos ingressantes de acordo com a variável de desfecho."}
g_t2(df,Nacionalidade,Situação)
```
A  tabela \@ref(tab:freqsnac) e o gráfico \@ref(fig:graphnac) apresentam em valores absolutos que a minoria dos alunos não-brasileiros evadem.

Para o teste de associação, temos que a estatística de teste é igual a `r round(as.numeric(e$Est),3)`, com p-valor menor que 0.01, ou seja, rejeita-se H0 com todos os níveis de significância, assim há evidência de que há associação entre as variáveis.

**Apoio Social vs Situação**
```{r freqsapoio}
freq_table2(df,`Apoio social`,Situação,"Tabela de contingência com os valores absolutos da variável Apoio social, considerando a variável de desfecho.")
```
```{r graphas,fig.cap = "Gráfico de barras para o Apoio Social dos ingressantes de acordo com a variável de desfecho."}
g_t2(df,`Apoio social`,Situação)

```
A  tabela \@ref(tab:freqsapoio) e o gráfico \@ref(fig:graphas) nos informa que evidentemente os estudantes que não possuem apoio social possuem maior retenção em valor absoluto (são $97\%$ dos ingressantes), mas veremos na próxima seção em termos relativos quanto é essa retenção.

Para o teste de associação, temos que a estatística de teste é igual a `r round((as.numeric(f$Est)),3)`, com p-valor menor que 0.01 , sendo assim, rejeitamos H0 a todos os níveis de significância, logo há evidência de que há associação entre as variáveis. Sendo uma tabela 2x2 podemos calcular a razão de prevalências que é: `r round(as.numeric(f$RPrevalencia),3)`,  com o  intervalo de confiança  [`r round(as.numeric(f$ICPrevalencia)[1],3)`;`r round(as.numeric(f$ICPrevalencia)[2],3)`], então podemos ver que 1 não está contido no intervalo, reforçando esta associação. O risco de retenção para os alunos de Psicologia que não possuem Apoio Social é em torno de 0,15 menor dos que possuem Apoio Social.

**Idade vs Situação**
```{r}
freq_table2(df,Idade,Situação,"Tabela de contingência com os valores absolutos da variável Idade, considerando a variável de desfecho.")
```
```{r graphidade,fig.cap = "Gráfico de barras para o Idade dos ingressantes de acordo com a variável de desfecho."}
g_t2(df,Idade,Situação)
```
Vemos que a faixa etária 15 a 20 anos possui maior retenção em valor absoluto, entretanto possui a maior quantidade de ingressantes, na proxima seção constataremos essa informação de forma mais detalhada.

Para o teste de associação, temos que a estatística de teste é igual a `r as.numeric(b$Est)` com o p valor igual a `r as.numeric(b$p)`, ou seja, não há nenhuma evidência de associação entre Idade e Situação. A correlação é igual a `r round(as.numeric(b$cor),3)`, então a relação linear entre as duas variáveis é muito baixa.


**Mobilidade vs Situação**
```{r freqsmob}
df %>% 
  subset(!is.na(Mobilidade)) %>%
freq_table2(Mobilidade,Situação,"Tabela de contingência com os valores absolutos da variável Mobilidade, considerando a variável de desfecho.")
```

```{r graphmob,fig.cap = "Gráfico de barras para o Mobilidade dos ingressantes de acordo com a variável de desfecho."}
df %>% 
  subset(!is.na(Mobilidade)) %>%
g_t2(Mobilidade,Situação)
```

Percebe-se na tabela \@ref(tab:freqsmob) e no gráfico \@ref(fig:graphmob) uma enorme retenção em valores absolutos dos ingressantes que são do mesmo estado que a universidade.

Para o teste de associação, temos que a estatística de teste é igual a `r round((as.numeric(a$Est)),3)`, com p-valor menor que 0.01 , sendo assim, rejeitamos H0 a todos os níveis de significância, logo há evidência de que há associação entre as variáveis. Sendo uma tabela 2x2 podemos calcular a razão de prevalências que é: `r round(as.numeric(a$RPrevalencia),3)`, podemos confirmar com o no intervalo de confiança  [`r round(as.numeric(a$ICPrevalencia)[1],3)`;`r round(as.numeric(a$ICPrevalencia)[2],3)`], então podemos ver que 1 não está contido no intervalo, reforçando esta associação.

\subsection{Tabelas Marginais}
Nesta seção, iremos visualizar através de tabelas , o efeito marginal das variáveis explicativas na variável de
desfecho

```{r turnosit}
freq_table2(df,Turno,Situação,"Proporção marginal do Turno na Situação do ingressante",1)
```
Observando a  tabela \@ref(tab:turnosit), constata-se que o turno que mais evade é o noturno, com proximidade do matutino, já o que menos evade é o turno integral. 
```{r sexosit}
freq_table2(df,Sexo,Situação,"Proporção marginal do Sexo na Situação do ingressante",1)
```

Apesar de termos aproximadamente $70\%$ de mulheres ingressantes e como vimos no gráfico \@ref(fig:graphsex) existe maior retenção de estudantes do sexo feminino em valores absolutos, entretanto em termos relativos podemos notar que não há diferença significativa na permanência no curso em relação ao sexo.

```{r nacsit}
freq_table2(df,Nacionalidade,Situação,"Proporção marginal da Nacionalidade na Situação do ingressante",1)
```
O gráfico \@ref(fig:graphnac) nos informou que em valores absolutos existe mais retenção entre os ingressantes de naturalidade Brasileira e também sabemos que temos poucas pessoas que vieram do exterior ou foram naturalizadas, mas com tabela \@ref(tab:nacsit)  notamos dentre os ingressantes  os  do exterior ou nacionalizados são os que possuem maior taxa de retenção no curso de psicologia.

```{r apoiosit}
freq_table2(df,`Apoio social`,Situação,"Proporção marginal de possuir Apoio Social na Situação do ingressante",1)
```
É observado que os ingressantes que possuem apoio social, mesmo sendo poucos, são os que possuem maior retenção. 

```{r idadesit}
freq_table2(df,Idade,Situação,"Proporção marginal da Idade na Situação do ingressante",1)

```
Em relação a taxa de permanência para a variável idade, é perceptivel na faixa etária de 65 a 90 anos uma proporção maior de evasão do que nas outras faixas etárias.

```{r mobsit}
df %>% 
  subset(!is.na(Mobilidade)) %>%
freq_table2(Mobilidade,Situação,"Proporção marginal da \"Mobilidade\" na Situação do ingressante",1)
```
Para estudantes de origem do mesmo estado onde está localizada a universidade, há uma taxa de evasão relativamente menor em comparação com estudantes de origem de estados diferentes.

```{r,include=F}


 df %>% 
  ggplot(aes(Sexo)) + 
  geom_bar(aes(fill=Sexo))+
  geom_text(stat='count',aes(label=..count..),vjust=-0.5)+
  theme_minimal()+
    labs(y="Quantidade",x=df %>%
         select(Sexo) %>% 
         names(),caption = "Fonte: Elaborado pelos autores")+
  scale_fill_brewer(palette="Set2")
```


[^1]:[Acesse o Censo aqui](https://www.gov.br/inep/pt-br/areas-de-atuacao/pesquisas-estatisticas-e-indicadores/censo-da-educacao-superior)

\section{Parte do vitor}
```{r}
Associar <-function(df,v1,tipo){
  ## tipo = o -> v1 é ordinal
  ## tipo = nd -> v1 é nominal dicotômica
  ## tipo = np -> v1 é nominal politômica
  df <- drop_na(df, {{v1}})
  if(tipo == "nd"){
    Qp <- df %>%
    select({{v1}},Situação) %>%
    table() %>%
    chisq.test(correct = F)
    Asso <- df %>%
      select({{v1}},Situação) %>%
      table() %>%
      epi.2by2(method = "cross.sectional")
    lista <- list('Est' = Qp$statistic, 'p' = Qp$p.value,
                  'RPrevalencia' = Asso$massoc.summary$est[1], 
                  'RChance' = Asso$massoc.summary$est[2], 
                  'ICPrevalencia' = c(Asso$massoc.summary$lower[1],
                                      Asso$massoc.summary$upper[1]),
                  'ICChance' = c(Asso$massoc.summary$lower[2],
                                      Asso$massoc.summary$upper[2]))
    return(lista)
  }
  
  if(tipo == "np"){
    Qp <- df %>%
    select({{v1}},Situação) %>%
    table() %>%
    chisq.test(correct = F)
    lista <- list('Est' = Qp$statistic, 'p' = Qp$p.value)
    return(lista)
  }
  
  if(tipo == "o"){
    y <- df %>% 
    mutate(Evadido=fct_recode(Situação,"0"="Retido","1"="Evadido")) %>%
      select(Evadido)
    x <- df %>%
      mutate(Idade=fct_recode(Idade,"1"="[15,18]", "2" = "(18,21]", "3" = "(21,25]", 
                              "4" = "(30,65]", "5" = "(65,90]")) %>%
      select(Idade)
    rac<-cor(x,y)
    n1<- length(y)
    qcs<-(n1-1)*rac^2
    p<-1-pchisq(qcs,1)  
    lista <- list('Est' = qcs, 'p' = p)
    return(lista)
  }}


```


\section{Modelagem}
Nesta seção, iremos ajustar uma regressão logística  e floresta aleatória, temos como variável a ser predita se o estudante vai estar retido no curso, ou vai evadir, utiliza-se como covariáveis o **sexo**, **nacionalidade** e **turno** do estudante. Como observado na tabela  \@ref(tab:freqsit) tem-se que as observações nas categorias da variável resposta são muito desequilibradas, portanto a acurácia (accuracy) não será uma boa métrica para a escolha dos hiperparâmetros, assim escolhe-se a área da curva roc como métrica.


```{r}
set.seed(4)
df=df %>% 
  select(-id,-Mobilidade) %>% 
  rename(Apoio=`Apoio social`,sit="Situação")
#df %>% 
  #filter(sit=="Retido") %>%   #opcional
  #group_by(sit) %>% 
  #slice_sample(n = 50000) %>% 
  #bind_rows(sit=.,df %>% filter(sit!="Retido"))
 
```

```{r mod}
set.seed(4)
df_split=initial_split(df,0.65,strata=sit)
df_train=training(df_split)
df_test=testing(df_split)
df_vf=vfold_cv(df_train,10)
df_boot=bootstraps(df_train,5)
```



```{r prp}
df_rec1=
  df_train %>% 
  recipe(sit~.) %>% 
   step_impute_mode(Mobilidade) %>% 
    step_dummy(all_nominal_predictors())    #modelo com tudo

df_rec2=df_train %>%
  recipe(sit~Sexo+Nacionalidade+Turno) %>% 
  step_dummy(all_nominal_predictors())

df_rec3=
  df_train %>%
  recipe(sit~Sexo+Turno) %>% 
  step_dummy(all_nominal_predictors())



```

```{r models}
df_log= logistic_reg() %>% 
  set_engine("glm")

df_rf=rand_forest(trees=500,mtry = 1,min_n = 36) %>% 
  set_mode("classification") %>% 
  set_engine("ranger",importance = "impurity")




```

```{r work}
df_work=workflow_set(list(si=df_rec2,t=df_rec3),
                     list(logistic=df_log,rf=df_rf),cross=T)
```

```{r}
set.seed(4)
doParallel::registerDoParallel(cores=2)

df_tuner=df_work %>% 
  workflow_map("tune_grid",
               resamples=df_boot,
               grid=15,
               metrics=metric_set(accuracy,roc_auc,specificity),verbose=T)


rank_results(df_tuner,rank_metric = "roc_auc") %>% view()




#Validação Cruzada
ft1=df_work %>% 
  extract_workflow("t_logistic") %>% 
  fit_resamples(df_vf,metrics=metric_set(accuracy,roc_auc,sensitivity,specificity),
                control = control_resamples(save_pred = TRUE)) %>% 
  collect_predictions()

ft2=df_work %>% 
  extract_workflow("t_rf") %>% 
  fit_resamples(df_vf,metrics=metric_set(accuracy,roc_auc,sensitivity,specificity),
                control = control_resamples(save_pred = TRUE)) %>% 
  collect_predictions()
  

  

```

```{r}


library(probably)
p=function(ft,p)
{  ft %>% 
    mutate(.pred_class = make_two_class_pred(.pred_Retido, 
                                             levels(ft$sit),
                                             threshold = p),
           .pred_class = factor(.pred_class, levels = levels(ft$sit)))
}
 

  
```



```{r}


#melhores hiperparâmetros usando curva roc como métrica
#brl=df_tuner %>% 
  # extract_workflow_set_result(id='si_logistic') %>% 
  #select_best(metric='roc_auc')
#brf=df_tuner %>% 
   #extract_workflow_set_result(id='si_rf') %>% 
  #select_best(metric='roc_auc')

 


fit2=df_tuner %>% 
  extract_workflow(id='si_rf') %>% 
  fit(df_train) 

fit1=df_tuner %>% 
  extract_workflow(id='si_logistic') %>% 
  fit(df_train)
 


fitt2=fit2 %>%  
  predict(new_data = df_test, type = "prob") %>% 
  bind_cols(df_test) %>%  
<<<<<<< HEAD
  p(0.813) 
  

=======
  p(0.807) 
>>>>>>> f585e7c305b9b2b54bfe9541c783b85f65788d64
  





fitt1=fit1 %>% 
  predict(new_data = df_test, type = "prob") %>% 
  bind_cols(df_test) %>% 
  p(0.808) 



#df_tuner %>% 
 # extract_workflow(id='si_logistic') %>% 
  #finalize_workflow(brl) %>% 
  #fit(df_train) %>% 
  #p(0.808)
  #accuracy(sit,.pred_class)


```


\subsection{Regressão Logística}
<<<<<<< HEAD
A regressão logistica foi necessário tomar como $0.808$ o ponto corte para a previsão ser da categória Retido, pois como já dito existem muitas observações para essa categoria. Ajustando no conjunto de treino obtem-se como única variáve não significativa a dummie Nacionalidade  Estrangeira que tem como p-valor associado `r round(tidy(fit1)[4,5],3)`
=======
Para a regressão logística foi necessário tomar como $0.808$ o ponto de corte, assim a previsão será da categoria Retido, pois como já dito existem muitas observações para essa categoria. Ajustando no conjunto de treino obtém-se como única variável não significativa a dummie Nacionalidade Estrangeira que tem como p-valor associado `r round(tidy(fit1)[4,5],3)`, o seguinte gráfico ilustra o ajuste sobre o modelo logístico, indicando que está bem ajustado inferencialmente, dado que todos resíduos estão dentro da banda simulada.
>>>>>>> f585e7c305b9b2b54bfe9541c783b85f65788d64

```{r}
hnp::hnp(fit1$fit$fit$fit)
```




<<<<<<< HEAD
Testando no conjunto de teste que possui `r nrow(df_test)` e deixando como padrão (0.5), a acurácia nos dados de teste é `r fitt1 %>% p(0.5) %>% accuracy(sit,.pred_class) %>% .[1,3] %>% round(3)`, poderia-se cair em tentanção e dizer que é um bom modelo, mas vejamos a matriz de confusão:
=======
Testando no conjunto de teste que possui `r nrow(df_test)` observações e deixando  o ponto de corte  padrão (0.5), a acurácia nos dados de teste é `r fitt1 %>% p(0.5) %>% accuracy(sit,.pred_class) %>% .[1,3] %>% round(3)`, poderia-se cair em tentanção e dizer que é um bom modelo, mas vejamos a matriz de confusão:
>>>>>>> f585e7c305b9b2b54bfe9541c783b85f65788d64

```{r}
fitt1 %>% 
  conf_table(0.5,"Matriz de Confusão com ponte de corte 0.5")

```
Todas as previsões foram Retido!!, ou seja o modelo está completamente viesado, e possui uma alta acurácia pois a maior parte dos estudantes permaneceu no curso. Então, é necessário ajustarmos esse ponto de corte, por otimização da área da curva roc ficamos com $0.808$ , assim obtém-se área sobre a curva roc `r fitt1 %>% p(0.808) %>% roc_auc(sit,.pred_Retido) %>% .[1,3] %>% round(3)`  e matriz de confusão

```{r}
fitt1 %>% 
  conf_table(0.808,"Matriz de Confusão com ponte de corte 0.808")
```
<<<<<<< HEAD
```{r}

```



\subsection{Random Forest}
O fenomeno de predizer todos como retido, também acontece na random forest, portanto o ponto de corte foi ajustado para 
$0.81$, assim obtendo-se area sobre a curva roc igual a `r fitt2 %>% p(0.81) %>% roc_auc(sit,.pred_Retido) %>% .[1,3] %>% round(3)`  e acurácia `r fitt2 %>% p(0.81) %>% accuracy(sit,.pred_class) %>% .[1,3] %>% round(3)`. O seguinte gráfico fornece  as variáveis mais importantes no conjunto de treino

```{r}
fit2 %>% 
extract_fit_parsnip() %>% 
  vip::vip(num_features = 10)
```


 
=======



>>>>>>> f585e7c305b9b2b54bfe9541c783b85f65788d64


\subsection{Random Forest}
O fenômeno de predizer todos como retido, também acontece na random forest, portanto o ponto de corte foi ajustado para $0.81$, utilizou-se também  a quantidade de árvores de decisão (ntree)  igual a  500, e por otimização  mtry igual a 1 e min_n igual a 36. O seguinte gráfico fornece  as variáveis mais importantes no conjunto de treino.

```{r}
fit2 %>% 
extract_fit_parsnip() %>% 
  vip::vip(num_features = 10)
```


O ajuste final considerando o conjunto de teste teve  área sobre a curva roc igual a `r fitt2 %>% p(0.81) %>% roc_auc(sit,.pred_Retido) %>% .[1,3] %>% round(3)`  e acurácia `r fitt2 %>% p(0.81) %>% accuracy(sit,.pred_class) %>% .[1,3] %>% round(3)`. A matriz de confusão da random forest é dada por :

```{r}
fitt2 %>% 
  conf_table(0.81,"Matriz de Confusão com ponte de corte 0.808")
```



\section{Conclusão}
De modo geral, para os ingressantes do curso de  Psicologia do ano de 2019, podemos notar que a maioria dos ingressantes foi do sexo feminino, mas em termos da variável de desfecho não existe uma diferença significativa na taxa de evasão em relação ao sexo. A taxa de evasão para as idades em classes só não se manteve em tornos de $80\%$ para a faixa etária de 65 a 90 anos que inclusive é a maior taxa de evasão de todos os níveis das variáveis estudadas. Para a variável turno do ingressante percebe-se que grande parte foi do noturno e em relação a variável de desfecho a maior taxa de permanência (retido) foi para os estudantes do turno integral. Apesar da variável resposta ser desbalanceada  em um primeiro momento  ambos os modelos  tiveram um desempenho semelhante : mediano.