cNORM
Danilo Assis Pereira
3/14/2021
Normas usando cNORM e ctree
A ideia é fazer uma partição dos dados baseados nas diferenças entre idade e escolaridade. Para este propósito, iremos utilizar uma árvore de decisão condicional, que possui a vantagem de fazer estas partições simultaneamente e em vários nÃveis hierárquicos.
Primeiro, vamos escrever as sintaxes para obtermos a média e o DP robustos, aparados em 20%, como sugerido por Rand Wilcox.
# Funcao pra calcular a media aparada em 20%
tmean <- function(x,tr=.2,na.rm=FALSE,STAND=NULL){
if(na.rm)x<-x[!is.na(x)]
val<-mean(x,tr)
val
}
# Funcao pra calcular o desvio padrao (DP) aparado em 20%
sd_trim <- function(x,trim=0.2, const=TRUE){
# trimmed sd, where x is a matrix (column-wise)
x <- as.matrix(x)
if (const){
if (trim==0.1){const <- 0.7892}
else if (trim==0.2){const <- 0.6615}
else {warning("Voce especificou a constante correta?")}
}
else{const <- 1}
m <- apply(x,2,mean,trim)
res <- x-rep(1,nrow(x))%*%t(m)
qu <- apply(abs(res),2,quantile,1-trim)
sdtrim <- apply(matrix(res[t(abs(t(res))<=qu)]^2,ncol=ncol(x),byrow=FALSE),2,sum)
sdtrim <- sqrt(sdtrim/((nrow(x)*(1-trim)-1)))/const
return(sdtrim)
}Banco de dados
Vamos agora selecionar o banco de dados.
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.5
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(haven)
# No PC Avell
# setwd("C:/Dropbox/Laboratorio/Pedro Brandao/PDCRS")
# No Mac M1
setwd("~/Dropbox/Laboratorio/Pedro Brandao/PDCRS")
# Le o banco .sav exportado pelo SurveyMonkey
dataset = read_sav("normas_pdcrs_2021.sav")
# Limpa os nomes das variaveis
dataset <- tibble::as_tibble(dataset, .name_repair = janitor::make_clean_names)
# Exporta para o formato .csv
write_csv(dataset, file ="dados.csv")
library(readr)
dados <- read_csv("dados.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double(),
## id = col_character(),
## data_nasc = col_character(),
## pesquisador = col_character(),
## cidade = col_character(),
## estado = col_character(),
## nome = col_character(),
## idade = col_character(),
## naturalidade = col_character(),
## estado_naturalidade = col_character(),
## data_coleta = col_character(),
## exclu03_recusa_tcle = col_logical(),
## exclu05_experiencia_teste = col_logical(),
## exclu06_linguaportuguesa = col_logical(),
## exclu07_etilismo = col_logical(),
## exclu08_problema_motor_mmss = col_logical(),
## exclu10_afasia = col_logical(),
## exclu12_tce_epilepsia = col_logical()
## )
## ℹ Use `spec()` for the full column specifications.dados %>% glimpse## Rows: 597
## Columns: 192
## $ id <chr> "2019001", "2019002", "20190…
## $ data_nasc <chr> "02/05/2002", "18/02/2004", …
## $ pesquisador <chr> "Tatiane H. Werle", "Tatiane…
## $ cidade <chr> "São Carlos", "Palmitos", "P…
## $ estado <chr> "SC", "SC", "SC", "DF", "SC"…
## $ nome <chr> "Gabriela Luft Werlans", "Br…
## $ idade <chr> "16", "14", "17", "16", "18"…
## $ genero <dbl> 2, 2, 1, 2, 1, 1, 2, 1, 1, 1…
## $ naturalidade <chr> "São Carlos", "Palmitos", "P…
## $ estado_naturalidade <chr> "SC", "SC", "SC", "DF", "SC"…
## $ regiao_brasil <dbl> 5, 5, 5, 3, 5, 5, 5, 3, 3, 3…
## $ data_coleta <chr> "12/12/2018", "13/12/2018", …
## $ escolaridade <dbl> 10, 8, 11, 10, 12, 10, 12, 9…
## $ faixa_escolaridade <dbl> 5, 2, 5, 4, 5, 4, 5, 4, 4, 6…
## $ criterio_exclusao <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ exclu01_faixa_etaria <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu02_psiquiatrico <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu03_recusa_tcle <lgl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu04_deficiencia <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu05_experiencia_teste <lgl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu06_linguaportuguesa <lgl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu07_etilismo <lgl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu08_problema_motor_mmss <lgl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu09_psicose <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu10_afasia <lgl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu11_prob_aprendizagem <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ exclu12_tce_epilepsia <lgl> NA, NA, NA, NA, NA, NA, NA, …
## $ fpt_cinco_pontos_total <dbl> 18, 12, 24, 23, 22, 25, 22, …
## $ fpt_cinco_pontos_persev <dbl> 0, 1, 0, 2, 0, 1, 0, 0, 2, 0…
## $ fluencia_animais <dbl> 18, 19, 27, 12, 14, 20, 20, …
## $ fluencia_frutas <dbl> 19, 13, 19, 11, 13, 10, 12, …
## $ fluencia_alternada_pares <dbl> 9, 7, 11, 6, 6, 7, 8, 7, 10,…
## $ fluencia_alternada_indiv_total <dbl> 18, 14, 22, 12, 12, 14, 16, …
## $ ccc_brown_peterson_total <dbl> 41, 44, 41, 43, 41, 36, 45, …
## $ ccc_brown_peterson_0seg <dbl> 15, 15, 15, 15, 15, 15, 15, …
## $ ccc_brown_peterson_9seg <dbl> 13, 9, 5, 13, 11, 5, 9, 15, …
## $ ccc_brown_peterson_18seg <dbl> 9, 10, 13, 5, 10, 8, 10, 11,…
## $ ccc_brown_peterson_36seg <dbl> 4, 10, 8, 10, 5, 8, 11, 9, 5…
## $ tmt_a <dbl> 52, 31, 15, 28, 23, 25, 31, …
## $ tmt_b <dbl> 64, 59, 41, 58, 61, 57, 49, …
## $ hads_a1_tenso_contraido <dbl> 1, 2, 0, 1, 1, 1, 2, 1, 2, 1…
## $ hads_d1_gosto_pelas_coisas <dbl> 1, 0, 0, 1, 1, 1, 2, 0, 0, 0…
## $ hads_a3_medo <dbl> 2, 3, 0, 2, 0, 0, 2, 0, 2, 0…
## $ hads_d4_risada <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 0, 0…
## $ hads_a5_preocupacoes <dbl> 3, 2, 1, 1, 1, 1, 2, 2, 3, 2…
## $ hads_d6_alegre <dbl> 0, 1, 0, 1, 2, 1, 2, 0, 1, 0…
## $ hads_a7_relaxado <dbl> 2, 2, 0, 1, 2, 1, 2, 1, 2, 1…
## $ hads_d8_lento_pensar <dbl> 1, 2, 0, 0, 1, 1, 3, 1, 2, 2…
## $ hads_a9_medo <dbl> 2, 3, 1, 1, 1, 1, 3, 1, 1, 0…
## $ hads_d10_interesse_aparencia <dbl> 0, 0, 0, 1, 0, 1, 1, 0, 0, 1…
## $ hads_a11_inquieto <dbl> 1, 3, 0, 0, 1, 2, 1, 0, 0, 0…
## $ hads_d12_animado <dbl> 1, 0, 0, 0, 0, 0, 3, 2, 0, 2…
## $ hads_a13_panico <dbl> 2, 1, 0, 0, 0, 0, 2, 0, 0, 0…
## $ hads_d14_prazer <dbl> 0, 2, 0, 1, 2, 0, 3, 3, 1, 1…
## $ pdcrs_1a_tent_01_luz <dbl> 1, 0, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_1a_tent_02_seda <dbl> 0, 1, 1, 1, 1, 1, 0, 1, 1, 1…
## $ pdcrs_1a_tent_03_areia <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 1…
## $ pdcrs_1a_tent_04_cilio <dbl> 0, 1, 0, 0, 0, 1, 1, 0, 1, 0…
## $ pdcrs_1a_tent_05_arroz <dbl> 0, 0, 0, 1, 1, 0, 0, 0, 0, 0…
## $ pdcrs_1a_tent_06_gravata <dbl> 1, 1, 1, 1, 1, 0, 1, 0, 1, 0…
## $ pdcrs_1a_tent_07_quadro <dbl> 1, 0, 0, 0, 0, 0, 0, 1, 0, 1…
## $ pdcrs_1a_tent_08_bicicleta <dbl> 0, 1, 0, 1, 0, 0, 0, 1, 1, 1…
## $ pdcrs_1a_tent_09_estrela <dbl> 1, 0, 0, 0, 0, 0, 0, 1, 1, 0…
## $ pdcrs_1a_tent_10_leao <dbl> 0, 0, 1, 0, 1, 1, 0, 1, 1, 1…
## $ pdcrs_1a_tent_11_anel <dbl> 1, 1, 1, 0, 0, 0, 0, 1, 1, 0…
## $ pdcrs_1a_tent_12_perfume <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 1, 0…
## $ pdcrs_2a_tent_01_luz <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_2a_tent_02_seda <dbl> 0, 0, 1, 1, 0, 1, 1, 1, 1, 1…
## $ pdcrs_2a_tent_03_areia <dbl> 1, 1, 0, 1, 1, 0, 1, 1, 0, 0…
## $ pdcrs_2a_tent_04_cilio <dbl> 1, 1, 1, 1, 0, 0, 0, 1, 1, 1…
## $ pdcrs_2a_tent_05_arroz <dbl> 0, 1, 1, 1, 0, 0, 1, 1, 0, 1…
## $ pdcrs_2a_tent_06_gravata <dbl> 1, 1, 1, 1, 0, 0, 1, 0, 1, 0…
## $ pdcrs_2a_tent_07_quadro <dbl> 1, 0, 1, 0, 1, 0, 0, 0, 0, 1…
## $ pdcrs_2a_tent_08_bicicleta <dbl> 1, 1, 0, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_2a_tent_09_estrela <dbl> 1, 0, 0, 1, 0, 1, 0, 1, 1, 0…
## $ pdcrs_2a_tent_10_leao <dbl> 0, 1, 1, 1, 0, 1, 1, 1, 1, 0…
## $ pdcrs_2a_tent_11_anel <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_2a_tent_12_perfume <dbl> 1, 1, 1, 0, 1, 0, 0, 0, 1, 1…
## $ pdcrs_3a_tent_01_luz <dbl> 1, 1, 1, 1, 1, 0, 1, 1, 1, 1…
## $ pdcrs_3a_tent_02_seda <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 1, 1…
## $ pdcrs_3a_tent_03_areia <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 0, 1…
## $ pdcrs_3a_tent_04_cilio <dbl> 1, 1, 1, 1, 0, 1, 0, 1, 1, 0…
## $ pdcrs_3a_tent_05_arroz <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 0, 0…
## $ pdcrs_3a_tent_06_gravata <dbl> 1, 1, 1, 1, 0, 1, 0, 0, 1, 1…
## $ pdcrs_3a_tent_07_quadro <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_3a_tent_08_bicicleta <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 1, 1…
## $ pdcrs_3a_tent_09_estrela <dbl> 1, 1, 0, 1, 1, 0, 0, 1, 1, 1…
## $ pdcrs_3a_tent_10_leao <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 0, 1…
## $ pdcrs_3a_tent_11_anel <dbl> 1, 1, 1, 1, 1, 0, 1, 0, 1, 0…
## $ pdcrs_3a_tent_12_perfume <dbl> 0, 1, 1, 1, 0, 0, 1, 1, 1, 1…
## $ pdcrs_s1_evoc_imediata <dbl> 12, 13, 12, 13, 10, 9, 9, 11…
## $ pdcrs_nomeacao_01_babador <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_02_vela <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_03_cereja <dbl> 1, 1, 1, 1, 1, 1, 1, 0, 1, 1…
## $ pdcrs_nomeacao_04_banqueta <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 0…
## $ pdcrs_nomeacao_05_ancora <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_06_tartaruga <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_07_pipa <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_08_aquario <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_09_lampada <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_10_violao <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_11_fivela <dbl> 1, 0, 1, 1, 0, 0, 1, 1, 1, 1…
## $ pdcrs_nomeacao_12_crina <dbl> 1, 1, 1, 0, 1, 1, 1, 0, 0, 0…
## $ pdcrs_nomeacao_13_anzol <dbl> 1, 1, 1, 1, 1, 1, 0, 1, 1, 1…
## $ pdcrs_nomeacao_14_chave_fenda <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_15_biombo <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1…
## $ pdcrs_nomeacao_16_alfinete <dbl> 1, 0, 1, 1, 0, 1, 1, 0, 0, 1…
## $ pdcrs_nomeacao_17_chocalho <dbl> 1, 1, 1, 0, 0, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_18_casco <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_19_extintor <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_nomeacao_20_trinco <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item01 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item02 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item03 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item04 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item05 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item06 <dbl> 1, 1, 1, 1, 1, 1, 1, 0, 1, 1…
## $ pdcrs_atencao_item07 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item08 <dbl> 1, 1, 1, 1, 1, 1, 1, 0, 0, 1…
## $ pdcrs_atencao_item09 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_atencao_item10 <dbl> 1, 0, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_mem_operac_item01 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_mem_operac_item02 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_mem_operac_item03 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_mem_operac_item04 <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1…
## $ pdcrs_mem_operac_item05 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_mem_operac_item06 <dbl> 1, 1, 1, 1, 1, 1, 1, 0, 1, 0…
## $ pdcrs_mem_operac_item07 <dbl> 1, 1, 1, 1, 1, 0, 1, 1, 0, 0…
## $ pdcrs_mem_operac_item08 <dbl> 1, 1, 1, 1, 1, 0, 0, 0, 1, 0…
## $ pdcrs_mem_operac_item09 <dbl> 1, 0, 0, 0, 0, 0, 1, 1, 1, 1…
## $ pdcrs_mem_operac_item10 <dbl> 0, 1, 0, 0, 0, 0, 0, 1, 0, 0…
## $ pdcrs_relogio_item01_parece_relog <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_relogio_item02_sem_setores <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_relogio_item03_simetria <dbl> 1, 1, 1, 1, 1, 1, 1, 0, 1, 0…
## $ pdcrs_relogio_item04_num_1a12 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_relogio_item05_sequencia_ok <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_relogio_item06_dois_pont <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 0…
## $ pdcrs_relogio_item07_setas <dbl> 1, 1, 0, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_relogio_item08_pont_horas_menor <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_relogio_item09_sem_palavra <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_relogio_item10_sem_num25 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item01_parece_relog <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item02_sem_setores <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item03_simetria <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 0…
## $ pdcrs_copia_relog_item04_num_1a12 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item05_sequencia_ok <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item06_dois_pont <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item07_setas <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item08_pont_horas_menor <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item09_sem_palavra <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_copia_relog_item10_sem_num25 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ pdcrs_evoc_tardia_01_luz <dbl> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1…
## $ pdcrs_evoc_tardia_02_seda <dbl> 1, 1, 1, 1, 0, 0, 1, 1, 1, 1…
## $ pdcrs_evoc_tardia_03_areia <dbl> 1, 1, 0, 1, 1, 0, 1, 1, 0, 1…
## $ pdcrs_evoc_tardia_04_cilio <dbl> 0, 1, 1, 1, 1, 1, 0, 1, 1, 0…
## $ pdcrs_evoc_tardia_05_arroz <dbl> 1, 1, 1, 1, 1, 1, 0, 0, 0, 0…
## $ pdcrs_evoc_tardia_06_gravata <dbl> 0, 1, 1, 1, 0, 1, 0, 0, 1, 1…
## $ pdcrs_evoc_tardia_07_quadro <dbl> 1, 1, 1, 0, 0, 1, 1, 1, 1, 1…
## $ pdcrs_evoc_tardia_08_bicicleta <dbl> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1…
## $ pdcrs_evoc_tardia_09_estrela <dbl> 1, 0, 1, 1, 0, 0, 1, 0, 0, 1…
## $ pdcrs_evoc_tardia_10_leao <dbl> 0, 0, 1, 0, 1, 1, 0, 1, 1, 0…
## $ pdcrs_evoc_tardia_11_anel <dbl> 1, 1, 0, 1, 1, 0, 1, 1, 0, 0…
## $ pdcrs_evoc_tardia_12_perfume <dbl> 1, 1, 0, 0, 0, 0, 0, 1, 0, 0…
## $ pdcrs_s8_fluencia_alternada <dbl> 13, 13, 13, 11, 13, 7, 13, 1…
## $ pdcrs_s9_fluencia_acoes <dbl> 16, 17, 24, 13, 13, 9, 15, 1…
## $ pdcfrs01lidardinheiro <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs02adm_contas <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs03plan_ferias <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs04recibos_faturas <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs05medicamentos <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs06ativ_cotidianas <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs07eletrodomesticos <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs08transp_publico <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs09resolv_imprevistos <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs10explicaroquediz <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs11entenderoquele <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcfrs12tel_celular <dbl> NA, NA, NA, NA, NA, NA, NA, …
## $ pdcrs_memoria1 <dbl> 6, 5, 5, 5, 5, 5, 4, 7, 9, 6…
## $ pdcrs_memoria2 <dbl> 8, 9, 9, 10, 6, 6, 8, 9, 9, …
## $ pdcrs_memoria3 <dbl> 11, 12, 11, 12, 9, 8, 6, 10,…
## $ pdcrs_s1_memoria_imediata <dbl> 11, 12, 11, 12, 9, 8, 8, 10,…
## $ pdcrs_s2_nomeacao <dbl> 20, 17, 19, 16, 16, 18, 18, …
## $ pdcrs_s3_atencao_sust <dbl> 10, 9, 10, 10, 10, 10, 10, 8…
## $ pdcrs_s4_memoria_operac <dbl> 9, 9, 8, 7, 8, 6, 8, 8, 8, 6…
## $ pdcrs_s5_desenho_relogio <dbl> 10, 10, 9, 10, 10, 10, 10, 9…
## $ pdcrs_s6_copia_relogio <dbl> 10, 10, 10, 10, 10, 10, 10, …
## $ pdcrs_s7_memoria_tardia <dbl> 7, 8, 9, 8, 4, 7, 6, 7, 7, 7…
## $ pdcrs_s8_fluencia_alternada_2 <dbl> 13, 13, 13, 11, 13, 7, 13, 1…
## $ pdcrs_s9_fluencia_acoes_2 <dbl> 16, 17, 24, 13, 13, 9, 15, 1…
## $ pdcrs_frontal_subcortical_score <dbl> 76, 78, 84, 71, 67, 57, 70, …
## $ pdcrs_posterior_cortical_score <dbl> 30, 27, 29, 26, 26, 28, 28, …
## $ pdcrs_total_score <dbl> 106, 105, 113, 97, 93, 85, 9…# Transforma a idade em variavel numerica
dados$idade <- as.numeric(dados$idade)## Warning: NAs introduzidos por coerção# Plota a distribuicao da idade
hist(dados$idade, breaks = 50)
# Plota a distribuicao da escolaridade
hist(dados$escolaridade, breaks = 50)
# Seleciona apenas as variaveis de interesse
df <- dados %>% select(pdcrs_total_score, idade, escolaridade) %>%
na.omit()
# Renomeia a variavel predita pra padronizar a sintaxe
names(df) <- c("score", "idade", "escolaridade")
# Transforma as variaveis integrais em numericas
df <- df %>% mutate(across(where(is.integer), as.numeric))
df <- as.data.frame(df)Arvore de decisao condicional
A árvore de decisão condicional pode gerar uma partição com overfit, mudando quando os dados são atualizados. Por isso, faremos uma cross-validation. Os resultados apresentados em cada nodo (mu e sd) poderão ser utilizados no site https://sapsi.shinyapps.io/bayesNeuro/ para se estimar o escore z e percentil, bem como calcular a probabidade do caso ser estatisticamente diferente do grupo controle.
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftlibrary(party)## Loading required package: grid## Loading required package: mvtnorm## Loading required package: modeltools## Loading required package: stats4## Loading required package: strucchange## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich##
## Attaching package: 'strucchange'## The following object is masked from 'package:stringr':
##
## boundary# Arvores condicionais com k-fold cross-validation (sem overfit)
set.seed(666)
model <- train(
score ~ idade + escolaridade,
data = df,
method = 'ctree2',
trControl = trainControl("cv", number = 5, classProbs = FALSE),
tuneGrid = expand.grid(maxdepth = 3, mincriterion = 0.95)
)
plot(model$finalModel)
# Identifica o nodo final para cada observacao
nodes <- where(model$finalModel)
# Calcula as estatisticas e modelos necessarios para a funcao da tabela normativa
sum_stats <- rbind(mu=tapply(df$score, nodes, tmean),
sd=tapply(df$score, nodes, sd_trim))
sum_stats## 4 5 7 8 11 12 14 15
## mu 69.83333 78.95455 90.79070 96.17073 66.69565 54.77778 75.87654 90.84848
## sd 13.24982 11.07585 12.65868 11.52433 9.84906 12.99519 16.84637 15.41901Funcao normativa
# Calcula a distribuicao empirica (ECDF) da variavel
ecdfs <- lapply(split(df, nodes), function(x)
ecdf(x$score))
# Regressao linear para cada nodo
reg_models <- lapply(split(df, nodes), function(x)
lm(score ~ idade + escolaridade, data=x))
# Define a funcao da tabela normativa
# O argumento 'case' precisa ser um data.frame com 3 colunas: 'score', 'escolaridade' e 'idade'
norm_table <- function(case, method='stats') {
if (!any(class(case) == 'data.frame') | (nrow(case) == 0) | (ncol(case) != 3) |
!all(names(case) %in% c('score', 'idade', 'escolaridade'))) {
stop('"case" precisa ser um objeto data.frame com pelo menos uma observacao e tres colunas, "score", "escolaridade" e "idade"')
}
node <- where(model$finalModel, case)
if (method == 'stats'){
out <- lapply(1:nrow(case), function(x) {
data.frame(
Escore=case$score[x],
Z=(case$score[x] - sum_stats[1, as.character(node[x])])/sum_stats[2, as.character(node[x])],
Percentil.Normal=pnorm((case$score[x] - sum_stats[1, as.character(node[x])])/sum_stats[2, as.character(node[x])]),
Percentil.EmpÃrico=ecdfs[[as.character(node[x])]](case$score[x])
)
})
do.call(rbind, out)
} else if (method == 'lm') {
out <- lapply(1:nrow(case), function(x) {
data.frame(
Escore=case$score[x],
Z=(case$score[x] - predict(reg_models[[as.character(node[x])]], case[x, ]))/
sigma(reg_models[[as.character(node[x])]]),
Percentil.Normal=pnorm((case$score[x] - predict(reg_models[[as.character(node[x])]], case[x, ]))/
sigma(reg_models[[as.character(node[x])]])),
Percentil.EmpÃrico=ecdf(residuals(reg_models[[as.character(node[x])]]))(case$score[x] - predict(reg_models[[as.character(node[x])]], case[x, ]))
)
})
do.call(rbind, out)
} else {
stop('Metodo deve ser "stats" ou "lm"')
}
}
# Testa a funcao
# A funcao funciona com observacoes multiplas simultaneamente!
test.df <- df[1:10,]
# Mostra os 10 primeiros casos
norm_table(test.df, method='stats')## Escore Z Percentil.Normal Percentil.EmpÃrico
## 1 106 1.2014919 0.8852198 0.9225352
## 2 105 1.1224947 0.8691739 0.8943662
## 3 113 1.4603248 0.9278996 0.9558824
## 4 97 0.4905173 0.6881161 0.7183099
## 5 93 -0.2751337 0.3916068 0.4117647
## 6 85 -0.4574487 0.3236743 0.3380282
## 7 98 0.1587310 0.5630596 0.5882353
## 8 94 0.2535258 0.6000690 0.6197183
## 9 103 0.5925956 0.7232741 0.7794118
## 10 88 -0.7089983 0.2391628 0.2352941norm_table(test.df, method='lm')## Escore Z Percentil.Normal Percentil.EmpÃrico
## 1 106 1.1488896 0.8746993 0.9225352
## 2 105 0.7960429 0.7869964 0.8169014
## 3 113 1.4545241 0.9270995 0.9338235
## 4 97 0.5117036 0.6955708 0.7042254
## 5 93 -0.2796098 0.3898885 0.3823529
## 6 85 -0.3378778 0.3677276 0.3098592
## 7 98 0.1481304 0.5588801 0.5441176
## 8 94 0.2246139 0.5888602 0.5704225
## 9 103 0.6090725 0.7287618 0.7573529
## 10 88 -0.7073499 0.2396745 0.2205882Calcula as normas usando cNORM
Dentro de cada nodo da árvore de decisão, obtida a partir dos dados score, idade e escolaridade, iremos obter a técnica de expansão de Taylor. Para mais informações consultar http://www.psychometrica.de/cNorm_en.html
library(cNORM)## Good morning star-shine, cNORM says 'Hello!'##
## Attaching package: 'cNORM'## The following object is masked from 'package:ggplot2':
##
## derivenode_mean <- as.vector(predict(model$finalModel, df))
node_val <- as.vector(where(model$finalModel, df))
names(node_mean) <- node_val
set.seed(666)
node_mean_cat <- as.numeric(factor(node_mean))
# Curvas percentilicas observadas e preditas
m_norm <- cnorm(raw = df$score, group = node_mean)## Warning in rankByGroup(raw = raw, group = group, scale = scale, weights =
## weights, : The dataset includes cases, whose percentile depends on less than
## 30 cases (minimum is 8). Please check the distribution of the cases over the
## grouping variable. The confidence of the norm scores is low in that part of the
## scale. Consider redividing the cases over the grouping variable. In cases of
## disorganized percentile curves after modeling, it might help to reduce the 'k'
## parameter.## Multiple R2 between raw score and explanatory variable: R2 = 0.4243
## Final solution: 5 terms
## R-Square Adj. = 0.974963
## Final regression model: raw ~ A3 + A4 + L1A2 + L2A4 + L4A4
## Regression function: raw ~ 15.28102932 + (-0.000234649196*A3) + (1.949358597e-06*A4) + (0.0004683356908*L1A2) + (-4.827287305e-10*L2A4) + (2.054914004e-14*L4A4)
## Raw Score RMSE = 2.85001##
## Use 'printSubset(model)' to get detailed information on the different solutions, 'plotPercentiles(model) to display percentile plot, plotSubset(model)' to inspect model fit.
# Cross-validation pra obter o N de termos que sera usado na linha abaixo
# cnorm.cv(m_norm$data)
# Calcula as normas cNORM usando o escore padrão WESCHLER (média 10 e DP 3)
m_norm <- cnorm(raw = df$score, group = node_mean, terms=4, scale=c(10, 3))## Warning in rankByGroup(raw = raw, group = group, scale = scale, weights =
## weights, : The dataset includes cases, whose percentile depends on less than
## 30 cases (minimum is 8). Please check the distribution of the cases over the
## grouping variable. The confidence of the norm scores is low in that part of the
## scale. Consider redividing the cases over the grouping variable. In cases of
## disorganized percentile curves after modeling, it might help to reduce the 'k'
## parameter.## Multiple R2 between raw score and explanatory variable: R2 = 0.4243
## User specified solution: 4 terms
## R-Square Adj. = 0.973838
## Final regression model: raw ~ L1 + L1A3 + L2A4 + L3A3
## Regression function: raw ~ 18.11080453 + (1.857949827*L1) + (1.439479168e-05*L1A3) + (-1.044256332e-08*L2A4) + (2.51423236e-08*L3A3)
## Raw Score RMSE = 2.91582##
## Use 'printSubset(model)' to get detailed information on the different solutions, 'plotPercentiles(model) to display percentile plot, plotSubset(model)' to inspect model fit.
# Mostra o R2 para cada numero de preditor
plotSubset(m_norm)
# Plota os escores observados vs. ajustados
plotRaw(m_norm)
# Plota os escores observados vs. ajustados pela norma
plotNorm(m_norm)
# Testa para o sujeito 1
predictNorm(raw = df[1,1][[1]], A=predict(model$finalModel, df[1, ])[[1]], model=m_norm)## [1] 13.75516# Testa para o sujeito 142
predictNorm(raw = df[142,1][[1]], A=predict(model$finalModel, df[142, ])[[1]], model=m_norm)## [1] 6.585025# Testa a funcao
normTable(predict(model$finalModel, df[142, ])[[1]], model = m_norm)## norm raw percentile
## 1 2.5 28.59539 0.6209665
## 2 2.8 29.78533 0.8197536
## 3 3.1 30.96316 1.0724110
## 4 3.4 32.12960 1.3903448
## 5 3.7 33.28538 1.7864421
## 6 4.0 34.43124 2.2750132
## 7 4.3 35.56789 2.8716560
## 8 4.6 36.69607 3.5930319
## 9 4.9 37.81650 4.4565463
## 10 5.2 38.92990 5.4799292
## 11 5.5 40.03702 6.6807201
## 12 5.8 41.13857 8.0756659
## 13 6.1 42.23528 9.6800485
## 14 6.4 43.32787 11.5069670
## 15 6.7 44.41709 13.5666061
## 16 7.0 45.50364 15.8655254
## 17 7.3 46.58827 18.4060125
## 18 7.6 47.67169 21.1855399
## 19 7.9 48.75464 24.1963652
## 20 8.2 49.83784 27.4253118
## 21 8.5 50.92202 30.8537539
## 22 8.8 52.00791 34.4578258
## 23 9.1 53.09623 38.2088578
## 24 9.4 54.18771 42.0740291
## 25 9.7 55.28308 46.0172163
## 26 10.0 56.38307 50.0000000
## 27 10.3 57.48839 53.9827837
## 28 10.6 58.59979 57.9259709
## 29 10.9 59.71799 61.7911422
## 30 11.2 60.84371 65.5421742
## 31 11.5 61.97768 69.1462461
## 32 11.8 63.12063 72.5746882
## 33 12.1 64.27329 75.8036348
## 34 12.4 65.43639 78.8144601
## 35 12.7 66.61064 81.5939875
## 36 13.0 67.79678 84.1344746
## 37 13.3 68.99554 86.4333939
## 38 13.6 70.20764 88.4930330
## 39 13.9 71.43382 90.3199515
## 40 14.2 72.67479 91.9243341
## 41 14.5 73.93128 93.3192799
## 42 14.8 75.20403 94.5200708
## 43 15.1 76.49375 95.5434537
## 44 15.4 77.80119 96.4069681
## 45 15.7 79.12705 97.1283440
## 46 16.0 80.47208 97.7249868
## 47 16.3 81.83700 98.2135579
## 48 16.6 83.22253 98.6096552
## 49 16.9 84.62941 98.9275890
## 50 17.2 86.05836 99.1802464
## 51 17.5 87.51010 99.3790335# Monta a funcao de predicao para cada nodo da arvore
tt <- lapply(sort(unique(node_mean)), function(x) normTable(x, model = m_norm, step = 1, minNorm = 1, maxNorm=19))
# Monta as tabelas normativas, uma para cada nodo da arvore
normTabs <- lapply(tt, function(x) data.frame(Norma=x$norm,
Escore.Inf=c(round(x$raw)),
Escore.Sup=c(round(x$raw[-1])-1, NA),
Percentil=x$percentile))Normas para cada nodo da árvore
As normas podem ser interpretadas usando a calculadora de conversão psicométrica http://www.psychometrica.de/normwertrechner_en.html
Faixas dos escores padrões: 1 a 4 = Muito baixo; 5 a 7 = Baixo; 8 a 13 = Médio; 14 a 16 = Alto; 17 a 19 = Muito alto.
Aqui é importante observar a inequalidade de Jensen: a média de uma variável aleatória transformada é sempre maior ou igual à transformação da média da variável aleatória. Como as distribuições não são simétricas, a média do resultado da aplicação do cNORM será, quase sempre, diferente da média empÃrica, mesmo a robusta. O teorema é uma lei da probabilidade que diz o seguinte: f(E(x)) <= E(f(x)). Como o cNORM faz uma transformação normal nos escores brutos para ajustar o modelo, o valor esperado das normas vai ser sempre maior que o valor esperado nos escores brutos.
# Mostra as tabelas normativas para cada nodo da árvore de decisão condicional, usando a escala padrão "Weschler" com média 10 e DP 3.
plot(model$finalModel)
normTabs## [[1]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 22 26 0.1349898
## 2 2 27 30 0.3830381
## 3 3 31 33 0.9815329
## 4 4 34 37 2.2750132
## 5 5 38 41 4.7790352
## 6 6 42 45 9.1211220
## 7 7 46 48 15.8655254
## 8 8 49 52 25.2492538
## 9 9 53 55 36.9441340
## 10 10 56 59 50.0000000
## 11 11 60 63 63.0558660
## 12 12 64 67 74.7507462
## 13 13 68 71 84.1344746
## 14 14 72 75 90.8788780
## 15 15 76 79 95.2209648
## 16 16 80 84 97.7249868
## 17 17 85 89 99.0184671
## 18 18 90 94 99.6169619
## 19 19 95 NA 99.8650102
##
## [[2]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 24 28 0.1349898
## 2 2 29 33 0.3830381
## 3 3 34 38 0.9815329
## 4 4 39 43 2.2750132
## 5 5 44 47 4.7790352
## 6 6 48 52 9.1211220
## 7 7 53 56 15.8655254
## 8 8 57 60 25.2492538
## 9 9 61 64 36.9441340
## 10 10 65 69 50.0000000
## 11 11 70 73 63.0558660
## 12 12 74 77 74.7507462
## 13 13 78 82 84.1344746
## 14 14 83 87 90.8788780
## 15 15 88 92 95.2209648
## 16 16 93 97 97.7249868
## 17 17 98 103 99.0184671
## 18 18 104 109 99.6169619
## 19 19 110 NA 99.8650102
##
## [[3]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 24 29 0.1349898
## 2 2 30 35 0.3830381
## 3 3 36 40 0.9815329
## 4 4 41 45 2.2750132
## 5 5 46 50 4.7790352
## 6 6 51 54 9.1211220
## 7 7 55 59 15.8655254
## 8 8 60 63 25.2492538
## 9 9 64 67 36.9441340
## 10 10 68 72 50.0000000
## 11 11 73 76 63.0558660
## 12 12 77 80 74.7507462
## 13 13 81 85 84.1344746
## 14 14 86 90 90.8788780
## 15 15 91 95 95.2209648
## 16 16 96 101 97.7249868
## 17 17 102 107 99.0184671
## 18 18 108 113 99.6169619
## 19 19 114 NA 99.8650102
##
## [[4]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 26 32 0.1349898
## 2 2 33 39 0.3830381
## 3 3 40 45 0.9815329
## 4 4 46 51 2.2750132
## 5 5 52 56 4.7790352
## 6 6 57 61 9.1211220
## 7 7 62 66 15.8655254
## 8 8 67 71 25.2492538
## 9 9 72 75 36.9441340
## 10 10 76 80 50.0000000
## 11 11 81 84 63.0558660
## 12 12 85 89 74.7507462
## 13 13 90 94 84.1344746
## 14 14 95 99 90.8788780
## 15 15 100 104 95.2209648
## 16 16 105 110 97.7249868
## 17 17 111 116 99.0184671
## 18 18 117 123 99.6169619
## 19 19 124 NA 99.8650102
##
## [[5]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 27 33 0.1349898
## 2 2 34 41 0.3830381
## 3 3 42 47 0.9815329
## 4 4 48 53 2.2750132
## 5 5 54 59 4.7790352
## 6 6 60 64 9.1211220
## 7 7 65 69 15.8655254
## 8 8 70 74 25.2492538
## 9 9 75 78 36.9441340
## 10 10 79 83 50.0000000
## 11 11 84 87 63.0558660
## 12 12 88 92 74.7507462
## 13 13 93 97 84.1344746
## 14 14 98 102 90.8788780
## 15 15 103 107 95.2209648
## 16 16 108 113 97.7249868
## 17 17 114 119 99.0184671
## 18 18 120 125 99.6169619
## 19 19 126 NA 99.8650102
##
## [[6]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 29 38 0.1349898
## 2 2 39 48 0.3830381
## 3 3 49 56 0.9815329
## 4 4 57 63 2.2750132
## 5 5 64 69 4.7790352
## 6 6 70 75 9.1211220
## 7 7 76 80 15.8655254
## 8 8 81 85 25.2492538
## 9 9 86 89 36.9441340
## 10 10 90 93 50.0000000
## 11 11 94 97 63.0558660
## 12 12 98 101 74.7507462
## 13 13 102 105 84.1344746
## 14 14 106 109 90.8788780
## 15 15 110 114 95.2209648
## 16 16 115 119 97.7249868
## 17 17 120 124 99.0184671
## 18 18 125 125 99.6169619
## 19 19 126 NA 99.8650102
##
## [[7]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 30 39 0.1349898
## 2 2 40 48 0.3830381
## 3 3 49 57 0.9815329
## 4 4 58 64 2.2750132
## 5 5 65 70 4.7790352
## 6 6 71 76 9.1211220
## 7 7 77 81 15.8655254
## 8 8 82 86 25.2492538
## 9 9 87 90 36.9441340
## 10 10 91 94 50.0000000
## 11 11 95 98 63.0558660
## 12 12 99 102 74.7507462
## 13 13 103 106 84.1344746
## 14 14 107 110 90.8788780
## 15 15 111 114 95.2209648
## 16 16 115 119 97.7249868
## 17 17 120 124 99.0184671
## 18 18 125 125 99.6169619
## 19 19 126 NA 99.8650102
##
## [[8]]
## Norma Escore.Inf Escore.Sup Percentil
## 1 1 32 43 0.1349898
## 2 2 44 53 0.3830381
## 3 3 54 62 0.9815329
## 4 4 63 70 2.2750132
## 5 5 71 77 4.7790352
## 6 6 78 83 9.1211220
## 7 7 84 88 15.8655254
## 8 8 89 93 25.2492538
## 9 9 94 96 36.9441340
## 10 10 97 100 50.0000000
## 11 11 101 103 63.0558660
## 12 12 104 106 74.7507462
## 13 13 107 109 84.1344746
## 14 14 110 112 90.8788780
## 15 15 113 115 95.2209648
## 16 16 116 118 97.7249868
## 17 17 119 122 99.0184671
## 18 18 123 125 99.6169619
## 19 19 126 NA 99.8650102ggparty
library(ggparty)## Loading required package: partykit## Loading required package: libcoin##
## Attaching package: 'partykit'## The following objects are masked from 'package:party':
##
## cforest, ctree, ctree_control, edge_simple, mob, mob_control,
## node_barplot, node_bivplot, node_boxplot, node_inner, node_surv,
## node_terminal, varimptr_tree <- lmtree(score ~ idade | escolaridade, data = df,
caseweights = FALSE)
plot(tr_tree)
ggparty(tr_tree) +
geom_edge() +
geom_edge_label() +
geom_node_splitvar() +
geom_node_plot(gglist = list(geom_point(aes(x = idade,
y = score),
alpha = 0.8),
theme_bw(base_size = 10)),
shared_axis_labels = TRUE,
legend_separator = TRUE,
# predict based on variable
predict = "idade",
# graphical parameters for geom_line of predictions
predict_gpar = list(col = "blue",
size = 1.2)
)
tr_tree <- lmtree(score ~ escolaridade | idade, data = df,
caseweights = FALSE)
ggparty(tr_tree) +
geom_edge() +
geom_edge_label() +
geom_node_splitvar() +
geom_node_plot(gglist = list(geom_point(aes(x = escolaridade,
y = score),
alpha = 0.8),
theme_bw(base_size = 10)),
shared_axis_labels = TRUE,
legend_separator = TRUE,
# predict based on variable
predict = "escolaridade",
# graphical parameters for geom_line of predictions
predict_gpar = list(col = "blue",
size = 1.2)
)