Prática R: Ler e Descrever Arquivo de Dados

MSP4112: Estatística Aplicada à Medicina

Bastão de Asclépio & Distribuição Normal

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))

options(warn=-1)
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(data.table, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(devtools, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(finalfit, warn.conflicts=FALSE))
suppressMessages(library(FSA, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(ggpubr, warn.conflicts=FALSE))
suppressMessages(library(gplots, warn.conflicts=FALSE))
suppressMessages(library(Hmisc, warn.conflicts=FALSE))
suppressMessages(library(htmltools, warn.conflicts=FALSE))
suppressMessages(library(kableExtra, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lattice, warn.conflicts=FALSE))
suppressMessages(library(lubridate, warn.conflicts=FALSE))
suppressMessages(library(microdatasus, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(read.dbc, warn.conflicts=FALSE)) # devtools::install_github("danicat/read.dbc")
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(remotes, warn.conflicts=FALSE))
suppressMessages(library(reshape2, warn.conflicts=FALSE))
suppressMessages(library(sjPlot, warn.conflicts=FALSE))
suppressMessages(library(summarytools, warn.conflicts=FALSE))

Material

• HTML de R Markdown em RPubs

Pensamento

Vídeo: Clima “atípico” em Angola deixa menina do tempo preocupada: temperaturas mínimas mais altas que as máximas!

https://www.youtube.com/watch?v=Jp0jJMPWl6g&t=2s

R

A linguagem R é uma linguagem de programação e ambiente de desenvolvimento voltado principalmente para análise estatística, manipulação de dados, visualização e criação de gráficos. Ela foi desenvolvida inicialmente por Ross Ihaka e Robert Gentleman na Universidade de Auckland, Nova Zelândia, em 2000.

O R é uma linguagem interpretada e de código aberto, o que significa que seu código-fonte está disponível livremente para que os usuários possam visualizar, modificar e distribuir. Essa característica da linguagem R encoraja uma comunidade colaborativa, levando a um grande ecossistema de pacotes e extensões criados por diversos desenvolvedores em todo o mundo.

Principais características da linguagem R:

1. Análise Estatística: R é amplamente utilizado em análise estatística e fornece uma grande variedade de funções e bibliotecas estatísticas para realizar cálculos, testes de hipóteses, regressão, análise de séries temporais e muitos outros procedimentos estatísticos.

2. Manipulação de Dados: O R oferece recursos avançados para importar, exportar e manipular dados, incluindo a capacidade de trabalhar com estruturas de dados como vetores, matrizes, data frames e listas.

3. Visualização de Dados: A linguagem R possui bibliotecas gráficas poderosas para criação de gráficos e visualizações estáticas e interativas, permitindo representar os dados de forma eficiente e atraente.

4. Programação Funcional: R é uma linguagem funcional, o que significa que ela suporta funções de primeira classe e operações de alto nível em funções.

5. Extensibilidade: Como mencionado anteriormente, R é altamente extensível por meio de pacotes. Há milhares de pacotes disponíveis para diferentes finalidades, permitindo que os usuários ampliem as funcionalidades da linguagem de acordo com suas necessidades.

6. Comunidade Ativa: A comunidade de usuários e desenvolvedores do R é ativa e vibrante. Fóruns online, grupos de discussão e conferências são realizadas regularmente para compartilhar conhecimentos, problemas e soluções.

A sintaxe do R é baseada em comandos em estilo inglês e é relativamente fácil de aprender, mesmo para pessoas que não têm uma formação em programação. Sua versatilidade e a ampla gama de pacotes disponíveis tornam o R uma das linguagens mais populares para análise de dados e estatística em diferentes campos, como ciência de dados, pesquisas acadêmicas, finanças, bioinformática, entre outros.

Até setembro de 2021, o Comprehensive R Archive Network (CRAN), o repositório oficial de pacotes R, continha mais de 17.000 pacotes disponíveis para download e uso. No entanto, é importante observar que a quantidade de pacotes R pode mudar ao longo do tempo, pois novos pacotes são desenvolvidos e adicionados regularmente, enquanto outros podem ser descontinuados ou removidos do repositório.

Além do CRAN, também existem outros repositórios de pacotes R, como o Bioconductor (especializado em bioinformática) e o GitHub, onde muitos desenvolvedores disponibilizam seus pacotes. Portanto, a quantidade total de pacotes R disponíveis pode ser ainda maior do que o número presente no CRAN.

Sumário

• Leitura e gravação de arquivo de dados
• Tipo de variável
  
• Estatística descritiva gráfica e numérica
• Estrutura do arquivo de dados com medidas repetidas: wide e long

Quatro problemas

1. Ler e descrever Biometria_FMUSP.xlsx;
2. Ler e transformar a estrutura de dados de NewDrug.xlsx de formato wide para long e vice-versa;
3. Ler dados em colunas e tabelas em script R;
4. Ler e descrever arquivo de dados grande DOBR2020.rds (documentação em Estrutura_SIM_para_CD.pdf).

Ler arquivo de dados em planilha Excel

Os dados da planilha Biometria_FMUSP.xlsx foram coletados pelos docentes dos estudantes do segundo ano de uma mesma disciplina do curso de Medicina da FMUSP em três anos consecutivos.

As variáveis do arquivo são:

• ID: idenficador do(a) estudante
  
• Ano da coleta dos dados: 1, 2, 3
  
• Turma: A, B
  
• Sexo: Feminino, Masculino
• Mao: Destro, Canhoto, Ambidestro
  
• TipoSang: A+, A-, …
• ABO: A, B, AB, O
• Rh: +, -
• AtivFisica: nível de atividade física
• Sedentarismo: Não, Sim
  
• MCT: massa corporal total (kg)
• Estatura: cm

Dados <- readxl::read_excel(path="Biometria_FMUSP.xlsx",
                            sheet="dados",
                            na=c("NA","na","nA","Na"))

Dados

tabela <- knitr::kable(head(Dados), format="html", caption="Biometria FMUSP")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

Biometria FMUSP

ID	Ano	Turma	Sexo	Mao	TipoSang	ABO	AtivFisica	Sedentarismo	MCT	Estatura	Rh
1	2	B	M	D	A+	A	baixa_intensidade	N	68	173
2	2	B	M	D	A+	A	atualmente_inativo	S	82	175
3	1	B	M	D	A+	A	media_intensidade	N	79	172
4	1	B	F	D	A+	A	media_intensidade	N	49	160
5	1	B	F	D	O-	O	atualmente_inativo	S	52	164
6	3	A	M	D	A+	A	alta_intensidade	N	73	175

tabela <- knitr::kable(tail(Dados), format="html", caption="Biometria FMUSP")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

Biometria FMUSP

ID	Ano	Turma	Sexo	Mao	TipoSang	ABO	AtivFisica	Sedentarismo	MCT	Estatura	Rh
544	3	B	F	C	NA	NA	alta_intensidade	N	54	157	NA
545	3	B	F	D	B+	B	alta_intensidade	N	47	157
546	3	B	F	D	B+	B	sempre_inativo	S	48	162
547	3	B	M	A	A+	A	alta_intensidade	N	70	175
548	3	B	M	D	O-	O	atualmente_inativo	S	67	171
549	3	B	F	D	A+	A	media_intensidade	N	45	160

Estrutura de dados

str(Dados)

tibble [549 × 12] (S3: tbl_df/tbl/data.frame)
 $ ID          : num [1:549] 1 2 3 4 5 6 7 8 9 10 ...
 $ Ano         : num [1:549] 2 2 1 1 1 3 3 2 3 3 ...
 $ Turma       : chr [1:549] "B" "B" "B" "B" ...
 $ Sexo        : chr [1:549] "M" "M" "M" "F" ...
 $ Mao         : chr [1:549] "D" "D" "D" "D" ...
 $ TipoSang    : chr [1:549] "A+" "A+" "A+" "A+" ...
 $ ABO         : chr [1:549] "A" "A" "A" "A" ...
 $ AtivFisica  : chr [1:549] "baixa_intensidade" "atualmente_inativo" "media_intensidade" "media_intensidade" ...
 $ Sedentarismo: chr [1:549] "N" "S" "N" "N" ...
 $ MCT         : num [1:549] 68 82 79 49 52 73 73 60 75 75 ...
 $ Estatura    : num [1:549] 173 175 172 160 164 175 175 173 182 182 ...
 $ Rh          : chr [1:549] "+" "+" "+" "+" ...

Pré-processamento

Em análise estatística de dados, o pré-processamento refere-se às etapas iniciais e fundamentais que os dados brutos passam antes de serem utilizados para análise estatística propriamente dita. Essas etapas são cruciais para garantir que os dados estejam limpos, organizados e em um formato adequado para que as análises subsequentes sejam precisas e confiáveis.

As principais etapas de pré-processamento incluem:

1. Coleta: Esta é a primeira etapa, onde os dados brutos são obtidos a partir de fontes diversas, como pesquisas, sensores, bancos de dados, etc.

2. Limpeza: Nesta etapa, os dados são revisados em busca de erros, valores ausentes ou inconsistentes. Dados ausentes podem ser tratados preenchendo-os com valores adequados ou removendo as observações afetadas, dependendo do caso. Outliers (valores extremos que não seguem o padrão dos demais dados) podem ser tratados ou removidos, caso tenham um efeito indesejado na análise.

3. Transformação e recodificação: Aqui, os dados podem ser transformados para atender aos requisitos de certos métodos estatísticos. Por exemplo, aplicar logaritmos a dados assimétricos ou normalizar dados para escala comparável.

4. Integração: Quando os dados vêm de várias fontes, é necessário combiná-los em um único conjunto de dados coerente e consistente.

5. Redução de dimensionalidade: Se os dados tiverem muitas variáveis, é possível utilizar técnicas de redução de dimensionalidade para simplificar a análise e visualização.

6. Discretização: Às vezes, é necessário transformar variáveis contínuas em categorias discretas para melhorar a interpretação ou usar certos métodos estatísticos específicos.

7. Normalização e padronização: Essas técnicas são aplicadas para garantir que as variáveis estejam em escalas comparáveis, especialmente quando diferentes variáveis têm unidades diferentes.

8. Seleção de características: Selecione as variáveis mais relevantes ou informativas para a análise, eliminando as menos relevantes.

Após o pré-processamento, os dados estarão prontos para serem utilizados em análises estatísticas, como regressão, análise de variância, análise de componentes principais, clustering, entre outras. O processo de pré-processamento pode variar dependendo da natureza dos dados e das questões de pesquisa específicas em análise estatística.

print(summary(Dados))

       ID           Ano           Turma               Sexo          
 Min.   :  1   Min.   :1.000   Length:549         Length:549        
 1st Qu.:138   1st Qu.:1.000   Class :character   Class :character  
 Median :275   Median :2.000   Mode  :character   Mode  :character  
 Mean   :275   Mean   :2.095                                        
 3rd Qu.:412   3rd Qu.:3.000                                        
 Max.   :549   Max.   :3.000                                        
                                                                    
     Mao              TipoSang             ABO             AtivFisica       
 Length:549         Length:549         Length:549         Length:549        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
 Sedentarismo            MCT            Estatura          Rh           
 Length:549         Min.   : 41.00   Min.   :120.0   Length:549        
 Class :character   1st Qu.: 56.00   1st Qu.:164.0   Class :character  
 Mode  :character   Median : 64.00   Median :172.0   Mode  :character  
                    Mean   : 66.76   Mean   :170.8                     
                    3rd Qu.: 73.00   3rd Qu.:177.0                     
                    Max.   :658.00   Max.   :195.0                     
                    NA's   :6        NA's   :6                         

print(sapply(Dados, class))

          ID          Ano        Turma         Sexo          Mao     TipoSang 
   "numeric"    "numeric"  "character"  "character"  "character"  "character" 
         ABO   AtivFisica Sedentarismo          MCT     Estatura           Rh 
 "character"  "character"  "character"    "numeric"    "numeric"  "character" 

Dados <- dplyr::mutate_if(Dados, is.character, as.factor)
Dados <- dplyr::mutate_at(Dados, dplyr::vars(ID), as.factor)
print(summary(Dados))

       ID           Ano        Turma   Sexo      Mao         TipoSang  
 1      :  1   Min.   :1.000   A:285   F:232   A   : 28   A+     :171  
 2      :  1   1st Qu.:1.000   B:264   M:317   C   : 37   O+     :132  
 3      :  1   Median :2.000                   D   :442   B+     : 65  
 4      :  1   Mean   :2.095                   NA's: 42   A-     : 49  
 5      :  1   3rd Qu.:3.000                              O-     : 45  
 6      :  1   Max.   :3.000                              (Other): 76  
 (Other):543                                              NA's   : 11  
   ABO                   AtivFisica  Sedentarismo      MCT        
 A   :220   alta_intensidade  :168   N:316        Min.   : 41.00  
 AB  : 42   atualmente_inativo:152   S:233        1st Qu.: 56.00  
 B   : 99   baixa_intensidade : 59                Median : 64.00  
 O   :177   media_intensidade : 99                Mean   : 66.76  
 NA's: 11   sempre_inativo    : 71                3rd Qu.: 73.00  
                                                  Max.   :658.00  
                                                  NA's   :6       
    Estatura        Rh     
 Min.   :120.0   -   :130  
 1st Qu.:164.0   +   :408  
 Median :172.0   NA's: 11  
 Mean   :170.8             
 3rd Qu.:177.0             
 Max.   :195.0             
 NA's   :6                 

print(sapply(Dados, class))

          ID          Ano        Turma         Sexo          Mao     TipoSang 
    "factor"    "numeric"     "factor"     "factor"     "factor"     "factor" 
         ABO   AtivFisica Sedentarismo          MCT     Estatura           Rh 
    "factor"     "factor"     "factor"    "numeric"    "numeric"     "factor" 

print(labelled::look_for(Dados[,-1]))

 pos variable     label col_type missing values            
 1   Ano          —     dbl      0                         
 2   Turma        —     fct      0       A                 
                                         B                 
 3   Sexo         —     fct      0       F                 
                                         M                 
 4   Mao          —     fct      42      A                 
                                         C                 
                                         D                 
 5   TipoSang     —     fct      11      A-                
                                         A+                
                                         AB-               
                                         AB+               
                                         B-                
                                         B+                
                                         O-                
                                         O+                
 6   ABO          —     fct      11      A                 
                                         AB                
                                         B                 
                                         O                 
 7   AtivFisica   —     fct      0       alta_intensidade  
                                         atualmente_inativo
                                         baixa_intensidade 
                                         media_intensidade 
                                         sempre_inativo    
 8   Sedentarismo —     fct      0       N                 
                                         S                 
 9   MCT          —     dbl      6                         
 10  Estatura     —     dbl      6                         
 11  Rh           —     fct      11      -                 
                                         +                 

tabela <- knitr::kable(head(Dados), format="html", caption="Biometria FMUSP")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

Biometria FMUSP

ID	Ano	Turma	Sexo	Mao	TipoSang	ABO	AtivFisica	Sedentarismo	MCT	Estatura	Rh
1	2	B	M	D	A+	A	baixa_intensidade	N	68	173
2	2	B	M	D	A+	A	atualmente_inativo	S	82	175
3	1	B	M	D	A+	A	media_intensidade	N	79	172
4	1	B	F	D	A+	A	media_intensidade	N	49	160
5	1	B	F	D	O-	O	atualmente_inativo	S	52	164
6	3	A	M	D	A+	A	alta_intensidade	N	73	175

tabela <- knitr::kable(tail(Dados), format="html", caption="Biometria FMUSP")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

Biometria FMUSP

ID	Ano	Turma	Sexo	Mao	TipoSang	ABO	AtivFisica	Sedentarismo	MCT	Estatura	Rh
544	3	B	F	C	NA	NA	alta_intensidade	N	54	157	NA
545	3	B	F	D	B+	B	alta_intensidade	N	47	157
546	3	B	F	D	B+	B	sempre_inativo	S	48	162
547	3	B	M	A	A+	A	alta_intensidade	N	70	175
548	3	B	M	D	O-	O	atualmente_inativo	S	67	171
549	3	B	F	D	A+	A	media_intensidade	N	45	160

print(class(Dados$AtivFisica))

[1] "factor"

print(levels(Dados$AtivFisica))

[1] "alta_intensidade"   "atualmente_inativo" "baixa_intensidade" 
[4] "media_intensidade"  "sempre_inativo"    

print(summary(subset(Dados, select=c(AtivFisica))))

              AtivFisica 
 alta_intensidade  :168  
 atualmente_inativo:152  
 baixa_intensidade : 59  
 media_intensidade : 99  
 sempre_inativo    : 71  

Dados$AtivFisica <- factor(Dados$AtivFisica,
                           levels=c("sempre_inativo",
                                    "atualmente_inativo",
                                    "baixa_intensidade",
                                    "media_intensidade",
                                    "alta_intensidade"))
print(class(Dados$AtivFisica))

[1] "factor"

print(levels(Dados$AtivFisica))

[1] "sempre_inativo"     "atualmente_inativo" "baixa_intensidade" 
[4] "media_intensidade"  "alta_intensidade"  

print(table(subset(Dados, select=AtivFisica)))

AtivFisica
    sempre_inativo atualmente_inativo  baixa_intensidade  media_intensidade 
                71                152                 59                 99 
  alta_intensidade 
               168 

Dados$ABO <- factor(Dados$ABO,
                   levels=c("A",
                            "O",
                            "B",
                            "AB"))
print(class(Dados$ABO))

[1] "factor"

print(levels(Dados$ABO))

[1] "A"  "O"  "B"  "AB"

print(table(subset(Dados, select=ABO)))

ABO
  A   O   B  AB 
220 177  99  42 

Dados$TipoSang <- factor(Dados$TipoSang,
                           levels=c("A+","O+","B+","AB+",
                                    "A-","O-","B-","AB-"))
print(class(Dados$TipoSang))

[1] "factor"

print(levels(Dados$TipoSang))

[1] "A+"  "O+"  "B+"  "AB+" "A-"  "O-"  "B-"  "AB-"

print(table(subset(Dados, select=TipoSang)))

TipoSang
 A+  O+  B+ AB+  A-  O-  B- AB- 
171 132  65  40  49  45  34   2 

# print(sjPlot::view_df(subset(Dados, select=c(-ID))))
print(summary(subset(Dados, select=-ID)))

      Ano        Turma   Sexo      Mao         TipoSang     ABO     
 Min.   :1.000   A:285   F:232   A   : 28   A+     :171   A   :220  
 1st Qu.:1.000   B:264   M:317   C   : 37   O+     :132   O   :177  
 Median :2.000                   D   :442   B+     : 65   B   : 99  
 Mean   :2.095                   NA's: 42   A-     : 49   AB  : 42  
 3rd Qu.:3.000                              O-     : 45   NA's: 11  
 Max.   :3.000                              (Other): 76             
                                            NA's   : 11             
              AtivFisica  Sedentarismo      MCT            Estatura    
 sempre_inativo    : 71   N:316        Min.   : 41.00   Min.   :120.0  
 atualmente_inativo:152   S:233        1st Qu.: 56.00   1st Qu.:164.0  
 baixa_intensidade : 59                Median : 64.00   Median :172.0  
 media_intensidade : 99                Mean   : 66.76   Mean   :170.8  
 alta_intensidade  :168                3rd Qu.: 73.00   3rd Qu.:177.0  
                                       Max.   :658.00   Max.   :195.0  
                                       NA's   :6        NA's   :6      
    Rh     
 -   :130  
 +   :408  
 NA's: 11  
           
           
           
           

Dados$MCT[Dados$MCT==658] <- NA
Dados$Estatura[Dados$Estatura==120] <- NA
Dados$IMC <- Dados$MCT/((Dados$Estatura/100)^2)
print(summary(subset(Dados, select=c(MCT, Estatura, IMC))))

      MCT            Estatura          IMC       
 Min.   : 41.00   Min.   :153.0   Min.   :14.53  
 1st Qu.: 56.00   1st Qu.:164.0   1st Qu.:20.03  
 Median : 64.00   Median :172.0   Median :21.79  
 Mean   : 65.67   Mean   :170.9   Mean   :22.32  
 3rd Qu.: 73.00   3rd Qu.:177.0   3rd Qu.:24.07  
 Max.   :125.00   Max.   :195.0   Max.   :40.35  
 NA's   :7        NA's   :7       NA's   :8      

# summarytools::view(summarytools::dfSummary(Dados))

saveRDS(Dados, "Biometria_FMUSP.rds")

Estratificação da amostra por sexo

Dados <- readRDS("Biometria_FMUSP.rds")
print(str(Dados))

tibble [549 × 13] (S3: tbl_df/tbl/data.frame)
 $ ID          : Factor w/ 549 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Ano         : num [1:549] 2 2 1 1 1 3 3 2 3 3 ...
 $ Turma       : Factor w/ 2 levels "A","B": 2 2 2 2 2 1 1 2 2 2 ...
 $ Sexo        : Factor w/ 2 levels "F","M": 2 2 2 1 1 2 2 2 2 2 ...
 $ Mao         : Factor w/ 3 levels "A","C","D": 3 3 3 3 3 3 3 3 3 3 ...
 $ TipoSang    : Factor w/ 8 levels "A+","O+","B+",..: 1 1 1 1 6 1 1 3 2 2 ...
 $ ABO         : Factor w/ 4 levels "A","O","B","AB": 1 1 1 1 2 1 1 3 2 2 ...
 $ AtivFisica  : Factor w/ 5 levels "sempre_inativo",..: 3 2 4 4 2 5 5 4 5 5 ...
 $ Sedentarismo: Factor w/ 2 levels "N","S": 1 2 1 1 2 1 1 1 1 1 ...
 $ MCT         : num [1:549] 68 82 79 49 52 73 73 60 75 75 ...
 $ Estatura    : num [1:549] 173 175 172 160 164 175 175 173 182 182 ...
 $ Rh          : Factor w/ 2 levels "-","+": 2 2 2 2 1 2 2 2 2 2 ...
 $ IMC         : num [1:549] 22.7 26.8 26.7 19.1 19.3 ...
NULL

Dados.F <- subset(Dados, Sexo=="F")
print(nrow(Dados.F))

[1] 232

Dados.M <- subset(Dados, Sexo=="M")
print(nrow(Dados.M))

[1] 317

Análise de dados faltantes

n.total <- nrow(Dados)
n.completo <- nrow(na.omit(Dados))
n.incompleto <- n.total - n.completo
cat("Numero de casos total = ", n.total, "\n", sep="")

Numero de casos total = 549

cat("Numero de casos completos = ", n.completo, 
    " (",round(100*n.completo/n.total,2),"%)\n", sep="")

Numero de casos completos = 488 (88.89%)

cat("Numero de casos incompletos = ", n.incompleto, 
    " (",round(100*n.incompleto/n.total,2),"%)\n", sep="")

Numero de casos incompletos = 61 (11.11%)

obs.falt <- sum(is.na(Dados))
obs.valid <- sum(!is.na(Dados))
obs.tot <- obs.falt + obs.valid
cat("Numero de observacoes validas = ", obs.valid, 
    " (",round(100*obs.valid/obs.tot,2),"%)\n", sep="")

Numero de observacoes validas = 7040 (98.64%)

cat("Numero de observacoes faltantes = ", obs.falt, 
    " (",round(100*obs.falt/obs.tot,2),"%)\n", sep="")

Numero de observacoes faltantes = 97 (1.36%)

Análise de dados faltantes com finalfit::missing_pattern

finalfit::missing_pattern(subset(Dados,select=-ID))

    Ano Turma Sexo AtivFisica Sedentarismo MCT Estatura IMC TipoSang ABO Rh Mao
488   1     1    1          1            1   1        1   1        1   1  1   1
42    1     1    1          1            1   1        1   1        1   1  1   0
11    1     1    1          1            1   1        1   1        0   0  0   1
1     1     1    1          1            1   1        0   0        1   1  1   1
1     1     1    1          1            1   0        1   0        1   1  1   1
6     1     1    1          1            1   0        0   0        1   1  1   1
      0     0    0          0            0   7        7   8       11  11 11  42
      
488  0
42   1
11   3
1    2
1    2
6    3
    97

Análise exploratória de dados (EDA)

A análise estatística exploratória (ou EDA - Exploratory Data Analysis) é uma abordagem de investigação de dados que visa entender a estrutura, padrões e características de um conjunto de dados, antes de realizar análises mais formais ou modelagem estatística. Ela é uma etapa essencial no processo de análise de dados, pois ajuda a obter insights preliminares, identificar tendências, detectar anomalias e formular hipóteses iniciais.

As principais características da análise estatística exploratória são:

1. Sumarização de dados: A EDA envolve a utilização de técnicas estatísticas descritivas para resumir os dados em estatísticas como média, mediana, desvio padrão, mínimo, máximo, quartis, histogramas e outras medidas relevantes. Essas estatísticas ajudam a entender a distribuição dos dados e suas principais características.

2. Visualização de dados: A EDA faz amplo uso de gráficos e visualizações para representar os dados de maneira compreensível. Gráficos como histogramas, box plots, gráficos de dispersão, gráficos de linha e gráficos de barras são usados para identificar padrões, relacionamentos e tendências nos dados.

3. Detecção de outliers e dados faltantes: Durante a análise exploratória, é essencial identificar valores extremos (outliers) e dados ausentes, entender suas possíveis causas e decidir como tratá-los, se necessário.

4. Análise de correlação: Explorar correlações entre variáveis ajuda a entender a relação entre diferentes características do conjunto de dados e pode fornecer insights valiosos sobre como elas estão relacionadas.

5. Análise de distribuição: Através de gráficos de distribuição, é possível verificar se os dados seguem uma distribuição específica, como normal, uniforme ou exponencial, por exemplo.

6. Identificação de padrões e tendências: A EDA permite identificar tendências temporais, sazonalidades, padrões cíclicos e outras estruturas nos dados que podem ser úteis na formulação de perguntas e hipóteses para análises posteriores.

A análise estatística exploratória é frequentemente conduzida usando ferramentas de software estatístico, como o R e Python. Essa etapa é fundamental para estabelecer uma base sólida para a análise posterior e para garantir que a interpretação dos resultados seja feita de forma mais informada e precisa.

psych::describe e table

# item name ,item number, nvalid, mean, sd,
# median, mad (Median Absolute Deviation), 
# min, max, skew, kurtosis, se
print(psych::describe(subset(Dados, select=c(MCT, Estatura))))

         vars   n   mean    sd median trimmed   mad min max range skew kurtosis
MCT         1 542  65.67 12.90     64   64.58 13.34  41 125    84 1.04     2.04
Estatura    2 542 170.94  8.92    172  170.80 10.38 153 195    42 0.12    -0.66
           se
MCT      0.55
Estatura 0.38

table(Dados$TipoSang)

 A+  O+  B+ AB+  A-  O-  B- AB- 
171 132  65  40  49  45  34   2 

round(prop.table(table(Dados$TipoSang)),2)

  A+   O+   B+  AB+   A-   O-   B-  AB- 
0.32 0.25 0.12 0.07 0.09 0.08 0.06 0.00 

tapply

# mean,median,25th and 75th quartiles,min,max
tapply(Dados$MCT, Dados$Sexo, summary)

$F
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  41.00   52.00   56.00   57.63   62.00  120.00       2 

$M
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  43.00   63.00   72.00   71.59   78.00  125.00       5 

psych::describeBy

# item name, item number, nvalid, mean, sd,
# median, mad (Median Absolute Deviation), 
# min, max, skew, kurtosis, se
print(psych::describeBy(subset(Dados,select=c(MCT, Estatura)),
                        group=list(Dados$Sexo),
                        mat=TRUE,
                        digits=2))

          item group1 vars   n   mean    sd median trimmed   mad min max range
MCT1         1      F    1 230  57.63  9.05     56   56.75  7.41  41 120    79
MCT2         2      M    1 312  71.59 12.09     72   70.88 11.12  43 125    82
Estatura1    3      F    2 229 164.34  6.38    163  163.94  5.93 153 186    33
Estatura2    4      M    2 313 175.77  7.26    175  175.88  7.41 155 195    40
           skew kurtosis   se
MCT1       2.43    12.04 0.60
MCT2       0.90     2.37 0.68
Estatura1  0.61     0.02 0.42
Estatura2 -0.13     0.21 0.41

print(psych::describeBy(subset(Dados,select=c(MCT, Estatura)),
                        group=list(Dados$Sexo, Dados$ABO),
                        mat=TRUE,
                        digits=2))

          item group1 group2 vars   n   mean    sd median trimmed   mad min max
MCT1         1      F      A    1  89  57.79  9.32   57.0   57.12  7.41  43 120
MCT2         2      M      A    1 126  73.21 13.06   72.0   72.16 10.38  50 125
MCT3         3      F      O    1  81  59.10 10.07   57.0   57.98  7.41  41 100
MCT4         4      M      O    1  96  71.38 11.48   73.0   71.12 10.38  43 120
MCT5         5      F      B    1  42  55.26  7.02   53.0   54.44  5.19  45  76
MCT6         6      M      B    1  56  68.38 10.48   65.0   67.61  7.41  49  95
MCT7         7      F     AB    1  14  56.93  5.69   56.5   57.00  6.67  48  65
MCT8         8      M     AB    1  27  71.93 11.33   72.0   71.65  8.90  50 102
Estatura1    9      F      A    2  88 164.59  5.83  164.0  164.39  5.93 153 181
Estatura2   10      M      A    2 126 175.17  6.79  175.0  175.08  7.41 160 192
Estatura3   11      F      O    2  81 164.33  6.90  163.0  163.80  7.41 154 186
Estatura4   12      M      O    2  96 176.98  7.12  177.0  177.14  7.41 156 193
Estatura5   13      F      B    2  42 163.43  6.84  162.5  162.76  6.67 153 185
Estatura6   14      M      B    2  57 175.26  7.73  175.0  175.23  5.93 160 195
Estatura7   15      F     AB    2  14 165.21  5.67  164.0  165.08  5.93 156 176
Estatura8   16      M     AB    2  27 175.41  8.95  177.0  175.83  8.90 155 195
          range  skew kurtosis   se
MCT1         77  3.29    19.99 0.99
MCT2         75  1.12     2.65 1.16
MCT3         59  1.69     4.67 1.12
MCT4         77  0.58     2.26 1.17
MCT5         31  1.03     0.52 1.08
MCT6         46  0.69     0.08 1.40
MCT7         17  0.06    -1.50 1.52
MCT8         52  0.47     0.64 2.18
Estatura1    28  0.39    -0.53 0.62
Estatura2    32  0.09    -0.35 0.61
Estatura3    32  0.69    -0.06 0.77
Estatura4    37 -0.36     0.81 0.73
Estatura5    32  0.91     0.72 1.06
Estatura6    35  0.04    -0.23 1.02
Estatura7    20  0.25    -1.05 1.52
Estatura8    40 -0.35     0.23 1.72

print(psych::describeBy(subset(Dados, select=c(MCT, Estatura)),
                        group=list(Dados$ABO, Dados$Sexo),
                        mat=TRUE,
                        digits=2))

          item group1 group2 vars   n   mean    sd median trimmed   mad min max
MCT1         1      A      F    1  89  57.79  9.32   57.0   57.12  7.41  43 120
MCT2         2      O      F    1  81  59.10 10.07   57.0   57.98  7.41  41 100
MCT3         3      B      F    1  42  55.26  7.02   53.0   54.44  5.19  45  76
MCT4         4     AB      F    1  14  56.93  5.69   56.5   57.00  6.67  48  65
MCT5         5      A      M    1 126  73.21 13.06   72.0   72.16 10.38  50 125
MCT6         6      O      M    1  96  71.38 11.48   73.0   71.12 10.38  43 120
MCT7         7      B      M    1  56  68.38 10.48   65.0   67.61  7.41  49  95
MCT8         8     AB      M    1  27  71.93 11.33   72.0   71.65  8.90  50 102
Estatura1    9      A      F    2  88 164.59  5.83  164.0  164.39  5.93 153 181
Estatura2   10      O      F    2  81 164.33  6.90  163.0  163.80  7.41 154 186
Estatura3   11      B      F    2  42 163.43  6.84  162.5  162.76  6.67 153 185
Estatura4   12     AB      F    2  14 165.21  5.67  164.0  165.08  5.93 156 176
Estatura5   13      A      M    2 126 175.17  6.79  175.0  175.08  7.41 160 192
Estatura6   14      O      M    2  96 176.98  7.12  177.0  177.14  7.41 156 193
Estatura7   15      B      M    2  57 175.26  7.73  175.0  175.23  5.93 160 195
Estatura8   16     AB      M    2  27 175.41  8.95  177.0  175.83  8.90 155 195
          range  skew kurtosis   se
MCT1         77  3.29    19.99 0.99
MCT2         59  1.69     4.67 1.12
MCT3         31  1.03     0.52 1.08
MCT4         17  0.06    -1.50 1.52
MCT5         75  1.12     2.65 1.16
MCT6         77  0.58     2.26 1.17
MCT7         46  0.69     0.08 1.40
MCT8         52  0.47     0.64 2.18
Estatura1    28  0.39    -0.53 0.62
Estatura2    32  0.69    -0.06 0.77
Estatura3    32  0.91     0.72 1.06
Estatura4    20  0.25    -1.05 1.52
Estatura5    32  0.09    -0.35 0.61
Estatura6    37 -0.36     0.81 0.73
Estatura7    35  0.04    -0.23 1.02
Estatura8    40 -0.35     0.23 1.72

aggregate

print(aggregate(subset(Dados, select=c(MCT, Estatura)),
                by=list(Dados$Sexo),
                FUN=mean,
                na.rm=TRUE))

  Group.1      MCT Estatura
1       F 57.63043 164.3406
2       M 71.58974 175.7732

print(aggregate(subset(Dados, select=c(MCT, Estatura)),
                by=list(Dados$Sexo, Dados$ABO),
                FUN=mean,
                na.rm=TRUE))

  Group.1 Group.2      MCT Estatura
1       F       A 57.78652 164.5909
2       M       A 73.20635 175.1746
3       F       O 59.09877 164.3333
4       M       O 71.37500 176.9792
5       F       B 55.26190 163.4286
6       M       B 68.37500 175.2632
7       F      AB 56.92857 165.2143
8       M      AB 71.92593 175.4074

print(aggregate(subset(Dados, select=c(MCT, Estatura)),
                by=list(Dados$ABO, Dados$Sexo),
                FUN=mean,
                na.rm=TRUE))

  Group.1 Group.2      MCT Estatura
1       A       F 57.78652 164.5909
2       O       F 59.09877 164.3333
3       B       F 55.26190 163.4286
4      AB       F 56.92857 165.2143
5       A       M 73.20635 175.1746
6       O       M 71.37500 176.9792
7       B       M 68.37500 175.2632
8      AB       M 71.92593 175.4074

aggregate por fórmula

print(aggregate(MCT ~ Sexo,
                FUN=median,
                na.rm=TRUE,
                data=Dados))

  Sexo MCT
1    F  56
2    M  72

print(aggregate(MCT ~ Sexo + ABO,
                FUN=median,
                na.rm=TRUE,
                data=Dados))

  Sexo ABO  MCT
1    F   A 57.0
2    M   A 72.0
3    F   O 57.0
4    M   O 73.0
5    F   B 53.0
6    M   B 65.0
7    F  AB 56.5
8    M  AB 72.0

Média e desvio-padrão: Parte A

Baseado em

• Higgins JPT, Li T, Deeks JJ (editors). Chapter 6: Choosing effect measures and computing estimates of effect. In: Higgins JPT, Thomas J, Chandler J, Cumpston M, Li T, Page MJ, Welch VA (editors). Cochrane Handbook for Systematic Reviews of Interventions version 6.3 (updated February 2022). Cochrane, 2022. Available from https:www.training.cochrane.org/handbook.

• ver Table 6.5.a em https://training.cochrane.org/handbook/current/chapter-06

• Suponha que não temos acesso aos dados brutos de MCT. Os dados disponíveis são apresentados na tabela abaixo:

group	n	mean	sd
F	230	57.63	9.05
M	312	71.59	12.09

• Como temos acesso aos dados brutos de MCT do arquivo Biometria_FMUSP.rds podemos obter a média (função mean) e o desvio-padrão (função sd, de standard deviation) dos indivíduos de cada um dos dois sexos que geraram essa tabela:

Dados <- readRDS("Biometria_FMUSP.rds")

Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")

n.F <- sum(!is.na(Dados.F$MCT))
mean.F <- mean(Dados.F$MCT,na.rm=TRUE)
sd.F <- sd(Dados.F$MCT,na.rm=TRUE)
n.M <- sum(!is.na(Dados.M$MCT))
mean.M <- mean(Dados.M$MCT,na.rm=TRUE)
sd.M <- sd(Dados.M$MCT,na.rm=TRUE)
cat("Sexo F: média=",mean.F,", desvio-padrão=",sd.F," (n=",n.F,")\n",sep="")

Sexo F: média=57.63043, desvio-padrão=9.052874 (n=230)

cat("Sexo M: média=",mean.M,", desvio-padrão=",sd.M," (n=",n.M,")\n",sep="")

Sexo M: média=71.58974, desvio-padrão=12.08723 (n=312)

• Quais são os valores da média e desvio-padrão amostrais dos 542 estudantes?

A média global é facilmente obtida fazendo-se uma média ponderada

m <- (230*57.63 + 312*71.59)/542
cat("Média = ", round(m,2), "\n", sep="")

Média = 65.67

Como temos os dados brutos, podemos conferir nosso resultado:

print(mean.All <- mean(Dados$MCT,na.rm=TRUE))

[1] 65.66605

Podemos fazer o mesmo para o desvio-padrão?

s <- (230*9.05 + 312*12.09)/542
cat("Desvio-padrão = ", round(s,2), "\n", sep="")

Desvio-padrão = 10.8

Conferindo nosso resultado, obtemos:

print(sd.All <- sd(Dados$MCT,na.rm=TRUE))

[1] 12.89833

Onde está o erro?

Repare que 12.9 é maior do que os dois desvios-padrão dos grupos (F: 9.05; M: 12.09).

Média e desvio-padrão: Parte B

• MCT

  n	mean	sd
  542	65.67	12.9

group	n	mean	sd
F	230	57.63	9.05
M	312	71.59	12.09

A média foi bem resolvida:

\[\Large \bar{x} = \dfrac{n_F \, \bar{x}_F + n_M \, \bar{x}_M}{n_F + n_M}\]

implementada em R como:

m <- (230*57.63 + 312*71.59)/542
cat("Media = ", round(m,2), "\n", sep="")

Media = 65.67

O desvio-padrão é mais complexo. Há algumas maneiras equivalentes para obtê-lo quando não dispomos dos dados brutos:

\[\Large s=\sqrt{\dfrac{(n_F-1)\,s_F^2+(n_M-1)\,s_M^2+\dfrac{n_F\;n_M}{n_F+n_M}\left(\bar{x}_F-\bar{x}_M\right)^2}{n_F+n_M-1}}\\ \large \dfrac{n_F\;n_M}{n_F+n_M}=\text{ média harmônica/2}\]

\[\Large s=\sqrt{\frac{n_F\,s_F^2+n_M\,s_M^2+n_F\,\left(\bar{x}_F-\bar{x}\right)^2+n_M\,\left(\bar{x}_M-\bar{x}\right)^2}{n_F+n_M}}\]

\[\Large s=\sqrt{\frac{n_F\,\left(s_F^2+\left(\bar{x}_F-\bar{x}\right)^2\right)+n_M\,\left(s_M^2+\left(\bar{x}_M-\bar{x}\right)^2\right)}{n_F+n_M}}\]

implementados em R, respectivamente, como:

dp <- sqrt(((230-1)*9.05^2 + (312-1)*12.09^2 + 
            (230*312/542)*(57.63-71.59)^2)/(542-1))
cat("Desvio-padrao = ", round(dp,2), "\n", sep="")

Desvio-padrao = 12.9

dp <- sqrt((230*9.05^2 + 312*12.09^2 + 
           230*(57.63-65.67)^2 + 312*(71.59-65.67)^2)/542)
cat("Desvio-padrao = ", round(dp,2), "\n", sep="")

Desvio-padrao = 12.9

dp <- sqrt( (230*(9.05^2+(57.63-65.67)^2)+
             312*(12.09^2+(71.59-65.67)^2)) / 542)
cat("Desvio-padrao = ", round(dp,2), "\n", sep="")

Desvio-padrao = 12.9

Como pode o desvio-padrão global ser maior que os de cada grupo? Podemos ver graficamente:

# obtem os dados brutos
Dados <- readRDS("Biometria_FMUSP.rds")
Dados.F <- subset(Dados, Sexo=="F")
Dados.M <- subset(Dados, Sexo=="M")
# tamanho dos grupos
n.All <- sum(!is.na(Dados$MCT))
n.F <- sum(!is.na(Dados.F$MCT))
n.M <- sum(!is.na(Dados.M$MCT))
# distribuicao densidade de probabilidade
d.All <- density(Dados$MCT,na.rm=TRUE)
d.F <- density(Dados.F$MCT,na.rm=TRUE)
d.M <- density(Dados.M$MCT,na.rm=TRUE)
# escala a densidade pelo tamanho dos grupos
d.All$y <- d.All$y * n.All
d.F$y <- d.F$y * n.F
d.M$y <- d.M$y * n.M
# exibe o grafico
x.min <- min(c(d.All$x,d.F$x,d.M$x),na.rm=TRUE)
x.max <- max(c(d.All$x,d.F$x,d.M$x),na.rm=TRUE)
y.min <- min(c(d.All$y,d.F$y,d.M$y),na.rm=TRUE)
y.max <- max(c(d.All$y,d.F$y,d.M$y),na.rm=TRUE)
plot(NA, main="Distribuição do MCT",
     xlab="MCT", ylab="Densidade * n",
     xlim=c(x.min,x.max), ylim=c(y.min,y.max),
     bty="n")
lines(d.All, lwd=2)
lines(d.F, col="red", lwd=2, lty=2)
lines(d.M, col="blue", lwd=2, lty=3)
legend("topright", 
       c("Global","Feminino","Masculino"), 
       col=c("black","red","blue"),
       lwd=2, 
       lty=c(1,2,3), 
       box.lwd=0, bg="transparent")  

Gráficos

Dotplot por plot & table

plot(Dados$Sexo, xlab="Sexo", ylab="Freq")

plot(Dados$TipoSang, xlab="Tipo Sanguineo", ylab="Freq")

plot(Dados$Sedentarismo, xlab="Sedentarismo", ylab="Freq")

plot(table(Dados$MCT), xlab="Massa corporal total (kg)", ylab="Freq")

plot(table(Dados$Estatura), xlab="Estatura (cm)", ylab="Freq")

Gráfico de setores

table(Dados$Sexo)

  F   M 
232 317 

pie(table(Dados$Sexo), 
    xlab="Sexo")

table(Dados$ABO)

  A   O   B  AB 
220 177  99  42 

pie(table(Dados$ABO), 
    xlab="ABO")

Tabela de contingência 2x2

print(Sedentarismo.Sexo <- xtabs(~Sedentarismo+Sexo, data=Dados))

            Sexo
Sedentarismo   F   M
           N 133 183
           S  99 134

# total por linha (1)
margin.table(Sedentarismo.Sexo,1)

Sedentarismo
  N   S 
316 233 

# total por coluna (2)
margin.table(Sedentarismo.Sexo,2)

Sexo
  F   M 
232 317 

# proporção de cada célula em relação ao total
round(proportions(Sedentarismo.Sexo),3)

            Sexo
Sedentarismo     F     M
           N 0.242 0.333
           S 0.180 0.244

# proporção por linha
round(proportions(Sedentarismo.Sexo,1),3)

            Sexo
Sedentarismo     F     M
           N 0.421 0.579
           S 0.425 0.575

# proporção por coluna
round(proportions(Sedentarismo.Sexo,2),3)

            Sexo
Sedentarismo     F     M
           N 0.573 0.577
           S 0.427 0.423

plot(Dados$Sedentarismo~Dados$Sexo, xlab="Sexo", ylab="Sedentarismo")

mosaicplot(~Sexo+Sedentarismo, data=Dados, color=FALSE)

gplots::balloonplot(t(Sedentarismo.Sexo), 
                    main ="Estudantes de Medicina\nFMUSP", 
                    xlab ="Sexo", 
                    ylab="Sedentarismo",
                    label=TRUE, 
                    show.margins=TRUE, 
                    show.zeros=TRUE, 
                    dotcolor="gray45")

barplot(Sedentarismo.Sexo,
        beside=TRUE, 
        legend.text=rownames(Sedentarismo.Sexo),
        ylab="Freq",
        xlab="Sexo x Sedentarismo")

barplot(proportions(Sedentarismo.Sexo),
        beside=TRUE, 
        legend.text=rownames(Sedentarismo.Sexo),
        ylab="Freq",
        xlab="Sexo x Sedentarismo")

sexo.ABO.freq <- as.data.frame(table(Dados$Sexo, Dados$ABO))
names(sexo.ABO.freq) <- c("Sexo", "ABO", "Freq")
ggpubr::ggbarplot(sexo.ABO.freq, 
                  x="ABO", 
                  y="Freq",
                  color="Sexo",
                  palette=c("gray45", "black"),
                  order=c("A", "O", "B", "AB"),
                  width=.7)

ggpubr::ggbarplot(sexo.ABO.freq, 
                  x="ABO", 
                  y="Freq",
                  color="Sexo",
                  palette=c("gray45", "black"),
                  order=c("A", "O", "B", "AB"),
                  position = ggplot2::position_dodge(),
                  width=.7)

Tabela de contingência 2x4

print(ABO.Sexo <- xtabs(~ABO+Sexo, data=Dados))

    Sexo
ABO    F   M
  A   91 129
  O   81  96
  B   42  57
  AB  14  28

round(proportions(ABO.Sexo),2)

    Sexo
ABO     F    M
  A  0.17 0.24
  O  0.15 0.18
  B  0.08 0.11
  AB 0.03 0.05

round(proportions(ABO.Sexo,1),2)

    Sexo
ABO     F    M
  A  0.41 0.59
  O  0.46 0.54
  B  0.42 0.58
  AB 0.33 0.67

round(proportions(ABO.Sexo,2),2)

    Sexo
ABO     F    M
  A  0.40 0.42
  O  0.36 0.31
  B  0.18 0.18
  AB 0.06 0.09

mosaicplot(~Sexo+ABO, data=Dados, color=FALSE)

gplots::balloonplot(t(ABO.Sexo), 
                    main ="Estudantes de Medicina\nFMUSP", 
                    xlab ="Sexo", 
                    ylab="ABO",
                    label=TRUE, show.margins=TRUE, 
                    show.zeros=TRUE, 
                    dotcolor="gray45")

barplot(ABO.Sexo,
        beside=TRUE, 
        legend.text=rownames(ABO.Sexo),
        ylab="Freq",
        xlab="Sexo x ABO")

barplot(proportions(ABO.Sexo),
        beside=TRUE, 
        legend.text=rownames(ABO.Sexo),
        ylab="Freq",
        xlab="Sexo x ABO")

Tabelas multidimensionais: três ou mais variáveis categóricas

ftable(Sexo + Sedentarismo ~ ABO, data=Dados)

    Sexo          F     M   
    Sedentarismo  N  S  N  S
ABO                         
A                51 40 83 46
O                45 36 56 40
B                26 16 27 30
AB                9  5 12 16

print(xtabs(~Sexo + Sedentarismo + ABO, data=Dados))

, , ABO = A

    Sedentarismo
Sexo  N  S
   F 51 40
   M 83 46

, , ABO = O

    Sedentarismo
Sexo  N  S
   F 45 36
   M 56 40

, , ABO = B

    Sedentarismo
Sexo  N  S
   F 26 16
   M 27 30

, , ABO = AB

    Sedentarismo
Sexo  N  S
   F  9  5
   M 12 16

mosaicplot(xtabs(~Sexo + Sedentarismo + ABO, data=Dados))

Boxplot

boxplot(Dados$MCT, horizontal=TRUE, 
        xlab="MCT (kg)")
rug(jitter(Dados$MCT))

boxplot(MCT~Sexo, data=Dados, horizontal=TRUE, 
        xlab="MCT (kg)")

boxplot(MCT~Sexo+Sedentarismo, data=Dados, horizontal=TRUE, 
        xlab="MCT (kg)")

boxplot(MCT~Sexo+Sedentarismo+Mao, data=Dados, horizontal=TRUE, 
        xlab="MCT (kg)", cex=0.6)

ggpubr::ggboxplot(data=Dados, 
                  y="MCT")

ggpubr::ggboxplot(data=Dados,
                  x="Sexo",
                  y="MCT", 
                  add="",
                  orientation="horizontal",
                  width=.7,
                  order=c("F", "M"))

ggpubr::ggboxplot(data=Dados,
                  x="Sexo",
                  y="MCT", 
                  add="",
                  orientation="horizontal",
                  width=.7,
                  order=c("M", "F"))

ggpubr::ggboxplot(data=Dados,
                  x="Sexo",
                  y="MCT", 
                  add="jitter",
                  orientation="horizontal",
                  width=.7,
                  order=c("F", "M"))

ggpubr::ggboxplot(data=Dados,
                  x="ABO",
                  y="MCT", 
                  add="",
                  orientation="horizontal",
                  width=.7,
                  select=c("A", "B", "O"),
                  order=c("O", "B", "A"))

Bagplot: boxplot bidimensional

bgp.F <- DescTools::PlotBag(Dados.F$Estatura, 
                            Dados.F$MCT,
                            main=paste("Feminino"),
                            xlab="Estatura (cm)",
                            ylab="Massa Corporal Total (kg)",
                            na.rm=TRUE,
                            show.bagpoints=FALSE,
                            show.looppoints=FALSE,
                            show.whiskers=FALSE,
                            show.baghull=TRUE,
                            col.loophull="white",
                            col.looppoints="black", 
                            col.baghull="white",
                            col.bagpoints="black",
                            add=FALSE,
                            cex=1)

[1] "Warning: NA elements have been removed!!"

print(outliers.F <- as.data.frame(bgp.F$pxy.outlier))

    x   y
1 165 100
2 165 100
3 158  78

for (o in 1:nrow(outliers.F))
{
  r.F <- which(Dados.F$Estatura==outliers.F$x[o] & 
               Dados.F$MCT==outliers.F$y[o])
  text(outliers.F$x[o],outliers.F$y[o], r.F, pos=1, cex=0.7)
}

bgp.M <- DescTools::PlotBag(Dados.M$Estatura, 
                          Dados.M$MCT,
                          main=paste("Masculino"),
                          xlab="Estatura (cm)",
                          ylab="Massa Corporal Total (kg)",
                          na.rm=TRUE,
                          show.bagpoints=FALSE,
                          show.looppoints=FALSE,
                          show.whiskers=FALSE,
                          show.baghull=TRUE,
                          col.loophull="white",
                          col.looppoints="black", 
                          col.baghull="white",
                          col.bagpoints="black",
                          add=FALSE,
                          cex=1)

[1] "Warning: NA elements have been removed!!"

print(outliers.M <- as.data.frame(bgp.M$pxy.outlier))

    x   y
1 165  94
2 193 120
3 184 125
4 176 125

for (o in 1:nrow(outliers.M))
{
  r.M <- which(Dados.M$Estatura==outliers.M$x[o] & 
                 Dados.M$MCT==outliers.M$y[o])
  text(outliers.M$x[o], outliers.M$y[o], r.M, pos=1, cex=0.7)
}

Histograma

• Silveira, PSP & Siqueira, JO (2022) Histogram lies about distribution shape and Pearson’s coefficient of variation lies about relative variability. The Quantitative Methods for Psychology 18(1). DOI 10.20982/tqmp.18.1.p091

• How to Lie with Histograms

  • https://www.rossidata.com/LinLog

hist(Dados$MCT)
rug(jitter(Dados$MCT))

Gráfico de densidade

car::densityPlot(~MCT, data=Dados)

car::densityPlot(MCT~Sexo, data=Dados, 
                 col=c("black","black"))

Note que a área de cada density plot é idêntica (igual a 1). Caso queiramos refletir o tamanho de cada grupo, uma alternativa é multiplicar os valores das ordenadas pelo tamanho do grupo (a área ficará igual ao tamanho do grupo).

# tamanho dos grupos
n.F <- sum(!is.na(Dados.F$MCT))
n.M <- sum(!is.na(Dados.M$MCT))
# distribuicao densidade de probabilidade
d.F <- density(Dados.F$MCT,na.rm=TRUE)
d.M <- density(Dados.M$MCT,na.rm=TRUE)
# escala a densidade pelo tamanho dos grupos
d.F$y <- d.F$y * n.F
d.M$y <- d.M$y * n.M
# exibe o grafico
x.min <- min(c(d.F$x,d.M$x),na.rm=TRUE)
x.max <- max(c(d.F$x,d.M$x),na.rm=TRUE)
y.min <- min(c(d.F$y,d.M$y),na.rm=TRUE)
y.max <- max(c(d.F$y,d.M$y),na.rm=TRUE)
plot(NA, main="Distribuição do MCT",
     xlab="MCT", ylab="Densidade * n",
     xlim=c(x.min,x.max), ylim=c(y.min,y.max),
     bty="n")
lines(d.F, col="violetred3", lwd=2, lty=1)
lines(d.M, col="dodgerblue2", lwd=3, lty=2)
legend("topright", 
       c("Feminino","Masculino"), 
       col=c("violetred3","dodgerblue2"),
       lwd=c(2,3), 
       lty=c(1,2), 
       box.lwd=0, bg="transparent")  

car::densityPlot(MCT~Sedentarismo, data=Dados.F, 
                 col=c("black","black"),
                 main="Feminino")

car::densityPlot(MCT~Sedentarismo, data=Dados.M, 
                 col=c("black","black"),
                 main="Masculino")

car::densityPlot(MCT~ABO, data=Dados.F,
                 main="Feminino")

car::densityPlot(MCT~ABO, data=Dados.M,
                 main="Masculino")

Gráfico de distribuição de probabilidade: ECDF plot

ggpubr::ggecdf(data=Dados,
               x="MCT",
               linetype="Sexo")

Gráfico de dispersão com sunflowerplot

Considere um scatterplot convencional:

pairs.total <- nrow(na.omit(data.frame(Dados.F[,c("MCT","Estatura")])))
cat("Existem ",pairs.total," pares de valores (Estatura,MCT) em Dados.F\n",sep="")

Existem 229 pares de valores (Estatura,MCT) em Dados.F

plot(MCT~Estatura,
     main="Estatura x MCT: Feminino",
     data=Dados.F, 
     pch=19,
     cex=0.6,
     col="black", 
     xlim=c(145,190),
     ylim=c(30,127),
     bty="n")
uniestaturas <- sort(unique(Dados.F$Estatura))
points.per.mct <- c()
for (e in uniestaturas)
{
  m <- sort(unique(Dados.F$MCT[Dados.F$Estatura==e]))
  text(e,m[1],length(m),pos=1,col="blue",cex=0.8)
  points.per.mct <- c(points.per.mct,length(m))
}

points.total <- sum(points.per.mct)
cat("Contamos ",points.total," pontos distintos no gráfico.\n",sep="")

Contamos 156 pontos distintos no gráfico.

Para facilitar, usei um loop com for para contar quantos valores distintos de MCT aparecem associados a cada Estatura, e coloquei os números em azul no próprio gráfico. Você pode conferir (se tiver paciência) que

2+1+3+9+7+5+3+16+4+7+10+8+12+7+9+6+5+8+2+8+6+3+5+3+3+1+1+1+1 = 156.

Os 229-156=73 pontos faltantes correspondem a indivíduos com pares de valores (Estatura,MCT) coincidentes, mas os deixamos de visualizar graficamente. Há duas maneiras para contornar este problema. O sunflowerplot é uma delas:

sunflowerplot(MCT~Estatura,
              main="Estatura x MCT: Feminino",
              data=Dados.F, 
              rotate=TRUE, 
              pch=1,
              size=.1,
              col="black", 
              seg.col="black",
              xlim=c(145,190),
              ylim=c(30,127),
              seg.lwd=.8)

sunflowerplot(MCT~Estatura,
              main="Estatura x MCT: Masculino",
              data=Dados.M, 
              rotate=TRUE, 
              pch=2,
              size=.1,
              col="black", 
              seg.col="black", 
              xlim=c(150,200),
              ylim=c(30,127),
              seg.lwd=.8)

sunflowerplot(MCT~Estatura,
              main="Estatura x MCT: Feminino & Masculino",
              data=Dados.F, 
              rotate=TRUE, 
              pch=1,
              size=.1,
              col="black", 
              seg.col="black", 
              xlim=c(145,200),
              ylim=c(30,127),
              seg.lwd=.8)
sunflowerplot(MCT~Estatura,
              data=Dados.M, 
              rotate=TRUE, 
              pch=2,
              size=.1,
              col="black", 
              seg.col="black", 
              seg.lwd=.8,
              add=TRUE)

Gráfico de dispersão com car::scatterplot

car::scatterplot(MCT~Estatura,
                 group=Dados$Sexo,
                 jitter=list(x=1, y=1),
                 regLine=FALSE, 
                 smooth=FALSE, 
                 ellipse=FALSE,
                 grid=FALSE,
                 col="black",
                 xlim=c(145,200),
                 ylim=c(30,127),
                 data=Dados)

Uma palavra sobre cores Observe que aparece um warning number of groups exceeds number of available colors colors are recycled Isto é causado por col=“black”, uma única cor, mas temos dois grupos, Feminino e Masculino. Se quiséssemos tudo em preto, como está, isto seria resolvido com col=c("black","black") No entanto, R pode usar todos os \(256^3=16777216\) (mais que 16 milhões de cores) disponíveis no sistema RGB (“red”, “green”, “blue”) utilizando a notação hexadecimal. Há 256 níveis de vermelho, verde ou azul, indo da ausência (0) à intensidade máxima (255, denotado como ff em hexadecimal). Preto é ausência das três cores e corresponde a “#000000” na sintaxe do R, que usa os três pares de valores para definir as intensidade com “#rrggbb”. Além de você poder colocar a cor que quiser, há algumas centenas de cores já disponíveis no R. Consulte ?colors para mais detalhes. Esta função fornece o seguinte: colors() [1] "white" "aliceblue" "antiquewhite" [4] "antiquewhite1" "antiquewhite2" "antiquewhite3" [7] "antiquewhite4" "aquamarine" "aquamarine1" [10] "aquamarine2" "aquamarine3" "aquamarine4" [13] "azure" "azure1" "azure2" [16] "azure3" "azure4" "beige" [19] "bisque" "bisque1" "bisque2" [22] "bisque3" "bisque4" "black" [25] "blanchedalmond" "blue" "blue1" [28] "blue2" "blue3" "blue4" [31] "blueviolet" "brown" "brown1" [34] "brown2" "brown3" "brown4" [37] "burlywood" "burlywood1" "burlywood2" [40] "burlywood3" "burlywood4" "cadetblue" [43] "cadetblue1" "cadetblue2" "cadetblue3" [46] "cadetblue4" "chartreuse" "chartreuse1" [49] "chartreuse2" "chartreuse3" "chartreuse4" [52] "chocolate" "chocolate1" "chocolate2" [55] "chocolate3" "chocolate4" "coral" [58] "coral1" "coral2" "coral3" [61] "coral4" "cornflowerblue" "cornsilk" [64] "cornsilk1" "cornsilk2" "cornsilk3" [67] "cornsilk4" "cyan" "cyan1" [70] "cyan2" "cyan3" "cyan4" [73] "darkblue" "darkcyan" "darkgoldenrod" [76] "darkgoldenrod1" "darkgoldenrod2" "darkgoldenrod3" [79] "darkgoldenrod4" "darkgray" "darkgreen" [82] "darkgrey" "darkkhaki" "darkmagenta" [85] "darkolivegreen" "darkolivegreen1" "darkolivegreen2" [88] "darkolivegreen3" "darkolivegreen4" "darkorange" [91] "darkorange1" "darkorange2" "darkorange3" [94] "darkorange4" "darkorchid" "darkorchid1" [97] "darkorchid2" "darkorchid3" "darkorchid4" [100] "darkred" "darksalmon" "darkseagreen" [103] "darkseagreen1" "darkseagreen2" "darkseagreen3" [106] "darkseagreen4" "darkslateblue" "darkslategray" [109] "darkslategray1" "darkslategray2" "darkslategray3" [112] "darkslategray4" "darkslategrey" "darkturquoise" [115] "darkviolet" "deeppink" "deeppink1" [118] "deeppink2" "deeppink3" "deeppink4" [121] "deepskyblue" "deepskyblue1" "deepskyblue2" [124] "deepskyblue3" "deepskyblue4" "dimgray" [127] "dimgrey" "dodgerblue" "dodgerblue1" [130] "dodgerblue2" "dodgerblue3" "dodgerblue4" [133] "firebrick" "firebrick1" "firebrick2" [136] "firebrick3" "firebrick4" "floralwhite" [139] "forestgreen" "gainsboro" "ghostwhite" [142] "gold" "gold1" "gold2" [145] "gold3" "gold4" "goldenrod" [148] "goldenrod1" "goldenrod2" "goldenrod3" [151] "goldenrod4" "gray" "gray0" [154] "gray1" "gray2" "gray3" [157] "gray4" "gray5" "gray6" [160] "gray7" "gray8" "gray9" [163] "gray10" "gray11" "gray12" [166] "gray13" "gray14" "gray15" [169] "gray16" "gray17" "gray18" [172] "gray19" "gray20" "gray21" [175] "gray22" "gray23" "gray24" [178] "gray25" "gray26" "gray27" [181] "gray28" "gray29" "gray30" [184] "gray31" "gray32" "gray33" [187] "gray34" "gray35" "gray36" [190] "gray37" "gray38" "gray39" [193] "gray40" "gray41" "gray42" [196] "gray43" "gray44" "gray45" [199] "gray46" "gray47" "gray48" [202] "gray49" "gray50" "gray51" [205] "gray52" "gray53" "gray54" [208] "gray55" "gray56" "gray57" [211] "gray58" "gray59" "gray60" [214] "gray61" "gray62" "gray63" [217] "gray64" "gray65" "gray66" [220] "gray67" "gray68" "gray69" [223] "gray70" "gray71" "gray72" [226] "gray73" "gray74" "gray75" [229] "gray76" "gray77" "gray78" [232] "gray79" "gray80" "gray81" [235] "gray82" "gray83" "gray84" [238] "gray85" "gray86" "gray87" [241] "gray88" "gray89" "gray90" [244] "gray91" "gray92" "gray93" [247] "gray94" "gray95" "gray96" [250] "gray97" "gray98" "gray99" [253] "gray100" "green" "green1" [256] "green2" "green3" "green4" [259] "greenyellow" "grey" "grey0" [262] "grey1" "grey2" "grey3" [265] "grey4" "grey5" "grey6" [268] "grey7" "grey8" "grey9" [271] "grey10" "grey11" "grey12" [274] "grey13" "grey14" "grey15" [277] "grey16" "grey17" "grey18" [280] "grey19" "grey20" "grey21" [283] "grey22" "grey23" "grey24" [286] "grey25" "grey26" "grey27" [289] "grey28" "grey29" "grey30" [292] "grey31" "grey32" "grey33" [295] "grey34" "grey35" "grey36" [298] "grey37" "grey38" "grey39" [301] "grey40" "grey41" "grey42" [304] "grey43" "grey44" "grey45" [307] "grey46" "grey47" "grey48" [310] "grey49" "grey50" "grey51" [313] "grey52" "grey53" "grey54" [316] "grey55" "grey56" "grey57" [319] "grey58" "grey59" "grey60" [322] "grey61" "grey62" "grey63" [325] "grey64" "grey65" "grey66" [328] "grey67" "grey68" "grey69" [331] "grey70" "grey71" "grey72" [334] "grey73" "grey74" "grey75" [337] "grey76" "grey77" "grey78" [340] "grey79" "grey80" "grey81" [343] "grey82" "grey83" "grey84" [346] "grey85" "grey86" "grey87" [349] "grey88" "grey89" "grey90" [352] "grey91" "grey92" "grey93" [355] "grey94" "grey95" "grey96" [358] "grey97" "grey98" "grey99" [361] "grey100" "honeydew" "honeydew1" [364] "honeydew2" "honeydew3" "honeydew4" [367] "hotpink" "hotpink1" "hotpink2" [370] "hotpink3" "hotpink4" "indianred" [373] "indianred1" "indianred2" "indianred3" [376] "indianred4" "ivory" "ivory1" [379] "ivory2" "ivory3" "ivory4" [382] "khaki" "khaki1" "khaki2" [385] "khaki3" "khaki4" "lavender" [388] "lavenderblush" "lavenderblush1" "lavenderblush2" [391] "lavenderblush3" "lavenderblush4" "lawngreen" [394] "lemonchiffon" "lemonchiffon1" "lemonchiffon2" [397] "lemonchiffon3" "lemonchiffon4" "lightblue" [400] "lightblue1" "lightblue2" "lightblue3" [403] "lightblue4" "lightcoral" "lightcyan" [406] "lightcyan1" "lightcyan2" "lightcyan3" [409] "lightcyan4" "lightgoldenrod" "lightgoldenrod1" [412] "lightgoldenrod2" "lightgoldenrod3" "lightgoldenrod4" [415] "lightgoldenrodyellow" "lightgray" "lightgreen" [418] "lightgrey" "lightpink" "lightpink1" [421] "lightpink2" "lightpink3" "lightpink4" [424] "lightsalmon" "lightsalmon1" "lightsalmon2" [427] "lightsalmon3" "lightsalmon4" "lightseagreen" [430] "lightskyblue" "lightskyblue1" "lightskyblue2" [433] "lightskyblue3" "lightskyblue4" "lightslateblue" [436] "lightslategray" "lightslategrey" "lightsteelblue" [439] "lightsteelblue1" "lightsteelblue2" "lightsteelblue3" [442] "lightsteelblue4" "lightyellow" "lightyellow1" [445] "lightyellow2" "lightyellow3" "lightyellow4" [448] "limegreen" "linen" "magenta" [451] "magenta1" "magenta2" "magenta3" [454] "magenta4" "maroon" "maroon1" [457] "maroon2" "maroon3" "maroon4" [460] "mediumaquamarine" "mediumblue" "mediumorchid" [463] "mediumorchid1" "mediumorchid2" "mediumorchid3" [466] "mediumorchid4" "mediumpurple" "mediumpurple1" [469] "mediumpurple2" "mediumpurple3" "mediumpurple4" [472] "mediumseagreen" "mediumslateblue" "mediumspringgreen" [475] "mediumturquoise" "mediumvioletred" "midnightblue" [478] "mintcream" "mistyrose" "mistyrose1" [481] "mistyrose2" "mistyrose3" "mistyrose4" [484] "moccasin" "navajowhite" "navajowhite1" [487] "navajowhite2" "navajowhite3" "navajowhite4" [490] "navy" "navyblue" "oldlace" [493] "olivedrab" "olivedrab1" "olivedrab2" [496] "olivedrab3" "olivedrab4" "orange" [499] "orange1" "orange2" "orange3" [502] "orange4" "orangered" "orangered1" [505] "orangered2" "orangered3" "orangered4" [508] "orchid" "orchid1" "orchid2" [511] "orchid3" "orchid4" "palegoldenrod" [514] "palegreen" "palegreen1" "palegreen2" [517] "palegreen3" "palegreen4" "paleturquoise" [520] "paleturquoise1" "paleturquoise2" "paleturquoise3" [523] "paleturquoise4" "palevioletred" "palevioletred1" [526] "palevioletred2" "palevioletred3" "palevioletred4" [529] "papayawhip" "peachpuff" "peachpuff1" [532] "peachpuff2" "peachpuff3" "peachpuff4" [535] "peru" "pink" "pink1" [538] "pink2" "pink3" "pink4" [541] "plum" "plum1" "plum2" [544] "plum3" "plum4" "powderblue" [547] "purple" "purple1" "purple2" [550] "purple3" "purple4" "red" [553] "red1" "red2" "red3" [556] "red4" "rosybrown" "rosybrown1" [559] "rosybrown2" "rosybrown3" "rosybrown4" [562] "royalblue" "royalblue1" "royalblue2" [565] "royalblue3" "royalblue4" "saddlebrown" [568] "salmon" "salmon1" "salmon2" [571] "salmon3" "salmon4" "sandybrown" [574] "seagreen" "seagreen1" "seagreen2" [577] "seagreen3" "seagreen4" "seashell" [580] "seashell1" "seashell2" "seashell3" [583] "seashell4" "sienna" "sienna1" [586] "sienna2" "sienna3" "sienna4" [589] "skyblue" "skyblue1" "skyblue2" [592] "skyblue3" "skyblue4" "slateblue" [595] "slateblue1" "slateblue2" "slateblue3" [598] "slateblue4" "slategray" "slategray1" [601] "slategray2" "slategray3" "slategray4" [604] "slategrey" "snow" "snow1" [607] "snow2" "snow3" "snow4" [610] "springgreen" "springgreen1" "springgreen2" [613] "springgreen3" "springgreen4" "steelblue" [616] "steelblue1" "steelblue2" "steelblue3" [619] "steelblue4" "tan" "tan1" [622] "tan2" "tan3" "tan4" [625] "thistle" "thistle1" "thistle2" [628] "thistle3" "thistle4" "tomato" [631] "tomato1" "tomato2" "tomato3" [634] "tomato4" "turquoise" "turquoise1" [637] "turquoise2" "turquoise3" "turquoise4" [640] "violet" "violetred" "violetred1" [643] "violetred2" "violetred3" "violetred4" [646] "wheat" "wheat1" "wheat2" [649] "wheat3" "wheat4" "whitesmoke" [652] "yellow" "yellow1" "yellow2" [655] "yellow3" "yellow4" "yellowgreen" Nos gráficos seguintes vou utilizar “orangered3” (o mesmo que “#cd3700”) para o sexo Feminino e “turquoise4” (“#00868b”) para Masculino.

car::scatterplot(MCT~Estatura,
                 group=Dados$Sexo, 
                 jitter=list(x=1, y=1),
                 regLine=FALSE, 
                 smooth=FALSE, 
                 boxplots=TRUE, 
                 ellipse=list(levels=c(0.68), 
                              robust=TRUE, 
                              fill=FALSE, 
                              fill.alpha=0.2),
                 grid=FALSE,
                 col=c("orangered3","turquoise4"), 
                 xlim=c(145,200),
                 ylim=c(30,127),
                 data=Dados)

car::scatterplot(MCT~Estatura,
                 group=Dados$Sexo, 
                 jitter=list(x=1, y=1),
                 regLine=FALSE, 
                 smooth=FALSE, 
                 ellipse=list(levels=c(0.68,0.999), 
                              robust=TRUE, 
                              fill=FALSE, 
                              fill.alpha=0.2),
                 grid=FALSE,
                 col=c("orangered3","turquoise4"), 
                 xlim=c(145,200),
                 ylim=c(30,127),
                 data=Dados)

car::scatterplot(MCT~Estatura,
                 group=Dados$Sexo, 
                 jitter=list(x=1, y=1),
                 regLine=TRUE, 
                 smooth=FALSE, 
                 ellipse=FALSE,
                 grid=FALSE,
                 col=c("orangered3","turquoise4"), 
                 xlim=c(145,200),
                 ylim=c(30,127),
                 data=Dados)

Gráfico de dispersão matricial

car::scatterplotMatrix(Dados[,c("Estatura",
                                "MCT")], 
                       groups=Dados$Sexo,
                       lower.panel=NULL,
                       regLine=TRUE, 
                       smooth=FALSE, 
                       by.groups=TRUE,
                       ellipse=list(levels=c(0.5), 
                                    robust=TRUE, 
                                    fill=FALSE),
                       col=c("orangered3","turquoise4"), 
                       cex=0.5,
                       cex.labels=1,
                       row1attop=TRUE)

oldop <- options()
options(warn=-1)
  GGally::ggpairs(subset(Dados, 
                       select=-c(ID,Ano,Turma,Mao,TipoSang,
                                 ABO,AtivFisica,IMC)), 
                ggplot2::aes(colour=Sexo))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

options(oldop) 

Intervalo de confiança: comparação de médias populacionais

alpha <- 0.05
gplots::plotmeans(Dados$MCT~Dados$Sexo,
                  main="Intervalo de confiança de 95%\ncom correção de Bonferroni",
                  col="black",
                  barcol="black",
                  p=1-alpha/length(unique(levels(Dados$Sexo))),
                  connect=FALSE)

gplots::plotmeans(Dados$MCT~Dados$Sedentarismo,
                  main="Intervalo de confiança de 95%\ncom correção de Bonferroni",
                  col="black",
                  barcol="black",
                  p=1-alpha/length(unique(levels(Dados$Sexo))),
                  connect=FALSE)

gplots::plotmeans(Dados$MCT~ interaction(Dados$Sexo, 
                                         Dados$Sedentarismo),
                  main="Intervalo de confiança de 95%\ncom correção de Bonferroni",
                  col="black",
                  barcol="black",
                  p=1-alpha/(length(unique(levels(Dados$Sexo)))*
                    length(unique(levels(Dados$Sedentarismo)))),
                  connect=FALSE)

Gráficos avançados

• ggplot2
  • https://ggplot2.tidyverse.org/
• Módulo V: Gráficos com ggplot2 por Bruno Lucian & Claudia Jakelline
  • https://bookdown.org/bruno_lucian_costa/CursoBasicoR/m5.html#gr%C3%A1ficos-com-ggplot2
• ggplot2 - Essentials
  • http://www.sthda.com/english/wiki/ggplot2-essentials
• The R Graph Gallery
  • https://r-graph-gallery.com/index.html
• Advanced plotting by Jeff Powell and Remko Duursma
  • http://www.hiercourse.com/advancedplot
  • http://www.hiercourse.com/docs/advanced_plotting.html
• Chapter 19. Advanced graphics with ggplot2
  • https://livebook.manning.com/book/r-in-action-second-edition/chapter-19/1
• Top 50 ggplot2 Visualizations - The Master List (With Full R Code)
  • http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html

Ler arquivo de dados com medidas repetidas

Ler e transformar a estrutura de dados de wide para long

• NewDrug.xls
  • Fonte dos dados: Norusis, M (2021) IBM SPSS Statistics 19 Statistical Procedures Companion. Chapter: GLM Repeated Measures. USA: Addison Wesley.

Este arquivo tem 12 pacientes, em dois grupos (nova droga e placebo), dada um deles com três medidas de resp e três de pulse.

Dados.wide <- readxl::read_xlsx("NewDrug.xlsx")
Dados.wide$drug <- factor(Dados.wide$drug)
tabela <- knitr::kable(Dados.wide, format="html", caption="New Drug: wide")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: wide

drug	resp1	resp2	resp3	pulse1	pulse2	pulse3
New Drug	3.4	3.3	3.3	2.2	2.1	2.1
New Drug	3.4	3.4	3.3	2.2	2.1	2.2
New Drug	3.3	3.4	3.4	2.3	2.4	2.3
New Drug	3.4	3.4	3.4	2.3	2.4	2.3
New Drug	3.3	3.4	3.3	2.2	2.2	2.4
New Drug	3.3	3.3	3.3	2.0	2.1	2.4
Placebo	3.3	3.3	3.3	2.8	2.9	2.7
Placebo	3.2	3.3	3.4	2.6	2.7	2.7
Placebo	3.2	3.2	3.2	2.7	2.9	2.7
Placebo	3.2	3.2	3.2	2.6	2.8	2.9
Placebo	3.2	3.3	3.3	2.7	2.8	2.9
Placebo	3.3	3.2	3.1	2.6	2.8	2.8

str(Dados.wide)

tibble [12 × 7] (S3: tbl_df/tbl/data.frame)
 $ drug  : Factor w/ 2 levels "New Drug","Placebo": 1 1 1 1 1 1 2 2 2 2 ...
 $ resp1 : num [1:12] 3.4 3.4 3.3 3.4 3.3 3.3 3.3 3.2 3.2 3.2 ...
 $ resp2 : num [1:12] 3.3 3.4 3.4 3.4 3.4 3.3 3.3 3.3 3.2 3.2 ...
 $ resp3 : num [1:12] 3.3 3.3 3.4 3.4 3.3 3.3 3.3 3.4 3.2 3.2 ...
 $ pulse1: num [1:12] 2.2 2.2 2.3 2.3 2.2 2 2.8 2.6 2.7 2.6 ...
 $ pulse2: num [1:12] 2.1 2.1 2.4 2.4 2.2 2.1 2.9 2.7 2.9 2.8 ...
 $ pulse3: num [1:12] 2.1 2.2 2.3 2.3 2.4 2.4 2.7 2.7 2.7 2.9 ...

descricao <- psych::describeBy(subset(Dados.wide,select=c(-drug)),
                               group=list(Dados.wide$drug),
                               mat=TRUE,
                               digits=2)
tabela <- knitr::kable(descricao, format="html", caption="New Drug: wide")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: wide

	item	group1	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
resp11	1	New Drug	1	6	3.35	0.05	3.35	3.35	0.07	3.3	3.4	0.1	0.00	-2.31	0.02
resp12	2	Placebo	1	6	3.23	0.05	3.20	3.23	0.00	3.2	3.3	0.1	0.54	-1.96	0.02
resp21	3	New Drug	2	6	3.37	0.05	3.40	3.37	0.00	3.3	3.4	0.1	-0.54	-1.96	0.02
resp22	4	Placebo	2	6	3.25	0.05	3.25	3.25	0.07	3.2	3.3	0.1	0.00	-2.31	0.02
resp31	5	New Drug	3	6	3.33	0.05	3.30	3.33	0.00	3.3	3.4	0.1	0.54	-1.96	0.02
resp32	6	Placebo	3	6	3.25	0.10	3.25	3.25	0.07	3.1	3.4	0.3	0.00	-1.57	0.04
pulse11	7	New Drug	4	6	2.20	0.11	2.20	2.20	0.07	2.0	2.3	0.3	-0.76	-0.92	0.04
pulse12	8	Placebo	4	6	2.67	0.08	2.65	2.67	0.07	2.6	2.8	0.2	0.48	-1.58	0.03
pulse21	9	New Drug	5	6	2.22	0.15	2.15	2.22	0.07	2.1	2.4	0.3	0.39	-2.00	0.06
pulse22	10	Placebo	5	6	2.82	0.08	2.80	2.82	0.07	2.7	2.9	0.2	-0.17	-1.54	0.03
pulse31	11	New Drug	6	6	2.28	0.12	2.30	2.28	0.15	2.1	2.4	0.3	-0.37	-1.62	0.05
pulse32	12	Placebo	6	6	2.78	0.10	2.75	2.78	0.07	2.7	2.9	0.2	0.25	-2.08	0.04

Dados.wide <- cbind(seq(1:nrow(Dados.wide)),Dados.wide)
names(Dados.wide) <- c("ID",names(Dados.wide[,2:ncol(Dados.wide)]))
Dados.wide$ID <- factor(Dados.wide$ID)
tabela <- knitr::kable(Dados.wide, format="html", caption="New Drug: wide")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: wide

ID	drug	resp1	resp2	resp3	pulse1	pulse2	pulse3
1	New Drug	3.4	3.3	3.3	2.2	2.1	2.1
2	New Drug	3.4	3.4	3.3	2.2	2.1	2.2
3	New Drug	3.3	3.4	3.4	2.3	2.4	2.3
4	New Drug	3.4	3.4	3.4	2.3	2.4	2.3
5	New Drug	3.3	3.4	3.3	2.2	2.2	2.4
6	New Drug	3.3	3.3	3.3	2.0	2.1	2.4
7	Placebo	3.3	3.3	3.3	2.8	2.9	2.7
8	Placebo	3.2	3.3	3.4	2.6	2.7	2.7
9	Placebo	3.2	3.2	3.2	2.7	2.9	2.7
10	Placebo	3.2	3.2	3.2	2.6	2.8	2.9
11	Placebo	3.2	3.3	3.3	2.7	2.8	2.9
12	Placebo	3.3	3.2	3.1	2.6	2.8	2.8

Dados.long <- Hmisc::reShape(Dados.wide,
                             id=c("ID"),
                             colvar=c("drug"),
                             base=c("resp", "pulse"),
                             reps=3,
                             timevar="time",
                             times=c(1,2,3))
str(Dados.long)

'data.frame':   36 obs. of  5 variables:
 $ time : num  1 2 3 1 2 3 1 2 3 1 ...
 $ ID   : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
 $ drug : Factor w/ 2 levels "New Drug","Placebo": 1 1 1 1 1 1 1 1 1 1 ...
 $ resp : num  3.4 3.3 3.3 3.4 3.4 3.3 3.3 3.4 3.4 3.4 ...
 $ pulse: num  2.2 2.1 2.1 2.2 2.1 2.2 2.3 2.4 2.3 2.3 ...

tabela <- knitr::kable(Dados.long, format="html", caption="New Drug: long")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: long

	time	ID	drug	resp	pulse
1 1	1	1	New Drug	3.4	2.2
1 2	2	1	New Drug	3.3	2.1
1 3	3	1	New Drug	3.3	2.1
2 1	1	2	New Drug	3.4	2.2
2 2	2	2	New Drug	3.4	2.1
2 3	3	2	New Drug	3.3	2.2
3 1	1	3	New Drug	3.3	2.3
3 2	2	3	New Drug	3.4	2.4
3 3	3	3	New Drug	3.4	2.3
4 1	1	4	New Drug	3.4	2.3
4 2	2	4	New Drug	3.4	2.4
4 3	3	4	New Drug	3.4	2.3
5 1	1	5	New Drug	3.3	2.2
5 2	2	5	New Drug	3.4	2.2
5 3	3	5	New Drug	3.3	2.4
6 1	1	6	New Drug	3.3	2.0
6 2	2	6	New Drug	3.3	2.1
6 3	3	6	New Drug	3.3	2.4
7 1	1	7	Placebo	3.3	2.8
7 2	2	7	Placebo	3.3	2.9
7 3	3	7	Placebo	3.3	2.7
8 1	1	8	Placebo	3.2	2.6
8 2	2	8	Placebo	3.3	2.7
8 3	3	8	Placebo	3.4	2.7
9 1	1	9	Placebo	3.2	2.7
9 2	2	9	Placebo	3.2	2.9
9 3	3	9	Placebo	3.2	2.7
10 1	1	10	Placebo	3.2	2.6
10 2	2	10	Placebo	3.2	2.8
10 3	3	10	Placebo	3.2	2.9
11 1	1	11	Placebo	3.2	2.7
11 2	2	11	Placebo	3.3	2.8
11 3	3	11	Placebo	3.3	2.9
12 1	1	12	Placebo	3.3	2.6
12 2	2	12	Placebo	3.2	2.8
12 3	3	12	Placebo	3.1	2.8

Temos dois data frames, nos formatos wide e long. Podemos juntá-los em um único objeto, neste caso uma lista, e salvá-los para uso futuro.

Dados.lista <- list(Dados.wide, Dados.long)
names(Dados.lista) <- c("wide", "long")
saveRDS(Dados.lista, "NewDrug.rds")

Uma vez salvo, você pode recuperar o arquivo quando precisar:

Dados <- readRDS("NewDrug.rds")
str(Dados)

List of 2
 $ wide:'data.frame':   12 obs. of  8 variables:
  ..$ ID    : Factor w/ 12 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
  ..$ drug  : Factor w/ 2 levels "New Drug","Placebo": 1 1 1 1 1 1 2 2 2 2 ...
  ..$ resp1 : num [1:12] 3.4 3.4 3.3 3.4 3.3 3.3 3.3 3.2 3.2 3.2 ...
  ..$ resp2 : num [1:12] 3.3 3.4 3.4 3.4 3.4 3.3 3.3 3.3 3.2 3.2 ...
  ..$ resp3 : num [1:12] 3.3 3.3 3.4 3.4 3.3 3.3 3.3 3.4 3.2 3.2 ...
  ..$ pulse1: num [1:12] 2.2 2.2 2.3 2.3 2.2 2 2.8 2.6 2.7 2.6 ...
  ..$ pulse2: num [1:12] 2.1 2.1 2.4 2.4 2.2 2.1 2.9 2.7 2.9 2.8 ...
  ..$ pulse3: num [1:12] 2.1 2.2 2.3 2.3 2.4 2.4 2.7 2.7 2.7 2.9 ...
 $ long:'data.frame':   36 obs. of  5 variables:
  ..$ time : num [1:36] 1 2 3 1 2 3 1 2 3 1 ...
  ..$ ID   : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
  ..$ drug : Factor w/ 2 levels "New Drug","Placebo": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ resp : num [1:36] 3.4 3.3 3.3 3.4 3.4 3.3 3.3 3.4 3.4 3.4 ...
  ..$ pulse: num [1:36] 2.2 2.1 2.1 2.2 2.1 2.2 2.3 2.4 2.3 2.3 ...

tabela <- knitr::kable(Dados$wide, format="html", caption="New Drug: wide")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: wide

ID	drug	resp1	resp2	resp3	pulse1	pulse2	pulse3
1	New Drug	3.4	3.3	3.3	2.2	2.1	2.1
2	New Drug	3.4	3.4	3.3	2.2	2.1	2.2
3	New Drug	3.3	3.4	3.4	2.3	2.4	2.3
4	New Drug	3.4	3.4	3.4	2.3	2.4	2.3
5	New Drug	3.3	3.4	3.3	2.2	2.2	2.4
6	New Drug	3.3	3.3	3.3	2.0	2.1	2.4
7	Placebo	3.3	3.3	3.3	2.8	2.9	2.7
8	Placebo	3.2	3.3	3.4	2.6	2.7	2.7
9	Placebo	3.2	3.2	3.2	2.7	2.9	2.7
10	Placebo	3.2	3.2	3.2	2.6	2.8	2.9
11	Placebo	3.2	3.3	3.3	2.7	2.8	2.9
12	Placebo	3.3	3.2	3.1	2.6	2.8	2.8

tabela <- knitr::kable(Dados$long, format="html", caption="New Drug: long")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: long

	time	ID	drug	resp	pulse
1 1	1	1	New Drug	3.4	2.2
1 2	2	1	New Drug	3.3	2.1
1 3	3	1	New Drug	3.3	2.1
2 1	1	2	New Drug	3.4	2.2
2 2	2	2	New Drug	3.4	2.1
2 3	3	2	New Drug	3.3	2.2
3 1	1	3	New Drug	3.3	2.3
3 2	2	3	New Drug	3.4	2.4
3 3	3	3	New Drug	3.4	2.3
4 1	1	4	New Drug	3.4	2.3
4 2	2	4	New Drug	3.4	2.4
4 3	3	4	New Drug	3.4	2.3
5 1	1	5	New Drug	3.3	2.2
5 2	2	5	New Drug	3.4	2.2
5 3	3	5	New Drug	3.3	2.4
6 1	1	6	New Drug	3.3	2.0
6 2	2	6	New Drug	3.3	2.1
6 3	3	6	New Drug	3.3	2.4
7 1	1	7	Placebo	3.3	2.8
7 2	2	7	Placebo	3.3	2.9
7 3	3	7	Placebo	3.3	2.7
8 1	1	8	Placebo	3.2	2.6
8 2	2	8	Placebo	3.3	2.7
8 3	3	8	Placebo	3.4	2.7
9 1	1	9	Placebo	3.2	2.7
9 2	2	9	Placebo	3.2	2.9
9 3	3	9	Placebo	3.2	2.7
10 1	1	10	Placebo	3.2	2.6
10 2	2	10	Placebo	3.2	2.8
10 3	3	10	Placebo	3.2	2.9
11 1	1	11	Placebo	3.2	2.7
11 2	2	11	Placebo	3.3	2.8
11 3	3	11	Placebo	3.3	2.9
12 1	1	12	Placebo	3.3	2.6
12 2	2	12	Placebo	3.2	2.8
12 3	3	12	Placebo	3.1	2.8

Ler e transformar a estrutura de dados de long para wide

Dados.wide.resp <- reshape2::dcast(Dados.long, 
                                   ID + drug ~ time, 
                                   value.var="resp")
names(Dados.wide.resp) <- c("ID","drug","resp1","resp2","resp3")
tabela <- knitr::kable(Dados.wide.resp, format="html", caption="New Drug: resp wide")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: resp wide

ID	drug	resp1	resp2	resp3
1	New Drug	3.4	3.3	3.3
2	New Drug	3.4	3.4	3.3
3	New Drug	3.3	3.4	3.4
4	New Drug	3.4	3.4	3.4
5	New Drug	3.3	3.4	3.3
6	New Drug	3.3	3.3	3.3
7	Placebo	3.3	3.3	3.3
8	Placebo	3.2	3.3	3.4
9	Placebo	3.2	3.2	3.2
10	Placebo	3.2	3.2	3.2
11	Placebo	3.2	3.3	3.3
12	Placebo	3.3	3.2	3.1

Dados.wide.pulse <- reshape2::dcast(Dados.long, 
                                    ID + drug ~ time, 
                                    value.var="pulse")
names(Dados.wide.pulse) <- c("ID","drug","pulse1","pulse2","pulse3")
tabela <- knitr::kable(Dados.wide.pulse, format="html", caption="New Drug: pulse wide")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: pulse wide

ID	drug	pulse1	pulse2	pulse3
1	New Drug	2.2	2.1	2.1
2	New Drug	2.2	2.1	2.2
3	New Drug	2.3	2.4	2.3
4	New Drug	2.3	2.4	2.3
5	New Drug	2.2	2.2	2.4
6	New Drug	2.0	2.1	2.4
7	Placebo	2.8	2.9	2.7
8	Placebo	2.6	2.7	2.7
9	Placebo	2.7	2.9	2.7
10	Placebo	2.6	2.8	2.9
11	Placebo	2.7	2.8	2.9
12	Placebo	2.6	2.8	2.8

Dados.wide <- cbind(Dados.wide.resp, 
                    Dados.wide.pulse[c("pulse1","pulse2","pulse3")])
tabela <- knitr::kable(Dados.wide, format="html", caption="New Drug: wide")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

New Drug: wide

ID	drug	resp1	resp2	resp3	pulse1	pulse2	pulse3
1	New Drug	3.4	3.3	3.3	2.2	2.1	2.1
2	New Drug	3.4	3.4	3.3	2.2	2.1	2.2
3	New Drug	3.3	3.4	3.4	2.3	2.4	2.3
4	New Drug	3.4	3.4	3.4	2.3	2.4	2.3
5	New Drug	3.3	3.4	3.3	2.2	2.2	2.4
6	New Drug	3.3	3.3	3.3	2.0	2.1	2.4
7	Placebo	3.3	3.3	3.3	2.8	2.9	2.7
8	Placebo	3.2	3.3	3.4	2.6	2.7	2.7
9	Placebo	3.2	3.2	3.2	2.7	2.9	2.7
10	Placebo	3.2	3.2	3.2	2.6	2.8	2.9
11	Placebo	3.2	3.3	3.3	2.7	2.8	2.9
12	Placebo	3.3	3.2	3.1	2.6	2.8	2.8

Ler dados dentro de script R

variáveis em colunas e casos em linhas

# https://rcompanion.org 
 
Dados <- read.table(header=TRUE, 
                    stringsAsFactors=TRUE, 
                    text="
 Medico   Conceito
 A        3
 A        5
 A        4
 A        4
 A        4
 A        4
 A        4
 A        4
 A        5
 A        5
 B        2
 B        4
 B        2
 B        2
 B        1
 B        2
 B        3
 B        2
 B        2
 B        3
")
Dados$Conceito <- factor(Dados$Conceito, ordered=TRUE)
levels(Dados$Conceito)

[1] "1" "2" "3" "4" "5"

class(Dados)

[1] "data.frame"

str(Dados)

'data.frame':   20 obs. of  2 variables:
 $ Medico  : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ Conceito: Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 3 5 4 4 4 4 4 4 5 5 ...

tabela <- knitr::kable(Dados, format="html", caption="Médico: conceito")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

Médico: conceito

Medico	Conceito
A	3
A	5
A	4
A	4
A	4
A	4
A	4
A	4
A	5
A	5
B	2
B	4
B	2
B	2
B	1
B	2
B	3
B	2
B	2
B	3

summary(Dados)

 Medico Conceito
 A:10   1:1     
 B:10   2:6     
        3:3     
        4:7     
        5:3     

FSA::Summarize(as.numeric(Conceito) ~ Medico,
               digits=3,
               data=Dados)

# A tibble: 2 × 9
  Medico     n  mean    sd   min    Q1 median    Q3   max
  <fct>  <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
1 A         10   4.2 0.632     3     4      4  4.75     5
2 B         10   2.3 0.823     1     2      2  2.75     4

lattice::histogram(~ Conceito | Medico,
                   data=Dados,
                   aspect="fill",
                   col="lightgray",
                   layout=c(1,2)) 

tabela <- xtabs(~ Medico + Conceito, data=Dados)
addmargins(tabela)

      Conceito
Medico  1  2  3  4  5 Sum
   A    0  0  1  6  3  10
   B    1  6  2  1  0  10
   Sum  1  6  3  7  3  20

round(prop.table(tabela, margin=1), 2)

      Conceito
Medico   1   2   3   4   5
     A 0.0 0.0 0.1 0.6 0.3
     B 0.1 0.6 0.2 0.1 0.0

spineplot(tabela)

mosaicplot(tabela, shade=TRUE)

saveRDS(object=Dados, file="MediConceito.rds")
rm(Dados)
Dados <- readRDS(file="MediConceito.rds")
class(Dados)

[1] "data.frame"

str(Dados)

'data.frame':   20 obs. of  2 variables:
 $ Medico  : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ Conceito: Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 3 5 4 4 4 4 4 4 5 5 ...

tabela bidimensional

Para este exemplo hipotético, os agricultores foram pesquisados sobre com que frequência eles utilizam uma determinada prática de manejo adequado. As respostas são organizadas de acordo com o tamanho da operação. Ambas as variáveis na tabela de contingência são categorias ordenadas.

tabela <- ("
                     Sempre 'Às Vezes' Nunca
   Hobbista          0       1         5
   'Pequeno negócio' 2       3         4
   Pequeno           4       4         4
   Médio             3       2         0
   Grande            2       0         0
")

Dados.tabela <- as.matrix(read.table(textConnection(tabela), 
                                     header=TRUE, 
                                     row.names=1))
names(dimnames(Dados.tabela)) <- c("Tamanho", "Frequência")
Dados.tabela <- as.table(Dados.tabela)
addmargins(Dados.tabela)

                 Frequência
Tamanho           Sempre Às.Vezes Nunca Sum
  Hobbista             0        1     5   6
  Pequeno negócio      2        3     4   9
  Pequeno              4        4     4  12
  Médio                3        2     0   5
  Grande               2        0     0   2
  Sum                 11       10    13  34

str(Dados.tabela)

 'table' int [1:5, 1:3] 0 2 4 3 2 1 3 4 2 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ Tamanho   : chr [1:5] "Hobbista" "Pequeno negócio" "Pequeno" "Médio" ...
  ..$ Frequência: chr [1:3] "Sempre" "Às.Vezes" "Nunca"

gplots::balloonplot(t(Dados.tabela), 
                    main ="Tabela de contingência\nDuas variáveis ordinais", 
                    xlab ="Freq prática manejo", 
                    ylab="Tamanho operação",
                    dotcolor="gray45",
                    show.zeros=TRUE,
                    show.margins=TRUE)

Dados <- DescTools::Untable(Dados.tabela)
print.data.frame(Dados)

           Tamanho Frequência
1  Pequeno negócio     Sempre
2  Pequeno negócio     Sempre
3          Pequeno     Sempre
4          Pequeno     Sempre
5          Pequeno     Sempre
6          Pequeno     Sempre
7            Médio     Sempre
8            Médio     Sempre
9            Médio     Sempre
10          Grande     Sempre
11          Grande     Sempre
12        Hobbista   Às.Vezes
13 Pequeno negócio   Às.Vezes
14 Pequeno negócio   Às.Vezes
15 Pequeno negócio   Às.Vezes
16         Pequeno   Às.Vezes
17         Pequeno   Às.Vezes
18         Pequeno   Às.Vezes
19         Pequeno   Às.Vezes
20           Médio   Às.Vezes
21           Médio   Às.Vezes
22        Hobbista      Nunca
23        Hobbista      Nunca
24        Hobbista      Nunca
25        Hobbista      Nunca
26        Hobbista      Nunca
27 Pequeno negócio      Nunca
28 Pequeno negócio      Nunca
29 Pequeno negócio      Nunca
30 Pequeno negócio      Nunca
31         Pequeno      Nunca
32         Pequeno      Nunca
33         Pequeno      Nunca
34         Pequeno      Nunca

Dados$Tamanho <- factor(Dados$Tamanho, 
                        ordered=TRUE,
                        levels=rownames(Dados.tabela))
Dados$Frequência <- factor(Dados$Frequência, 
                           ordered=TRUE,
                           levels=colnames(Dados.tabela))

str(Dados)

'data.frame':   34 obs. of  2 variables:
 $ Tamanho   : Ord.factor w/ 5 levels "Hobbista"<"Pequeno negócio"<..: 2 2 3 3 3 3 4 4 4 5 ...
 $ Frequência: Ord.factor w/ 3 levels "Sempre"<"Às.Vezes"<..: 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "out.attrs")=List of 2
  ..$ dim     : Named int [1:2] 5 3
  .. ..- attr(*, "names")= chr [1:2] "Tamanho" "Frequência"
  ..$ dimnames:List of 2
  .. ..$ Tamanho   : chr [1:5] "Tamanho=Hobbista" "Tamanho=Pequeno negócio" "Tamanho=Pequeno" "Tamanho=Médio" ...
  .. ..$ Frequência: chr [1:3] "Frequência=Sempre" "Frequência=Às.Vezes" "Frequência=Nunca"

summary(Dados)

            Tamanho      Frequência
 Hobbista       : 6   Sempre  :11  
 Pequeno negócio: 9   Às.Vezes:10  
 Pequeno        :12   Nunca   :13  
 Médio          : 5                
 Grande         : 2                

xtabs(~Tamanho+Frequência, data=Dados)

                 Frequência
Tamanho           Sempre Às.Vezes Nunca
  Hobbista             0        1     5
  Pequeno negócio      2        3     4
  Pequeno              4        4     4
  Médio                3        2     0
  Grande               2        0     0

round(prop.table(Dados.tabela,
           margin=NULL),2)

                 Frequência
Tamanho           Sempre Às.Vezes Nunca
  Hobbista          0.00     0.03  0.15
  Pequeno negócio   0.06     0.09  0.12
  Pequeno           0.12     0.12  0.12
  Médio             0.09     0.06  0.00
  Grande            0.06     0.00  0.00

spineplot(Dados.tabela)

mosaicplot(Dados.tabela, shade=TRUE)

Da mesma forma que fizemos no exemplo acima, podemos salvar os dados, remover da memória e recuperar a lista para mostrar que o trabalho fica preservado:

DadosLista <- list(Dados.tabela, Dados)
names(DadosLista) <- c("Dados.tabela", "Dados")
saveRDS(object=DadosLista, file="PraticaManejoLista.rds")
rm(DadosLista)
DadosLista <- readRDS(file="PraticaManejoLista.rds")
DadosLista$Dados.tabela

                 Frequência
Tamanho           Sempre Às.Vezes Nunca
  Hobbista             0        1     5
  Pequeno negócio      2        3     4
  Pequeno              4        4     4
  Médio                3        2     0
  Grande               2        0     0

DadosLista$Dados

# A tibble: 34 × 2
   Tamanho         Frequência
   <ord>           <ord>     
 1 Pequeno negócio Sempre    
 2 Pequeno negócio Sempre    
 3 Pequeno         Sempre    
 4 Pequeno         Sempre    
 5 Pequeno         Sempre    
 6 Pequeno         Sempre    
 7 Médio           Sempre    
 8 Médio           Sempre    
 9 Médio           Sempre    
10 Grande          Sempre    
# ℹ 24 more rows

tabela tridimensional

Tabela de contingência 2x2 estratificada com relação entre tabagismo e sobrevivência em 20 anos (1974-1994) em 1.134 mulheres adultas do Reino Unido.

• Delineamento: coorte
• Fonte: APPLETON, DR et al. (1996) Ignoring a covariate: An example of Simpson’s paradox. The American Statistician, 50(4): 340-1.

# A aspa inicial TEM que começar na primeira coluna da linha
# e os espaçamentos distintos dos dessa tabela podem causar 
# problemas na geração da tabela horizontalizada (ftable).
tabela3D <- (
"         FaixaEtaria 18-24 25-34 35-44 45-54 55-64 65-74 75+
Tabagista Desfecho 
Sim       Morta       2     3     14    27    51    29   13  
          Viva        53    121   95    103   64    7    0 
Nao       Morta       1     5     7     12    40    101  64
          Viva        61    152   114   66    81    28   0
")
Dados.tabela3D <- as.table(read.ftable(textConnection(tabela3D)))
ftable(Dados.tabela3D) # Display a flattened table (tabela horizontalizada)

                   FaixaEtaria 18-24 25-34 35-44 45-54 55-64 65-74 75+
Tabagista Desfecho                                                    
Sim       Morta                    2     3    14    27    51    29  13
          Viva                    53   121    95   103    64     7   0
Nao       Morta                    1     5     7    12    40   101  64
          Viva                    61   152   114    66    81    28   0

mosaicplot(Dados.tabela3D, shade=TRUE)

saveRDS(object=Dados.tabela3D, file="Tabagismo3D.rds")
rm(Dados.tabela3D)
Dados.tabela3D <- readRDS(file="Tabagismo3D.rds")
Dados.tabela2D <- margin.table(Dados.tabela3D, margin=c(1,2))
addmargins(Dados.tabela2D)

         Desfecho
Tabagista Morta Viva  Sum
      Sim   139  443  582
      Nao   230  502  732
      Sum   369  945 1314

Dados.tabela.margin <- addmargins(Dados.tabela3D, margin=c(1,2))
Dados.tabela.margin

, , FaixaEtaria = 18-24

         Desfecho
Tabagista Morta Viva Sum
      Sim     2   53  55
      Nao     1   61  62
      Sum     3  114 117

, , FaixaEtaria = 25-34

         Desfecho
Tabagista Morta Viva Sum
      Sim     3  121 124
      Nao     5  152 157
      Sum     8  273 281

, , FaixaEtaria = 35-44

         Desfecho
Tabagista Morta Viva Sum
      Sim    14   95 109
      Nao     7  114 121
      Sum    21  209 230

, , FaixaEtaria = 45-54

         Desfecho
Tabagista Morta Viva Sum
      Sim    27  103 130
      Nao    12   66  78
      Sum    39  169 208

, , FaixaEtaria = 55-64

         Desfecho
Tabagista Morta Viva Sum
      Sim    51   64 115
      Nao    40   81 121
      Sum    91  145 236

, , FaixaEtaria = 65-74

         Desfecho
Tabagista Morta Viva Sum
      Sim    29    7  36
      Nao   101   28 129
      Sum   130   35 165

, , FaixaEtaria = 75+

         Desfecho
Tabagista Morta Viva Sum
      Sim    13    0  13
      Nao    64    0  64
      Sum    77    0  77

DadosLista <- list(Dados.tabela3D, Dados.tabela2D)
names(DadosLista) <- c("Dados.tabela3D", "Dados.tabela2D")
saveRDS(object=DadosLista, file="DadosLista.rds")
rm(DadosLista)
tabelas <- readRDS(file="DadosLista.rds")
ftable(tabelas$Dados.tabela3D)

                   FaixaEtaria 18-24 25-34 35-44 45-54 55-64 65-74 75+
Tabagista Desfecho                                                    
Sim       Morta                    2     3    14    27    51    29  13
          Viva                    53   121    95   103    64     7   0
Nao       Morta                    1     5     7    12    40   101  64
          Viva                    61   152   114    66    81    28   0

ftable(tabelas$Dados.tabela2D)

          Desfecho Morta Viva
Tabagista                    
Sim                  139  443
Nao                  230  502

Ler arquivo de dados grande

• DOBR2020.dbc: Declarações de Óbitos (DO) de BR em 2020 do Sistema de Informações de Mortalidade (SIM) do DATASUS: https://datasus.saude.gov.br/transferencia-de-arquivos/

A documentação sobre as variáveis está no arquivo Estrutura_SIM_para_CD.pdf.

As definições das variáveis e de seus níveis estão em Convenções SIM.

Anteriormente, era possível baixar diretamente do DATASUS o arquivo com dados do Brasil sem separação por unidade da federação (UF). Atualmente, é possível baixar os dados por UF.

O arquivo DOBR2020.dbc usado aqui não tem separação por UF.

Outra maneira de ler arquivo de microdados do DATASUS é a seguinte:

# https://medium.com/@danielly.bx/realizando-download-de-dados-p%C3%BAblicos-do-datasus-com-o-pacote-microdatasus-do-r-218cf4181e47
# https://github.com/rfsaldanha/microdatasus

## fetch_datasus()
# Sistema de Informação sobre Mortalidade (SIM);
# Sistema de Informação sobre Nascidos Vivos (SINASC);
# Sistema de Informações Hospitalares (SIH, SIH-RD (AIH reduzida));
# Cadastro Nacional de Estabelecimentos de Saúde (CNES, CNES-PF (apenas arquivos de profissionais));
# Sistema de Informações Ambulatoriais (SIA);
# Sistema de Informação de Agravos de Notificação (SINAN), incluindo o SINAN-Dengue, SINAN-Chikingunya, SINAN-Zika e SINAN-Malária.

## fetch_sigtab()
# Sistema de Gerenciamento da Tabela de Procedimentos, Medicamentos e OPM do SUS (SIGTAP)

# library(remotes)
# remotes::install_github("rfsaldanha/microdatasus")
# library(microdatasus)

# SIM - Sistema de Informaçõess de Mortalidade
# DO - Declarações de Óbitos 
# 2020 - BR
# 1.556.824 linhas x 88 colunas

# info_sys <- c("SIM-DO")
# uf <- c("all") # all = todas as UF do Brasil
# # Transforma .dbc em dataframe
# DOBR2020 <- microdatasus::fetch_datasus(year_start = 2020, 
#                                           year_end = 2020, 
#                                           month_start = 12, 
#                                           month_end = 12, 
#                                           uf = uf, 
#                                           information_system = info_sys)
# str(DOBR2020)
# saveRDS(DOBR2020, "DOBR2020.rds")

Função para obter a idade de óbito

# precisaremos desta funcao adiante para obter a idade de óbito
age <- function(dob, age.death, units="years", floor=TRUE) {
                calc.age <- lubridate::interval(dob, age.death)/
                            lubridate::duration(num = 1, units=units)
                if (floor) return(as.integer(floor(calc.age)))
                return(calc.age) 
}

Dicionário de dados

Ler DOBR2020.dbc

# suppressMessages(library(devtools, warn.conflicts=FALSE))
# devtools::install_github("danicat/read.dbc")
# suppressMessages(library(read.dbc, warn.conflicts=FALSE))

Dados <- read.dbc::read.dbc("DOBR2020.dbc")
sapply(Dados,class)

    ORIGEM   TIPOBITO    DTOBITO  HORAOBITO    NATURAL CODMUNNATU     DTNASC 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
     IDADE       SEXO    RACACOR     ESTCIV        ESC    ESC2010 SERIESCFAL 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
      OCUP  CODMUNRES    LOCOCOR   CODESTAB ESTABDESCR CODMUNOCOR   IDADEMAE 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
    ESCMAE ESCMAE2010 SERIESCMAE    OCUPMAE QTDFILVIVO QTDFILMORT   GRAVIDEZ 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
SEMAGESTAC   GESTACAO      PARTO OBITOPARTO       PESO TPMORTEOCO  OBITOGRAV 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
OBITOPUERP  ASSISTMED      EXAME   CIRURGIA  NECROPSIA     LINHAA     LINHAB 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
    LINHAC     LINHAD    LINHAII   CAUSABAS     CB_PRE COMUNSVOIM DTATESTADO 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
 CIRCOBITO   ACIDTRAB      FONTE NUMEROLOTE      TPPOS DTINVESTIG CAUSABAS_O 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
DTCADASTRO  ATESTANTE STCODIFICA CODIFICADO VERSAOSIST  VERSAOSCB   FONTEINV 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
 DTRECEBIM   ATESTADO DTRECORIGA   CAUSAMAT ESCMAEAGR1 ESCFALAGR1 STDOEPIDEM 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
  STDONOVA    DIFDATA NUDIASOBCO NUDIASOBIN   DTCADINV TPOBITOCOR   DTCONINV 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
    FONTES TPRESGINFO TPNIVELINV  NUDIASINF   DTCADINF MORTEPARTO  DTCONCASO 
  "factor"   "factor"   "factor"   "factor"   "factor"   "factor"   "factor" 
 FONTESINF   ALTCAUSA   CONTADOR 
  "factor"   "factor"   "factor" 

Exemplo de uma inconsistência

Este arquivo tem diversas inconsistências com sua documentação. Por exemplo, para IDADE dizem:

  Idade: composto de dois subcampos.
  - O primeiro, de 1 dígito, indica a unidade da idade (se 1 = minuto, se 2 = hora, se 3 = mês, se 4 = ano, se = 5 idade maior que 100 anos).
  - O segundo, de dois dígitos, indica a quantidade de unidades: 
    Idade menor de 1 hora: subcampo varia de 01 e 59 (minutos); 
    De 1 a 23 Horas: subcampo varia de 01 a 23 (horas); 
    De 24 horas e 29 dias: subcampo varia de 01 a 29 (dias);
    De 1 a menos de 12 meses completos: subcampo varia de 01 a 11 (meses); 
    Anos - subcampo varia de 00 a 99;
  - 9 - ignorado

Vamos verificar o primeiro dígito:

# como fator é ruim de manipular
Dados$IDADE <- as.character(Dados$IDADE)

# códigos existentes
todos <- table(Dados$IDADE)
names(todos)

  [1] "001" "002" "003" "004" "005" "006" "007" "008" "009" "010" "011" "012"
 [13] "013" "014" "015" "016" "017" "018" "019" "020" "021" "022" "023" "024"
 [25] "025" "026" "027" "028" "029" "030" "031" "032" "033" "034" "035" "036"
 [37] "037" "038" "039" "040" "041" "042" "043" "044" "045" "046" "047" "048"
 [49] "049" "050" "051" "052" "053" "054" "055" "056" "057" "058" "059" "100"
 [61] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110" "111" "112"
 [73] "113" "114" "115" "116" "117" "118" "119" "120" "121" "122" "123" "200"
 [85] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210" "211" "212"
 [97] "213" "214" "215" "216" "217" "218" "219" "220" "221" "222" "223" "224"
[109] "225" "226" "227" "228" "229" "301" "302" "303" "304" "305" "306" "307"
[121] "308" "309" "310" "311" "400" "401" "402" "403" "404" "405" "406" "407"
[133] "408" "409" "410" "411" "412" "413" "414" "415" "416" "417" "418" "419"
[145] "420" "421" "422" "423" "424" "425" "426" "427" "428" "429" "430" "431"
[157] "432" "433" "434" "435" "436" "437" "438" "439" "440" "441" "442" "443"
[169] "444" "445" "446" "447" "448" "449" "450" "451" "452" "453" "454" "455"
[181] "456" "457" "458" "459" "460" "461" "462" "463" "464" "465" "466" "467"
[193] "468" "469" "470" "471" "472" "473" "474" "475" "476" "477" "478" "479"
[205] "480" "481" "482" "483" "484" "485" "486" "487" "488" "489" "490" "491"
[217] "492" "493" "494" "495" "496" "497" "498" "499" "500" "501" "502" "503"
[229] "504" "505" "506" "507" "508" "509" "510" "511" "512" "513" "514" "515"
[241] "516" "517" "518" "519" "520" "525" "999"

# verificamos o primeiro dígito
# O primeiro, de 1 dígito, indica a unidade da idade
#  (se 1 = minuto, se 2 = hora, se 3 = mês, se 4 = ano,
#   se = 5 idade maior que 100 anos).
# 9 - ignorado

table(substr(Dados$IDADE,1,1))

      0       1       2       3       4       5       9 
   2575    5137   15082    8638 1511205   11800    2387 

Ao que será que 0 corresponde?

Mesmo assim, vamos verificar o que aparece como segundo e terceiro dígitos:

#  1 = minuto, 2 = hora, 3 = mês, 4 = ano, 5 idade maior que 100 anos
# - O segundo, de dois dígitos, indica a quantidade de unidades:
#   (1) Idade menor de 1 hora: subcampo varia de 01 e 59 (minutos);
#   (2) De 1 a 23 Horas: subcampo varia de 01 a 23 (horas);
#   (?) De 24 horas e 29 dias: subcampo varia de 01 a 29 (dias);
#   (3) De 1 a menos de 12 meses completos: subcampo varia de 01 a 11 (meses);
#   (4) Anos - subcampo varia de 00 a 99;
#   (5) > 100 anos

for(i in c(0:5,9))
{
  cat("\nPrimeiro dígito = ",i," (",sep="")
  if(i == 0){cat("não documentado")}
  if(i == 1){cat("minuto, deve variar de 01 a 59")}
  if(i == 2){cat("hora, deve variar de 01 a 23")}
  if(i == 3){cat("mês, deve variar de 01 a 11")}
  if(i == 4){cat("ano, deve variar de 00 a 99")}
  if(i == 5){cat("> 100 anos)")}
  if(i == 9){cat("não informado)")}
  cat(")\n")
  Dadostmp <- Dados$IDADE[substr(Dados$IDADE,1,1)==as.character(i)]
  tabela <- table(substr(Dadostmp,2,3))
  tabela <- data.frame(tabela)
  names(tabela) <- c("Posições [2,3]","Frequência")
  print.data.frame(tabela)
}

Primeiro dígito = 0 (não documentado)
   Posições [2,3] Frequência
1              01        225
2              02         50
3              03         51
4              04         20
5              05        141
6              06         16
7              07         22
8              08         19
9              09         18
10             10        182
11             11         19
12             12         23
13             13         23
14             14         20
15             15        119
16             16         21
17             17         26
18             18         22
19             19         18
20             20        155
21             21         33
22             22         34
23             23         22
24             24         19
25             25         81
26             26         26
27             27         29
28             28         30
29             29         24
30             30        204
31             31         21
32             32         28
33             33         26
34             34         19
35             35         66
36             36         17
37             37         33
38             38         20
39             39         13
40             40        106
41             41         21
42             42         23
43             43         22
44             44         26
45             45         74
46             46         18
47             47         25
48             48         28
49             49         14
50             50         75
51             51         19
52             52         13
53             53         20
54             54         16
55             55         58
56             56         27
57             57         25
58             58         15
59             59         15

Primeiro dígito = 1 (minuto, deve variar de 01 a 59)
   Posições [2,3] Frequência
1              00         53
2              01       1474
3              02        752
4              03        437
5              04        329
6              05        246
7              06        226
8              07        155
9              08        172
10             09        166
11             10        132
12             11        129
13             12        149
14             13         94
15             14         94
16             15         75
17             16         66
18             17         54
19             18         56
20             19         49
21             20         40
22             21         38
23             22         43
24             23        108

Primeiro dígito = 2 (hora, deve variar de 01 a 23)
   Posições [2,3] Frequência
1              00        336
2              01       2674
3              02       2098
4              03       1495
5              04       1107
6              05        870
7              06        708
8              07        619
9              08        484
10             09        457
11             10        438
12             11        348
13             12        323
14             13        312
15             14        263
16             15        274
17             16        222
18             17        201
19             18        198
20             19        183
21             20        183
22             21        179
23             22        155
24             23        170
25             24        158
26             25        151
27             26        128
28             27        131
29             28        126
30             29         91

Primeiro dígito = 3 (mês, deve variar de 01 a 11)
   Posições [2,3] Frequência
1              01       2879
2              02       1384
3              03        988
4              04        722
5              05        601
6              06        502
7              07        400
8              08        345
9              09        304
10             10        284
11             11        229

Primeiro dígito = 4 (ano, deve variar de 00 a 99)
    Posições [2,3] Frequência
1               00          7
2               01       1951
3               02       1173
4               03        789
5               04        673
6               05        624
7               06        523
8               07        495
9               08        451
10              09        489
11              10        565
12              11        588
13              12        636
14              13        857
15              14       1234
16              15       1732
17              16       2451
18              17       3401
19              18       4190
20              19       4868
21              20       5306
22              21       5305
23              22       5462
24              23       5358
25              24       5397
26              25       5425
27              26       5363
28              27       5183
29              28       5110
30              29       5044
31              30       5483
32              31       5841
33              32       5911
34              33       6240
35              34       6522
36              35       6852
37              36       6990
38              37       7707
39              38       8238
40              39       8380
41              40       8823
42              41       9218
43              42       9597
44              43       9736
45              44      10225
46              45      10648
47              46      11172
48              47      11903
49              48      12374
50              49      12792
51              50      13997
52              51      14785
53              52      15977
54              53      16813
55              54      17757
56              55      18999
57              56      20268
58              57      21171
59              58      21646
60              59      22689
61              60      24474
62              61      25278
63              62      26749
64              63      27200
65              64      28436
66              65      29051
67              66      30439
68              67      30949
69              68      31736
70              69      31878
71              70      33134
72              71      32751
73              72      33927
74              73      33515
75              74      32910
76              75      33239
77              76      34197
78              77      34293
79              78      33233
80              79      33534
81              80      36715
82              81      35022
83              82      34032
84              83      33866
85              84      34839
86              85      31939
87              86      29034
88              87      27930
89              88      27590
90              89      26378
91              90      24936
92              91      21315
93              92      19182
94              93      16729
95              94      13918
96              95      11619
97              96       9026
98              97       7227
99              98       5407
100             99       4174

Primeiro dígito = 5 (> 100 anos))
   Posições [2,3] Frequência
1              00       3270
2              01       2325
3              02       1849
4              03       1296
5              04        963
6              05        679
7              06        404
8              07        300
9              08        209
10             09        157
11             10        105
12             11         52
13             12         47
14             13         50
15             14         55
16             15         16
17             16          9
18             17          7
19             18          2
20             19          2
21             20          2
22             25          1

Primeiro dígito = 9 (não informado))
  Posições [2,3] Frequência
1             99       2387

Calculando a idade

Desistimos de entender esta variável, então arrumaremos datas de nascimento e óbito para recalcular a idade ao óbito:

Dados$DTOBITO <- as.Date(as.character(Dados$DTOBITO),
                         format="%d%m%Y")
Dados$DTNASC <- as.Date(as.character(Dados$DTNASC),
                        format="%d%m%Y")
Dados$IdadeObito <- age(Dados$DTNASC, Dados$DTOBITO)
# verificando os valores encontrados
table(Dados$IdadeObito)

    0     1     2     3     4     5     6     7     8     9    10    11    12 
31240  1949  1169   789   672   622   521   495   451   490   564   586   635 
   13    14    15    16    17    18    19    20    21    22    23    24    25 
  856  1233  1729  2448  3397  4176  4853  5292  5298  5453  5341  5397  5395 
   26    27    28    29    30    31    32    33    34    35    36    37    38 
 5360  5176  5110  5031  5434  5830  5906  6229  6517  6809  6984  7701  8229 
   39    40    41    42    43    44    45    46    47    48    49    50    51 
 8369  8762  9210  9593  9723 10222 10607 11169 11893 12367 12785 13925 14776 
   52    53    54    55    56    57    58    59    60    61    62    63    64 
15976 16796 17755 18961 20274 21167 21631 22664 24451 25266 26729 27175 28435 
   65    66    67    68    69    70    71    72    73    74    75    76    77 
29021 30432 30919 31730 31869 33111 32725 33925 33509 32894 33217 34196 34270 
   78    79    80    81    82    83    84    85    86    87    88    89    90 
33218 33508 36713 35012 34019 33838 34839 31920 29030 27909 27590 26369 24920 
   91    92    93    94    95    96    97    98    99   100   101   102   103 
21301 19176 16725 13909 11610  9019  7227  5406  4165  3266  2325  1846  1294 
  104   105   106   107   108   109   110   111   112   113   114   115   116 
  963   679   403   300   209   157   105    52    47    50    55    16     9 
  117   118   119   120   125  2018 
    7     2     2     2     1     1 

# vamos encontrar o registro com idades muito altas
r.aux <- which(Dados$IdadeObito > 130)
Dados[r.aux,c("DTNASC","DTOBITO")]

# A tibble: 1 × 2
  DTNASC     DTOBITO   
  <date>     <date>    
1 0001-12-11 2020-04-11

# e eliminar o registro com dado impossível
Dados$IdadeObito[r.aux] <- NA
# eliminando colunas que estão vazias
summary(Dados[, c(19, 47, 74, 81, 85)])

 ESTABDESCR      CB_PRE        NUDIASOBIN     NUDIASINF      FONTESINF     
 NA's:1556824   NA's:1556824   NA's:1556824   NA's:1556824   NA's:1556824  

Dados <- subset(Dados,
                select=-c(19, 47, 74, 81, 85))
# um jeito de fazer summary turbinado, bom para arquivo grande
sjPlot::view_df(Dados)

Data frame: Dados

ID	Name	Label	Values	Value Labels
1	ORIGEM			1 2
2	TIPOBITO			2
3	DTOBITO
4	HORAOBITO			0000 0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 <… truncated>
5	NATURAL			004 008 010 013 015 016 019 020 021 023 025 026 027 028 029 <… truncated>
6	CODMUNNATU			110000 110001 110002 110003 110004 110005 110006 110007 110008 110009 110010 110011 110012 110013 110014 <… truncated>
7	DTNASC
8	IDADE			<output omitted>
9	SEXO			0 1 2
10	RACACOR			1 2 3 4 5
11	ESTCIV			1 2 3 4 5 9
12	ESC			1 2 3 4 5 9
13	ESC2010			0 1 2 3 4 5 9
14	SERIESCFAL			1 2 3 4 5 6 7 8
15	OCUP			010105 010110 010115 010205 010210 010215 010305 010310 010315 020105 020110 020115 020205 020305 020310 <… truncated>
16	CODMUNRES			110000 110001 110002 110003 110004 110005 110006 110007 110008 110009 110010 110011 110012 110013 110014 <… truncated>
17	LOCOCOR			1 2 3 4 5 6 9
18	CODESTAB			0000001 0000018 0000019 0000030 0000035 0000044 0000046 0000057 0000058 0000059 0000063 0000066 0000084 0000094 0000098 <… truncated>
19	CODMUNOCOR			110001 110002 110003 110004 110005 110006 110007 110008 110009 110010 110011 110012 110013 110014 110015 <… truncated>
20	IDADEMAE			08 10 11 12 13 14 15 16 17 18 19 20 21 22 23 <… truncated>
21	ESCMAE			1 2 3 4 5 9
22	ESCMAE2010			0 1 2 3 4 5 9
23	SERIESCMAE			1 2 3 4 5 6 7 8
24	OCUPMAE			020305 021105 021205 021210 031105 031205 031210 111220 111405 111410 111415 121005 121010 122515 123105 <… truncated>
25	QTDFILVIVO			00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 <… truncated>
26	QTDFILMORT			00 01 02 03 04 05 06 07 08 09 10 11 12 20 23 <… truncated>
27	GRAVIDEZ			1 2 3 9
28	SEMAGESTAC			0 1 10 11 12 13 14 15 16 17 18 19 2 20 21 <… truncated>
29	GESTACAO			1 2 3 4 5 6 9
30	PARTO			1 2 9
31	OBITOPARTO			3 9
32	PESO			0000 0001 0002 0003 0008 0009 0018 0023 0025 0028 0050 0052 0064 0066 0075 <… truncated>
33	TPMORTEOCO			1 2 3 4 5 8 9
34	OBITOGRAV			1 2 9
35	OBITOPUERP			1 2 3 9
36	ASSISTMED			1 2 9
37	EXAME			1 2 9
38	CIRURGIA			1 2 9
39	NECROPSIA			1 2 9
40	LINHAA			A020 A021 A029 A041 A042 A044 A047 A047E871 A048 A049 A049A419 A058 A059 A060 *A065 <… truncated>
41	LINHAB			A014E069 A020 A021 A029 A040 A041 A043 A044 A047 A048 A049 A049E86X A049K567 A049N390 *A052 <… truncated>
42	LINHAC			A010 A020 A021 A029 A039 A040 A041 A044 A046 A047 A047D849 A048 A049 A04X A051 <… truncated>
43	LINHAD			A029 A044 A047 A048 A049 A049K566 A051 A052 A059 A069 A079 A084 A085 A092 A09X <… truncated>
44	LINHAII			A020 A039I10X A047 A047C751 A047D473 A047E149 A047I10X A047I739 A047L010 A047N179 A047N189 A047N289 A047P700 A047Y96X A048Y96X <… truncated>
45	CAUSABAS			A010 A014 A020 A021 A039 A040 A041 A042 A043 A044 A046 A047 A048 A049 A051 <… truncated>
46	COMUNSVOIM			110002 110004 110010 110012 110013 110015 110020 110025 110028 110030 110032 110033 110045 110170 110180 <… truncated>
47	DTATESTADO			01012020 01012021 01022020 01022021 01032020 01032021 01042020 01042021 01052020 01052021 01062020 01062021 01072020 01072021 01082020 <… truncated>
48	CIRCOBITO			1 2 3 4 9
49	ACIDTRAB			1 2 9
50	FONTE			1 2 3 4 9
51	NUMEROLOTE			20190012 20190013 20190015 20190016 20190018 20190020 20190022 20190031 20190032 20190033 20190054 20200001 20200002 20200003 20200004 <… truncated>
52	TPPOS			N S
53	DTINVESTIG			01012020 01012021 01012022 01022020 01022021 01022022 01032020 01032021 01032022 01042020 01042021 01052020 01052021 01062020 01062021 <… truncated>
54	CAUSABAS_O			A009 A014 A020 A021 A039 A040 A041 A042 A043 A044 A046 A047 A048 A049 A051 <… truncated>
55	DTCADASTRO			01012002 01012004 01012006 01012007 01012008 01012009 01012011 01012014 01012020 01012021 01022020 01022021 01022022 01032020 01032021 <… truncated>
56	ATESTANTE			1 2 3 4 5
57	STCODIFICA			N S
58	CODIFICADO			N S
59	VERSAOSIST			2…0 2.2.03 3.2.00 3.2.01 3.2.02 3.2.30 3.2.50
60	VERSAOSCB			2.2 3.2 3.3 3.4 4.880
61	FONTEINV			1 2 3 4 5 6 7 8 9
62	DTRECEBIM			01022021 01022022 01032021 01032022 01042020 01042021 01052020 01052021 01062020 01062021 01072020 01072021 01082020 01082021 01092020 <… truncated>
63	ATESTADO			A039 I10 A153 A162 A162 B200 A169 A419 R688 U049 A46 A46 U049 B182 B182 K746 B24 B24 I10 B24 I10 E149 F03 B24 I671 *B24 J418 N63 <… truncated>
64	DTRECORIGA			01022021 01022022 01032021 01042020 01042021 01052020 01052021 01062020 01062021 01072020 01072021 01082020 01082021 01092020 01092021 <… truncated>
65	CAUSAMAT			O930 O931 O934 O936 O937 O938 O939
66	ESCMAEAGR1			00 01 02 03 04 05 06 07 08 09 10 11 12
67	ESCFALAGR1			00 01 02 03 04 05 06 07 08 09 10 11 12
68	STDOEPIDEM			0 1
69	STDONOVA			0 1
70	DIFDATA			000 001 002 003 004 005 006 007 008 009 010 011 012 013 014 <… truncated>
71	NUDIASOBCO			0 1 10 100 101 102 103 104 105 106 107 108 109 11 110 <… truncated>
72	DTCADINV			01012021 01022021 01022022 01032021 01042020 01042021 01052020 01062020 01062021 01072020 01072021 01082020 01092020 01092021 01102020 <… truncated>
73	TPOBITOCOR			1 2 3 4 5 6 7 8 9
74	DTCONINV			01012020 01012021 01022020 01022021 01022022 01032020 01032021 01032022 01042020 01042021 01052020 01052021 01062020 01062021 01072020 <… truncated>
75	FONTES			SSSSSS SSSSSX SSSSXS SSSSXX SSSXSX SSSXXS SSSXXX SSXSSX SSXSXS SSXSXX SSXXSS SSXXSX SSXXXS SSXXXX SXSSSX <… truncated>
76	TPRESGINFO			1 2 3
77	TPNIVELINV			E M R
78	DTCADINF			01022021 01022022 01032021 01042020 01042021 01052020 01052021 01062020 01062021 01072020 01072021 01082020 01092020 01092021 01102020 <… truncated>
79	MORTEPARTO			1 2 3 9
80	DTCONCASO			01012021 01022020 01022021 01022022 01032020 01032021 01042020 01042021 01052020 01062020 01062021 01072020 01072021 01082020 01082021 <… truncated>
81	ALTCAUSA			1 2
82	CONTADOR			1 10 100 1000 10000 100000 1000000 1000001 1000002 1000003 1000004 1000005 1000006 1000007 1000008 <… truncated>
83	IdadeObito		range: 0-125

Recodificação

Seguindo o manual, substituímos os códigos numéricos nas colunas que usaremos (o mesmo procedimento pode ser feito para todas as outras que quiserem usar):

Dados$TIPOBITO <- factor(Dados$TIPOBITO,
                         levels=c(1,2),
                         labels=c("Fetal",
                                  "Não fetal"))
Dados$LOCOCOR <- factor(Dados$LOCOCOR,
                        levels=c(1:6,9),
                        labels=c("Hospital",
                                 "Outros estabelecimentos de saúde",
                                 "Domicílio",
                                 "Via pública",
                                 "Outros",
                                 "Aldeia indígena",
                                 "Ignorado"))
Dados$SEXO <- factor(Dados$SEXO,
                     levels=0:2,
                     labels=c("Ignorado",
                              "Masculino",
                              "Feminino"))
Dados$RACACOR <- factor(Dados$RACACOR,
                        levels=1:5,
                        labels=c("Branca",
                                 "Preta",
                                 "Amarela",
                                 "Parda",
                                 "Indígena"))
Dados$ESTCIV <- factor(Dados$ESTCIV,
                       levels=c(1:5,9),
                       labels=c("Solteiro",
                                "Casado",
                                "Viúvo",
                                "Separado judicialmente/ divorciado",
                                "União Estável",
                                "Ignorado"))
Dados$ESC2010 <- factor(Dados$ESC2010,
                        levels=c(0:5,9),
                        labels=c("Sem",
                                 "Fundamental I",
                                 "Fundamental II",
                                 "Médio",
                                 "Superior incompleto",
                                 "Superior completo",
                                 "Ignorado"))
Dados$CIRCOBITO <- factor(Dados$CIRCOBITO,
                         levels=c(1:4,9),
                         labels=c("Acidente",
                                  "Suicídio",
                                  "Homicídio",
                                  "Outros",
                                  "Ignorado"))
Dados$ACIDTRAB <- factor(Dados$ACIDTRAB,
                         levels=c(1,2,9),
                         labels=c("Sim",
                                  "Não",
                                  "Ignorado"))
Dados$FONTE <- factor(Dados$FONTE,
                         levels=c(1,2,3,4,9),
                         labels=c("Ocorrência policial",
                                  "Hospital",
                                  "Família",
                                  "Outra",
                                  "Ignorado"))

Salvando a versão recodificada

Salvamos o arquivo modificado em formato .rds

saveRDS(Dados, "DOBR2020.rds")

Este trecho só precisa ser feito uma vez. Com os dados organizados, agora é possível trabalhar a partir do arquivo DOBR2020.rds:

Lendo o arquivo recodificado

Lemos o arquivo preparado, e separamos apenas algumas colunas para exemplificar a análise:

Dados <- readRDS("DOBR2020.rds")
# separando as variáveis de interesse
Variaveis <- c("TIPOBITO","IdadeObito", "SEXO",
               "RACACOR", "ESTCIV", "ESC2010","PESO",
               "LOCOCOR","CIRCOBITO","ACIDTRAB","FONTE")
Dados <- subset(Dados, select=Variaveis)
# verificando a integridade do arquivo
sjPlot::view_df(Dados)

Data frame: Dados

ID	Name	Label	Values	Value Labels
1	TIPOBITO			Fetal Não fetal
2	IdadeObito		range: 0-125
3	SEXO			Ignorado Masculino Feminino
4	RACACOR			Branca Preta Amarela Parda Indígena
5	ESTCIV			Solteiro Casado Viúvo Separado judicialmente/ divorciado União Estável Ignorado
6	ESC2010			Sem Fundamental I Fundamental II Médio Superior incompleto Superior completo Ignorado
7	PESO			0000 0001 0002 0003 0008 0009 0018 0023 0025 0028 0050 0052 0064 0066 0075 <… truncated>
8	LOCOCOR			Hospital Outros estabelecimentos de saúde Domicílio Via pública Outros Aldeia indígena Ignorado
9	CIRCOBITO			Acidente Suicídio Homicídio Outros Ignorado
10	ACIDTRAB			Sim Não Ignorado
11	FONTE			Ocorrência policial Hospital Família Outra Ignorado

Análise de dados faltantes

n.total <- nrow(Dados)
n.completo <- nrow(na.omit(Dados))
n.incompleto <- n.total - n.completo
cat("Numero de casos total = ", n.total, "\n", sep="")

Numero de casos total = 1556824

cat("Numero de casos completos = ", n.completo, 
    " (",round(100*n.completo/n.total,2),"%)\n", sep="")

Numero de casos completos = 0 (0%)

cat("Numero de casos incompletos = ", n.incompleto, 
    " (",round(100*n.incompleto/n.total,2),"%)\n", sep="")

Numero de casos incompletos = 1556824 (100%)

obs.falt <- sum(is.na(Dados))
obs.valid <- sum(!is.na(Dados))
obs.tot <- obs.falt + obs.valid
cat("Numero de observacoes validas = ", obs.valid, 
    " (",round(100*obs.valid/obs.tot,2),"%)\n", sep="")

Numero de observacoes validas = 10989406 (64.17%)

cat("Numero de observacoes faltantes = ", obs.falt, 
    " (",round(100*obs.falt/obs.tot,2),"%)\n", sep="")

Numero de observacoes faltantes = 6135658 (35.83%)

# usando o modo texto (o gráfico fica ruim de ver)
missing <- finalfit::missing_pattern(Dados, plot=FALSE)
contagem <- rownames(missing)
missing <- data.frame(missing)
missing <- cbind(contagem,missing)
names(missing)[ncol(missing)] <- "missing"
rownames(missing) <- NULL
tabela <- knitr::kable(missing, format="html", caption="Análise de dados faltantes")
tabela <- kableExtra::kable_styling(tabela, bootstrap_options = "striped", full_width = FALSE)
tabela <- kableExtra::row_spec(tabela, 0, background = "#b2dfee")
tabela

Análise de dados faltantes

contagem	TIPOBITO	SEXO	LOCOCOR	IdadeObito	RACACOR	ESTCIV	ESC2010	CIRCOBITO	FONTE	ACIDTRAB	PESO	missing
39266	1	1	1	1	1	1	1	1	1	1	0	1
67637	1	1	1	1	1	1	1	1	1	0	0	2
1948	1	1	1	1	1	1	1	1	0	1	0	2
22597	1	1	1	1	1	1	1	1	0	0	0	3
7	1	1	1	1	1	1	1	0	1	0	0	3
1263555	1	1	1	1	1	1	1	0	0	0	0	4
1026	1	1	1	1	1	1	0	1	1	1	0	2
2294	1	1	1	1	1	1	0	1	1	0	0	3
110	1	1	1	1	1	1	0	1	0	1	0	3
1972	1	1	1	1	1	1	0	1	0	0	0	4
48361	1	1	1	1	1	1	0	0	0	0	0	5
343	1	1	1	1	1	0	1	1	1	1	0	2
711	1	1	1	1	1	0	1	1	1	0	0	3
22	1	1	1	1	1	0	1	1	0	1	0	3
247	1	1	1	1	1	0	1	1	0	0	0	4
1	1	1	1	1	1	0	1	0	1	0	0	4
13618	1	1	1	1	1	0	1	0	0	0	0	5
223	1	1	1	1	1	0	0	1	1	1	1	2
1181	1	1	1	1	1	0	0	1	1	1	0	3
175	1	1	1	1	1	0	0	1	1	0	1	3
1083	1	1	1	1	1	0	0	1	1	0	0	4
13	1	1	1	1	1	0	0	1	0	1	1	3
63	1	1	1	1	1	0	0	1	0	1	0	4
120	1	1	1	1	1	0	0	1	0	0	1	4
605	1	1	1	1	1	0	0	1	0	0	0	5
25342	1	1	1	1	1	0	0	0	0	0	1	5
20773	1	1	1	1	1	0	0	0	0	0	0	6
395	1	1	1	1	0	1	1	1	1	1	0	2
627	1	1	1	1	0	1	1	1	1	0	0	3
25	1	1	1	1	0	1	1	1	0	1	0	3
273	1	1	1	1	0	1	1	1	0	0	0	4
20869	1	1	1	1	0	1	1	0	0	0	0	5
45	1	1	1	1	0	1	0	1	1	1	0	3
107	1	1	1	1	0	1	0	1	1	0	0	4
3	1	1	1	1	0	1	0	1	0	1	0	4
63	1	1	1	1	0	1	0	1	0	0	0	5
2603	1	1	1	1	0	1	0	0	0	0	0	6
43	1	1	1	1	0	0	1	1	1	1	0	3
71	1	1	1	1	0	0	1	1	1	0	0	4
5	1	1	1	1	0	0	1	1	0	1	0	4
63	1	1	1	1	0	0	1	1	0	0	0	5
3575	1	1	1	1	0	0	1	0	0	0	0	6
6	1	1	1	1	0	0	0	1	1	1	1	3
90	1	1	1	1	0	0	0	1	1	1	0	4
4	1	1	1	1	0	0	0	1	1	0	1	4
211	1	1	1	1	0	0	0	1	1	0	0	5
17	1	1	1	1	0	0	0	1	0	1	0	5
2	1	1	1	1	0	0	0	1	0	0	1	5
221	1	1	1	1	0	0	0	1	0	0	0	6
2193	1	1	1	1	0	0	0	0	0	0	1	6
8275	1	1	1	1	0	0	0	0	0	0	0	7
120	1	1	1	0	1	1	1	1	1	1	0	2
447	1	1	1	0	1	1	1	1	1	0	0	3
83	1	1	1	0	1	1	1	1	0	0	0	4
1006	1	1	1	0	1	1	1	0	0	0	0	5
4	1	1	1	0	1	1	0	1	1	1	0	3
26	1	1	1	0	1	1	0	1	1	0	0	4
4	1	1	1	0	1	1	0	1	0	0	0	5
61	1	1	1	0	1	1	0	0	0	0	0	6
3	1	1	1	0	1	0	1	1	1	1	0	3
6	1	1	1	0	1	0	1	1	1	0	0	4
2	1	1	1	0	1	0	1	1	0	0	0	5
15	1	1	1	0	1	0	1	0	0	0	0	6
1	1	1	1	0	1	0	0	1	1	1	1	3
59	1	1	1	0	1	0	0	1	1	1	0	4
2	1	1	1	0	1	0	0	1	1	0	1	4
400	1	1	1	0	1	0	0	1	1	0	0	5
2	1	1	1	0	1	0	0	1	0	1	0	5
1	1	1	1	0	1	0	0	1	0	0	1	5
139	1	1	1	0	1	0	0	1	0	0	0	6
132	1	1	1	0	1	0	0	0	0	0	1	6
382	1	1	1	0	1	0	0	0	0	0	0	7
21	1	1	1	0	0	1	1	1	1	1	0	3
102	1	1	1	0	0	1	1	1	1	0	0	4
25	1	1	1	0	0	1	1	1	0	0	0	5
123	1	1	1	0	0	1	1	0	0	0	0	6
1	1	1	1	0	0	1	0	1	1	1	0	4
6	1	1	1	0	0	1	0	1	1	0	0	5
4	1	1	1	0	0	1	0	1	0	0	0	6
6	1	1	1	0	0	1	0	0	0	0	0	7
5	1	1	1	0	0	0	1	1	1	0	0	5
3	1	1	1	0	0	0	1	1	0	0	0	6
4	1	1	1	0	0	0	1	0	0	0	0	7
31	1	1	1	0	0	0	0	1	1	1	0	5
209	1	1	1	0	0	0	0	1	1	0	0	6
1	1	1	1	0	0	0	0	1	0	1	0	6
67	1	1	1	0	0	0	0	1	0	0	0	7
31	1	1	1	0	0	0	0	0	0	0	1	7
211	1	1	1	0	0	0	0	0	0	0	0	8
	0	0	0	3745	40636	81002	118961	1411143	1439835	1511757	1528579	6135658

Estatística descritiva

pie(table(Dados$CIRCOBITO), main="Tipo de óbito por causa externa")

xtabs(~ Dados$CIRCOBITO)

Dados$CIRCOBITO
 Acidente  Suicídio Homicídio    Outros  Ignorado 
    61751     13167     45687      4692     20384 

round(prop.table(xtabs(~ Dados$CIRCOBITO)),2)

Dados$CIRCOBITO
 Acidente  Suicídio Homicídio    Outros  Ignorado 
     0.42      0.09      0.31      0.03      0.14 

pie(table(Dados$FONTE), main="Fonte de informação do óbito por causa externa")

xtabs(~ Dados$FONTE)

Dados$FONTE
Ocorrência policial            Hospital             Família               Outra 
              76362               10956                7476               11355 
           Ignorado 
              10840 

round(prop.table(xtabs(~ Dados$FONTE)),2)

Dados$FONTE
Ocorrência policial            Hospital             Família               Outra 
               0.65                0.09                0.06                0.10 
           Ignorado 
               0.09 

plot(xtabs(~ Dados$FONTE + Dados$CIRCOBITO))

tabela <- xtabs(~ Dados$FONTE + Dados$CIRCOBITO)
print(tabela)

                     Dados$CIRCOBITO
Dados$FONTE           Acidente Suicídio Homicídio Outros Ignorado
  Ocorrência policial    31703     8363     29647   2431     4214
  Hospital                7098      542      1783    366     1167
  Família                 4412     1006      1608    107      342
  Outra                   5209     1157      3761    550      676
  Ignorado                2853      502      2163    119     5202

gplots::balloonplot(t(tabela), main ="Declarações de Óbito\nBrasil 2020", 
                    xlab ="Óbito (causa externa)", 
                    ylab="Fonte de informação",
                    label=TRUE, show.margins=TRUE, 
                    show.zeros=TRUE, 
                    dotcolor="gray75")

Estudando a idade do óbito

# car::densityPlot(Dados$IdadeObito~Dados$SEXO) # demora demais!
print(psych::describeBy(Dados$IdadeObito, Dados$SEXO, mat=1, digits=2))

    item    group1 vars      n  mean    sd median trimmed   mad min max range
X11    1  Ignorado    1    284 21.69 33.99      0   16.15  0.00   0 120   120
X12    2 Masculino    1 871339 63.02 21.27     67   65.11 19.27   0 125   125
X13    3  Feminino    1 681456 70.11 20.34     74   72.63 17.79   0 120   120
     skew kurtosis   se
X11  1.08    -0.57 2.02
X12 -0.88     0.47 0.02
X13 -1.30     2.02 0.02

plot(table(Dados$IdadeObito))

boxplot(Dados$IdadeObito~Dados$SEXO, horizontal=TRUE)

ggplot2::ggplot(Dados, 
                ggplot2::aes(IdadeObito, 
                             fill=SEXO, 
                             colour=SEXO)) +
  ggplot2::geom_density(alpha=0.2) +
  ggplot2::theme_bw()

O número de Ignorados é muito pequeno. Podemos corrigir pelo tamanho dos grupos:

# tamanho dos grupos
descricao <- psych::describeBy(Dados$IdadeObito, Dados$SEXO, mat=1, digits=2)
r.i <- which(descricao$group1=="Ignorado")
r.F <- which(descricao$group1=="Feminino")
r.M <- which(descricao$group1=="Masculino")
n.i <- descricao$n[r.i]
n.F <- descricao$n[r.F]
n.M <- descricao$n[r.M]
# distribuicao densidade de probabilidade
d.i <- density(Dados$IdadeObito[Dados$SEXO=="Ignorado"],na.rm=TRUE)
d.F <- density(Dados$IdadeObito[Dados$SEXO=="Feminino"],na.rm=TRUE)
d.M <- density(Dados$IdadeObito[Dados$SEXO=="Masculino"],na.rm=TRUE)
# escala a densidade pelo tamanho dos grupos
d.i$y <- d.i$y * n.i
d.F$y <- d.F$y * n.F
d.M$y <- d.M$y * n.M
# exibe o grafico
x.min <- min(c(d.i$x,d.F$x,d.M$x),na.rm=TRUE)
x.max <- max(c(d.i$x,d.F$x,d.M$x),na.rm=TRUE)
y.min <- min(c(d.i$y,d.F$y,d.M$y),na.rm=TRUE)
y.max <- max(c(d.i$y,d.F$y,d.M$y),na.rm=TRUE)
plot(NA, main="Distribuição do MCT",
     xlab="Idade do Óbito", ylab="Densidade * n",
     xlim=c(x.min,x.max), ylim=c(y.min,y.max),
     bty="n")
lines(d.i, col="gray60", lwd=1)
lines(d.F, col="orangered3", lwd=2, lty=2)
lines(d.M, col="turquoise4", lwd=2, lty=3)
legend("topright", 
       c("Ignorado","Feminino","Masculino"), 
       col=c("gray60","orangered3","turquoise4"),
       lwd=c(1,2,2), 
       lty=c(1,2,3), 
       box.lwd=0, bg="transparent")  

Alternativas para verificar e descrever o conteúdo das variáveis

Uma outra maneira de ver os dados em conjunto (requer o ambiente de um navegador):

summarytools::st_options(lang="pt")
section_title <- "**DOBR2020**"
summarytools::define_keywords(title.freq=section_title,
                             freq="n")
summarytools::view(summarytools::dfSummary(Dados), 
                 file="./image/Summary01.html")

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary01.html

O resultado está disponível neste link.

summarytools::view(summarytools::freq(Dados[,"SEXO"],
                                      plain.ascii=FALSE,
                                      style="rmarkdown"), 
                 file="./image/Summary02.html")

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary02.html

O resultado está disponível neste link.

summarytools::view(summarytools::freq(Dados[,"SEXO"],
                                      plain.ascii=FALSE,
                                      style="rmarkdown",
                                      report.nas=FALSE), 
                 file="./image/Summary03.html")

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary03.html

O resultado está disponível neste link.

summarytools::view(summarytools::freq(Dados[,"SEXO"],
                                      plain.ascii=FALSE,
                                      style="rmarkdown",
                                      report.nas=FALSE,
                                      totals=FALSE,
                                      cumul=FALSE,
                                      headings=FALSE), 
                 file="./image/Summary04.html")

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary04.html

O resultado está disponível neste link.

summarytools::view(summarytools::freq(Dados,
                                      plain.ascii=FALSE,
                                      style="rmarkdown"), 
                 file="./image/Summary05.html")

Variable(s) ignored: IdadeObito

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary05.html

Output file appended: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary05.html

O resultado está disponível neste link.

summarytools::view(summarytools::freq(Dados,
                                      plain.ascii=FALSE,
                                      style="rmarkdown",
                                      order="freq",
                                      rows=1:5), 
                 file="./image/Summary06.html")

There are only 1 rows to show; higher numbers will be ignored

There are only 3 rows to show; higher numbers will be ignored
There are only 3 rows to show; higher numbers will be ignored

Variable(s) ignored: IdadeObito

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary06.html

Output file appended: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary06.html

O resultado está disponível neste link.

summarytools::view(summarytools::ctable(Dados$SEXO,
                                        Dados$RACACOR,
                                        prop='n',
                                        totals=FALSE,
                                        headings=TRUE,
                                        useNA="no"), 
                 file="./image/Summary07.html")

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary07.html

O resultado está disponível neste link.

summarytools::view(summarytools::ctable(Dados$SEXO,
                                        Dados$RACACOR,
                                        prop="t"), 
                 file="./image/Summary08.html")

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary08.html

O resultado está disponível neste link.

summarytools::view(summarytools::descr(Dados,
                                       transpose=TRUE,
                                       headings=TRUE), 
                 file="./image/Summary09.html")

Non-numerical variable(s) ignored: TIPOBITO, SEXO, RACACOR, ESTCIV, ESC2010, PESO, LOCOCOR, CIRCOBITO, ACIDTRAB, FONTE

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary09.html

O resultado está disponível neste link.

summarytools::view(summarytools::descr(Dados,
                                       stats=c("mean", "sd","n.valid"),
                                       transpose=TRUE,
                                       headings=TRUE), 
                 file="./image/Summary10.html")

Non-numerical variable(s) ignored: TIPOBITO, SEXO, RACACOR, ESTCIV, ESC2010, PESO, LOCOCOR, CIRCOBITO, ACIDTRAB, FONTE

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary10.html

O resultado está disponível neste link.

summarytools::view(stby(Dados,
                        INDICES=list(Dados$SEXO),
                        FUN=summarytools::descr,
                        stats="common",
                        transpose=TRUE), 
                 file="./image/Summary11.html")

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Non-numerical variable(s) ignored: TIPOBITO, SEXO, RACACOR, ESTCIV, ESC2010, PESO, LOCOCOR, CIRCOBITO, ACIDTRAB, FONTE

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary11.html

Non-numerical variable(s) ignored: TIPOBITO, SEXO, RACACOR, ESTCIV, ESC2010, PESO, LOCOCOR, CIRCOBITO, ACIDTRAB, FONTE

Output file appended: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary11.html

O resultado está disponível neste link.

summarytools::view(stby(data=list(Dados$SEXO, Dados$RACACOR),
                        INDICES=list(Dados$ESTCIV),
                        FUN=summarytools::ctable,
                        prop="n",
                        totals=FALSE,
                        headings=TRUE,
                        useNA="no"), 
                 file="./image/Summary12.html")

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary12.html

Output file appended: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary12.html

O resultado está disponível neste link.

summarytools::view(stby(data=list(Dados$SEXO, Dados$RACACOR),
                        INDICES=list(Dados$ESTCIV),
                        FUN=summarytools::ctable,
                        prop="t"), 
                 file="./image/Summary13.html")

NA detected in grouping variable(s); consider using useNA = TRUE

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Error in match.call(f, call): ... used in a situation where it does not exist

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary13.html

Output file appended: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary13.html

O resultado está disponível neste link.

summarytools::view(Dados %>%
                     dplyr::group_by(SEXO) %>%
                     summarytools::descr(stats=c("mean", "sd","n.valid"),
                                         transpose=TRUE,
                                         headings=TRUE), 
                 file="./image/Summary14.html")

Error in table(names(candidates))[["tested"]]: índice fora dos limites

Non-numerical variable(s) ignored: TIPOBITO, RACACOR, ESTCIV, ESC2010, PESO, LOCOCOR, CIRCOBITO, ACIDTRAB, FONTE

Output file written: D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Pratica_R\image\Summary14.html

O resultado está disponível neste link.

Epílogo

https://www.mitsloanreview.com.br/post/cuidado-para-nao-se-afogar-em-numeros

Uma (muito) breve introdução ao R

O R é uma linguagem de programação especialmente desenvolvida para análise estatística e visualização de dados. É muito utilizado por estatísticos, cientistas de dados e pesquisadores para realizar análises complexas e gerar gráficos informativos a partir de conjuntos de dados.

Operador de Atribuição

O operador de atribuição no R é representado pelo símbolo <-. Ele é usado para atribuir valores a variáveis. Por exemplo:

x <- 8
nome <- "João"

Isso cria duas variáveis, x e nome, com os valores atribuídos a elas. O valor 8 é numérico. Palavras são escritas entre aspas, indicando que são strings (algo como um cordão, uma fita, uma corrente, um conjunto encadeado de letras); 8 não é “8”- este último é a string “8” (uma corrente com apenas um elo), tratável como “a”, “b” ou “c”.

Nomes de variáveis precisam iniciar com uma letra, seguidas de outras letras (melhor não usar letras acentuadas ou cedilha), dígitos e ou separadores como ‘.’ e ’_’. São nomes válidos valor, codigo.carro, desvio_padrao. Letras maiúsculas e minúsculas são distinguidas: idade, e Idade são variáveis diferentes. Use nomes que tenham significado para você.

Para ver os conteúdos das variáveis criadas acima, basta usar seu nome seguido de [enter] (na Console):

x

[1] 8

nome

[1] "João"

vetores

Variáveis de um mesmo tipo podem ser agrupadas em vetores, criados com a função c (para lembrar de combine; voltaremos às funções adiante):

numeros <- c(1,2,3)
numeros

[1] 1 2 3

nomes <- c("João","Pedro","Alberto","José","Paulo")
nomes

[1] "João"    "Pedro"   "Alberto" "José"    "Paulo"  

Quando você usa o mesmo nome de uma variável, o novo conteúdo substitui o antigo. Neste exemplo todos os números de 1 a 100 são armazenados na variável numeros

numeros <- 1:100
numeros

  [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
 [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
 [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
 [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
 [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
 [91]  91  92  93  94  95  96  97  98  99 100

Esta propriedade de substituição permite que o valor de uma variável seja atualizada a qualquer momento. Por exemplo, para incrementar o valor de x de 1 unidade:

x

[1] 8

x <- x + 1
x

[1] 9

A linha x <- x + 1 (lida como \(x=x+1\)) não faz sentido matemático. Isto é computação, e este comando deve ser lido como (compute \(x+1\) e atribua o resultado a \(x\)).

operadores

No exemplo que acabamos de ver, + é um operador de soma. Existem muitos outros, incluindo os mais fundamentais: - (subtração), * (multiplicação) e / (divisão). Operam, também, sobre vetores. Se quiséssemos a metade dos valores que estão no vetor numeros:

numeros <- numeros / 2

o que atualiza cada um de seus valores:

numeros

  [1]  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5
 [16]  8.0  8.5  9.0  9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0
 [31] 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0 20.5 21.0 21.5 22.0 22.5
 [46] 23.0 23.5 24.0 24.5 25.0 25.5 26.0 26.5 27.0 27.5 28.0 28.5 29.0 29.5 30.0
 [61] 30.5 31.0 31.5 32.0 32.5 33.0 33.5 34.0 34.5 35.0 35.5 36.0 36.5 37.0 37.5
 [76] 38.0 38.5 39.0 39.5 40.0 40.5 41.0 41.5 42.0 42.5 43.0 43.5 44.0 44.5 45.0
 [91] 45.5 46.0 46.5 47.0 47.5 48.0 48.5 49.0 49.5 50.0

data frames

Um data frame é uma estrutura de dados fundamental no R, similar a uma planilha no Excel. Ele é usado para armazenar dados tabulares, com linhas e colunas. Cada coluna em um data frame pode ter um tipo de dado diferente. Para criar um data frame (a função é data.frame):

Nome <- c("Ana", "Carlos", "Maria", "Pedro", "Maristela")
Idade <- c(25, 30, 28, 25, 32)
Nota <- c(8.5, 7.2, 9.0, 4.5, 3.0)
dados <- data.frame(Nome, Idade, Nota)

Nesse exemplo, dados é um data frame com três colunas: “Nome”, “Idade” e “Nota”. Observe:

dados

# A tibble: 5 × 3
  Nome      Idade  Nota
  <chr>     <dbl> <dbl>
1 Ana          25   8.5
2 Carlos       30   7.2
3 Maria        28   9  
4 Pedro        25   4.5
5 Maristela    32   3  

endereçamento de linha e coluna

Para acessar elementos específicos em um data frame, você pode utilizar o índice da linha e o nome da coluna. Por exemplo:

primeira_nota <- dados[1, "Nota"]
idade_carlos <- dados[2, "Idade"]

A variável primeira_nota recebeu o valor da coluna Nota que existia na primeira linha e idade_carlos o conteúdo da coluna Idade da segunda linha. Confira:

primeira_nota

[1] 8.5

idade_carlos

[1] 30

Como colocamos o Carlos na segunda linha, o nome da variável idade_carlos está ok. Cuidado para escolher nomes de variáveis; o computador não o critica se fizer o seguinte:

idade_carlos <- dados[3, "Idade"]
idade_carlos

[1] 28

O que faz com que o conteúdo de idade_carlos seja a idade da Maria. O computador não ficará confuso (mas você tem grande chance de se atrapalhar).

endereçamento de linha ou coluna inteiras

Para ver uma linha (e.g., a segunda linha), basta não informar as colunas:

dados[2,]

# A tibble: 1 × 3
  Nome   Idade  Nota
  <chr>  <dbl> <dbl>
1 Carlos    30   7.2

e para ver uma coluna inteira (e.g., a segunda coluna), basta não informar as linhas:

dados[,2]

[1] 25 30 28 25 32

uma alternativa é usar o nome da coluna, com $:

dados$Idade

[1] 25 30 28 25 32

é também possível atribuir o conteúdo desta coluna a uma variável (que será um vetor)

idades <- dados$Idade
idades

[1] 25 30 28 25 32

funções

As funções são blocos de código que realizam ações específicas. O R oferece muitas funções embutidas e você também pode criar as suas próprias. Por exemplo, a função mean() calcula a média de um vetor numérico, que lhe é fornecida como um parâmetro:

notas <- c(7.5, 8.0, 6.8, 9.2, 5.5)
media <- mean(x=notas)
media

[1] 7.4

media <- mean(notas)
media

[1] 7.4

Muitas vezes, quando coletamos dados, há valores faltantes (missing values) que são denotados em R como NA (not available, não disponível). Por exemplo, não conseguiremos ver a média se fizermos o seguinte:

notas <- c(7.5, 8.0, 6.8, NA, 5.5)
media <- mean(notas)
media

[1] NA

Isto faz sentido. A média de um conjunto de números contendo um ou mais deles desconhecidos é, também, desconhecida.

No entanto, funções em R recebem parâmetros que você pode consultar lendo a documentação. Neste exemplo, digite ?mean ou help(mean) na Console. Verá que existe o parâmetro na.rm (nome do parâmetro escolhido para fazer lembrar de NA remove), cujo default é FALSE. Para calcular a média removendo NA, caso algum exista, experimente:

media <- mean(x=notas, na.rm=TRUE)
media

[1] 6.95

media <- mean(na.rm=TRUE, x=notas)
media

[1] 6.95

media <- mean(notas, na.rm=TRUE)
media

[1] 6.95

Cada função tem seus parâmetros. No exemplo de dados, suponha que eu quisesse pegar a nota e a idade do Carlos, mas não soubesse em que linha do data frame ele está registrado. Para isso usaremos a função which (‘em qual’ para quem fala português) e a função cat (concatenating) para exibir o resultado na tela com um formato mais bonito:

# o operador == é para comparar
# guardo em linha o resultado de em qual linha o Nome é igual a Carlos?
alguem <- "Carlos"
linha <- which(dados$Nome==alguem)
nota_indiv <- dados$Nota[linha]
idade_indiv <- dados$Idade[linha]
cat(alguem," obteve nota ", nota_indiv,
    " quando tinha ",idade_indiv," anos de idade.\n",
    sep="")

Carlos obteve nota 7.2 quando tinha 30 anos de idade.

Note que os parâmetros de uma função ficam entre parênteses, separados por vírgulas. Troque alguem <- “Carlos” por alguem <- “Maria” ou por qualquer outro nome existente em dados$Nome para relatar suas respectivas notas e idades.

No caso de cat, para melhor legibilidade do código, a escrevemos em várias linhas. As strings e as variáveis são concatenadas antes de exibir e as separações foram feitas com vazio (sep=""). Para ver o efeito deste parâmetro, veja as diferenças aproveitando nosso vetor numeros:

numeros

  [1]  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5
 [16]  8.0  8.5  9.0  9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0
 [31] 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0 20.5 21.0 21.5 22.0 22.5
 [46] 23.0 23.5 24.0 24.5 25.0 25.5 26.0 26.5 27.0 27.5 28.0 28.5 29.0 29.5 30.0
 [61] 30.5 31.0 31.5 32.0 32.5 33.0 33.5 34.0 34.5 35.0 35.5 36.0 36.5 37.0 37.5
 [76] 38.0 38.5 39.0 39.5 40.0 40.5 41.0 41.5 42.0 42.5 43.0 43.5 44.0 44.5 45.0
 [91] 45.5 46.0 46.5 47.0 47.5 48.0 48.5 49.0 49.5 50.0

cat(numeros)

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12 12.5 13 13.5 14 14.5 15 15.5 16 16.5 17 17.5 18 18.5 19 19.5 20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 25.5 26 26.5 27 27.5 28 28.5 29 29.5 30 30.5 31 31.5 32 32.5 33 33.5 34 34.5 35 35.5 36 36.5 37 37.5 38 38.5 39 39.5 40 40.5 41 41.5 42 42.5 43 43.5 44 44.5 45 45.5 46 46.5 47 47.5 48 48.5 49 49.5 50

cat(numeros, sep=",")

0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12,12.5,13,13.5,14,14.5,15,15.5,16,16.5,17,17.5,18,18.5,19,19.5,20,20.5,21,21.5,22,22.5,23,23.5,24,24.5,25,25.5,26,26.5,27,27.5,28,28.5,29,29.5,30,30.5,31,31.5,32,32.5,33,33.5,34,34.5,35,35.5,36,36.5,37,37.5,38,38.5,39,39.5,40,40.5,41,41.5,42,42.5,43,43.5,44,44.5,45,45.5,46,46.5,47,47.5,48,48.5,49,49.5,50

cat(numeros, sep=", ")

0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5, 15, 15.5, 16, 16.5, 17, 17.5, 18, 18.5, 19, 19.5, 20, 20.5, 21, 21.5, 22, 22.5, 23, 23.5, 24, 24.5, 25, 25.5, 26, 26.5, 27, 27.5, 28, 28.5, 29, 29.5, 30, 30.5, 31, 31.5, 32, 32.5, 33, 33.5, 34, 34.5, 35, 35.5, 36, 36.5, 37, 37.5, 38, 38.5, 39, 39.5, 40, 40.5, 41, 41.5, 42, 42.5, 43, 43.5, 44, 44.5, 45, 45.5, 46, 46.5, 47, 47.5, 48, 48.5, 49, 49.5, 50

Outro detalhe é a sequência especial \n para quebra de linha.

cat("Esta é uma string","\n","e esta é outra, após a quebra de linha\n", sep="")

Esta é uma string
e esta é outra, após a quebra de linha

Convém terminar com uma quebra para garantir que o próximo comando não fique grudado no anterior. Observe:

cat("Esta é uma string qualquer")
x <- 9
cat("A raiz quadrada de ",x," = ",sqrt(x), sep="")

Esta é uma string qualquerA raiz quadrada de 9 = 3

e compare o bom hábito de marcar onde as linhas devem quebrar com:

cat("Esta é uma string qualquer\n")
x <- 9
cat("A raiz quadrada de ",x," = ",sqrt(x),"\n", sep="")

Esta é uma string qualquer
A raiz quadrada de 9 = 3

gráficos

Existem muitas funções para criar gráficos em R (voltaremos a isso adiante). (quando chamadas, no ambiente do RStudio, aparecem na aba Plots). Uma função do r-base muito usada é plot Por exemplo:

# Criação de um gráfico de dispersão das idades em relação às notas
plot(dados$Idade, dados$Nota, main="Idade vs. Nota", xlab="Idade", ylab="Nota")

comentários

Comentários em R são linhas que não são executadas, mas servem para você anotar o que está fazendo para referência futura. Basta iniciar a linha com # para que o R não a interprete. Por exemplo:

# Criação de um data frame com informações pessoais
Nome <- c("Ana", "Carlos", "Maria", "Pedro", "Maristela")
Idade <- c(25, 30, 28, 25, 32)
Nota <- c(8.5, 7.2, 9.0, 4.5, 3.0)
dados <- data.frame(Nome, Idade, Nota)
# Exibe o conteúdo de dados
print(dados)

       Nome Idade Nota
1       Ana    25  8.5
2    Carlos    30  7.2
3     Maria    28  9.0
4     Pedro    25  4.5
5 Maristela    32  3.0

um exemplo

Vamos juntar tudo em um exemplo simples. Digamos que você queira calcular a média das idades das pessoas aprovadas (com nota igual ou acima de 5) e reprovadas (abaixo de 5) no nosso data frame dados:

# Criação de um data frame com informações pessoais
Nome <- c("Ana", "Carlos", "Maria", "Pedro", "Maristela")
Idade <- c(25, 30, 28, 25, 32)
Nota <- c(8.5, 7.2, 9.0, 4.5, 3.0)
dados <- data.frame(Nome, Idade, Nota)
# Exibe o conteúdo de dados
print(dados)

       Nome Idade Nota
1       Ana    25  8.5
2    Carlos    30  7.2
3     Maria    28  9.0
4     Pedro    25  4.5
5 Maristela    32  3.0

# separa as notas acima e abaixo de 5
notas_aprov <- dados$Nota[dados$Nota >= 5]
media_aprov <- mean(notas_aprov)
notas_reprov <- dados$Nota[dados$Nota < 5]
media_reprov <- mean(notas_reprov)
cat("Entre os ",length(notas_aprov),
    " aprovados com nota maior ou igual a 5, a nota média foi ",
    round(media_aprov,3),"\n", sep="")

Entre os 3 aprovados com nota maior ou igual a 5, a nota média foi 8.233

cat("Entre os ",length(notas_reprov),
    " reprovados com nota inferior a 5, a nota média foi ",
    round(media_reprov,3),"\n", sep="")

Entre os 2 reprovados com nota inferior a 5, a nota média foi 3.75

# Criação de um gráfico de dispersão das idades em relação às notas
plot(dados$Idade, dados$Nota, 
     main="Idade vs. Nota", xlab="Idade", ylab="Nota",
     xlim=c(20,35), ylim=c(0,10))
# linha horizontal pontilhada em y=5
abline(h=5, lty=2)
# color em azul para aprovados e vermelho para reprovados
points(dados$Idade[dados$Nota>=5],dados$Nota[dados$Nota>=5],pch=19,col="blue")
points(dados$Idade[dados$Nota<5],dados$Nota[dados$Nota<5],pch=19,col="red")
# percorremos linha a linha do data frame, colocando os nomes das pessoas
for(linha in 1:nrow(dados))
{
  text(x=dados$Idade[linha],y=dados$Nota[linha],
       dados$Nome[linha], pos=3)
}

Neste exemplo, primeiro filtramos as idades das pessoas com nota igual ou acima de 5 e depois calculamos a média dessas idades. Depois fizemos o mesmo para quem teve nota abaixo de 5. Finalmente, criamos um gráfico relacionando suas idades com suas notas.

Aqui usamos mais parâmetros e mais funções (algumas ainda não tinham sido apresentadas): print, length, round, abline, points, for, nrow e text.

para saber mais

Isso é uma introdução básica ao R, abordando os conceitos de atribuição, data frames, endereçamento de linha e coluna, e funções. A partir desses fundamentos, você pode explorar mais recursos avançados e complexos da linguagem à medida que avança no aprendizado.

Além de ter o hábito de consultar a documentação dos pacotes (vistos adiante) e funções no próprio R, que trazem exemplos, explore a Web. Há vários tutoriais disponíveis. Ache exemplos e repita-os passo a passo, para ganhar familiaridade com a linguagem.

ChatGPT e GPTs específicos para R da versão ChatGPT 4o auxiliam na produção de script R.