Introdução à Análise de Dados e Modelação Estatística

Lab#1 – Análise exploratória de dados

Autor

Paulo Soares

Data de Publicação

16 de maio de 2025

1 Infeções hospitalares

O primeiro conjunto de dados que vamos explorar é uma amostra de vários hospitais selecionados de um conjunto inicial de 338 hospitais nos EUA incluídos num estudo sobre a ocorrência de infeções hospitalares no período de 1975-76.

O ficheiro encontra-se disponível em:

https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/hospitals.txt (5.2KB).

Cada linha de esse ficheiro tem um número de identificação e fornece informação sobre 8 outras variáveis para um único hospital. As nove colunas são:

Coluna Variável Descrição
1 Identification number —
2 Length of stay Average length of stay of all patients in hospital (in days)
3 Age Average age of patients (in years)
4 Infection risk Average estimated probability of acquiring infection in hospital (in percent)
5 Routine culturing ratio Ratio of number of cultures performed to number of patients without signs or symptoms of hospital-acquired infection, times 100
6 Routine chest X-ray ratio Ratio of number of X-rays performed to number of patients without signs or symptoms of pneumonia, times 100
7 Number of beds Average number of beds in hospital during study period
8 Medical school affiliation 1=Yes, 2=No
9 Region Geographic region, where: 1=NE, 2=NC, 3=S, 4=W

1.1 Leitura e arrumação

Importe o ficheiro de dados para o R e inspecione a data frame criada. Melhore a estrutura dos dados, reparando possíveis incorreções.

Utilize os seguintes nomes para as variáveis: “stay”, “age”, “risk”, “cultures”, “xrays”, “beds”, “medschool” e “region” (excluindo a primeira coluna que é, neste caso, irrelevante).

library(dplyr)
url <- "https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/hospitals.txt"
hosp <- read.table(url)
# Alternativa
# hosp <- readr::read_table(url, col_names = FALSE)

str(hosp)
'data.frame':   113 obs. of  9 variables:
 $ V1: int  1 2 3 4 5 6 7 8 9 10 ...
 $ V2: num  7.13 8.82 8.34 8.95 11.2 ...
 $ V3: chr  "55.7" "58.2" "56.9" "53.7" ...
 $ V4: num  4.1 1.6 2.7 5.6 5.7 5.1 4.6 5.4 4.3 6.3 ...
 $ V5: num  9 3.8 8.1 18.9 34.5 21.9 16.7 60.5 24.4 29.6 ...
 $ V6: num  39.6 51.7 74 122.8 88.9 ...
 $ V7: int  279 80 107 147 180 150 186 640 182 85 ...
 $ V8: int  2 2 2 2 2 2 2 1 2 2 ...
 $ V9: int  4 2 3 4 1 2 3 2 3 1 ...
summary(hosp)
       V1            V2             V3                  V4             V5      
 Min.   :  1   Min.   : 6.70   Length:113         Min.   :1.30   Min.   : 1.6  
 1st Qu.: 29   1st Qu.: 8.34   Class :character   1st Qu.:3.70   1st Qu.: 8.4  
 Median : 57   Median : 9.42   Mode  :character   Median :4.40   Median :14.1  
 Mean   : 57   Mean   : 9.65                      Mean   :4.35   Mean   :15.8  
 3rd Qu.: 85   3rd Qu.:10.47                      3rd Qu.:5.20   3rd Qu.:20.3  
 Max.   :113   Max.   :19.56                      Max.   :7.80   Max.   :60.5  
       V6              V7            V8             V9      
 Min.   : 39.6   Min.   : 29   Min.   :1.00   Min.   :1.00  
 1st Qu.: 69.5   1st Qu.:106   1st Qu.:2.00   1st Qu.:2.00  
 Median : 82.3   Median :186   Median :2.00   Median :2.00  
 Mean   : 81.6   Mean   :252   Mean   :1.85   Mean   :2.36  
 3rd Qu.: 94.1   3rd Qu.:312   3rd Qu.:2.00   3rd Qu.:3.00  
 Max.   :133.5   Max.   :835   Max.   :2.00   Max.   :4.00  
hosp$V1 <- NULL

names(hosp) <- c("stay", "age", "risk", "cultures", "xrays", "beds", "medschool",
                 "region")

# hosp <- read.table(url, na.strings = "x")
hosp$age <- as.numeric(hosp$age)
Warning: NAs introduced by coercion
hosp$medschool <- (hosp$medschool == 1)

hosp$region <- as.factor(hosp$region)
levels(hosp$region) <- c('NE', 'NC', 'S', 'W')

str(hosp)
'data.frame':   113 obs. of  8 variables:
 $ stay     : num  7.13 8.82 8.34 8.95 11.2 ...
 $ age      : num  55.7 58.2 56.9 53.7 56.5 50.9 57.8 45.7 48.2 NA ...
 $ risk     : num  4.1 1.6 2.7 5.6 5.7 5.1 4.6 5.4 4.3 6.3 ...
 $ cultures : num  9 3.8 8.1 18.9 34.5 21.9 16.7 60.5 24.4 29.6 ...
 $ xrays    : num  39.6 51.7 74 122.8 88.9 ...
 $ beds     : int  279 80 107 147 180 150 186 640 182 85 ...
 $ medschool: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ region   : Factor w/ 4 levels "NE","NC","S",..: 4 2 3 4 1 2 3 2 3 1 ...

1.2 Análise exploratória

  1. Em quantas observações falta o valor da variável age?

    summary(hosp$age)
       Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   38.8    50.8    52.9    53.2    56.1    65.9       3 
    sum(is.na(hosp$age))
    [1] 3
    hosp |> filter(is.na(age)) |> nrow()
    [1] 3

  2. Determine a mediana e o MAD da variável risk e a média e o desvio-padrão da variável age por região.

    hosp |> summarise(mediana_risk = median(risk),
                  mad_risk = mad(risk),
                  media_age = mean(age, na.rm = TRUE),
                  dp_age = sd(age, na.rm = TRUE),
                  .by = region)
      region mediana_risk mad_risk media_age dp_age
1      W         4.45    0.593     53.37  5.223
2     NC         4.40    1.112     51.73  4.498
3      S         4.20    1.631     53.62  4.351
4     NE         4.85    1.334     54.06  4.043

  3. Quantos hospitais há em cada região? E quais são as proporções de hospitais universitários nas 4 regiões?

    summary(hosp$region)
    NE NC  S  W 
28 32 37 16 
    with(hosp, table(region))
    region
NE NC  S  W 
28 32 37 16 
    hosp |>
  summarise(n = n(), .by = region) |>
  arrange(region)
      region  n
1     NE 28
2     NC 32
3      S 37
4      W 16
    tab <- with(hosp, table(medschool, region))
tab_prop <- prop.table(tab, margin = 2)
tab_prop
             region
medschool      NE      NC       S       W
    FALSE 0.82143 0.78125 0.91892 0.87500
    TRUE  0.17857 0.21875 0.08108 0.12500
    hosp |>
  summarise(n = n(), .by = c(medschool, region) ) |>
  mutate(freq = n / sum(n), .by = region) |> 
  arrange(region)
      medschool region  n    freq
1     FALSE     NE 23 0.82143
2      TRUE     NE  5 0.17857
3     FALSE     NC 25 0.78125
4      TRUE     NC  7 0.21875
5     FALSE      S 34 0.91892
6      TRUE      S  3 0.08108
7     FALSE      W 14 0.87500
8      TRUE      W  2 0.12500
    as.data.frame(tab_prop)
      medschool region    Freq
1     FALSE     NE 0.82143
2      TRUE     NE 0.17857
3     FALSE     NC 0.78125
4      TRUE     NC 0.21875
5     FALSE      S 0.91892
6      TRUE      S 0.08108
7     FALSE      W 0.87500
8      TRUE      W 0.12500

  4. Determine o número de hospitais não universitários na região Sul em que o risco de infeção é superior ao valor médio dessa variável para todos os hospitais.

    hosp |>
  filter(!medschool & region == "S" & risk > mean(risk)) |>
  nrow()
    [1] 12

  5. Produza uma tabela de dimensões \(2\times 4\) que permita comparar o risco médio entre hospitais universitários e não universitários em cada uma das 4 regiões.

        hosp |>
      summarise(risk = mean(risk), .by = c(region, medschool)) |>
      arrange(region) |> 
      tidyr::pivot_wider(names_from = region, values_from = risk)
    # A tibble: 2 × 5
  medschool    NE    NC     S     W
  <lgl>     <dbl> <dbl> <dbl> <dbl>
1 FALSE       4.7  4.33  3.76  4.39
2 TRUE        5.6  4.63  5.87  4.3

  6. Agrupe os valores da variável stay num número de classes razoável.

    ncla <- nclass.Sturges(hosp$stay)
ncla
    [1] 8
    range(hosp$stay)
    [1]  6.70 19.56
    delta <- (max(hosp$stay) - min(hosp$stay)) / ncla
brk <- seq(min(hosp$stay), max(hosp$stay), by = delta)
brk
    [1]  6.700  8.308  9.915 11.523 13.130 14.738 16.345 17.953 19.560
    table(cut(hosp$stay, breaks = brk, include.lowest = TRUE))
    
 [6.7,8.31] (8.31,9.91] (9.91,11.5] (11.5,13.1] (13.1,14.7] (14.7,16.3] 
         28          44          30           7           2           0 
  (16.3,18]   (18,19.6] 
          1           1

  7. Construa uma nova data frame que contenha todas as variáveis originais mas em que os valores das variáveis numéricas não inteiras estejam estandardizados.

    hosp |>
  mutate(across(where(is.double), scale)) ->
  nova_df
head(nova_df)
          stay     age    risk cultures   xrays beds medschool region
1 -1.31749  0.5669 -0.1901  -0.6637 -2.1705  279     FALSE      W
2 -0.43334  1.1235 -2.0545  -1.1718 -1.5456   80     FALSE     NC
3 -0.68446  0.8341 -1.2341  -0.7517 -0.3939  107     FALSE      S
4 -0.36533  0.1216  0.9286   0.3036  2.1262  147     FALSE      W
5  0.81178  0.7450  1.0032   1.8278  0.3755  180     FALSE     NE
6  0.05843 -0.5018  0.5557   0.5967  0.7938  150     FALSE     NC
    mean(nova_df$stay)
    [1] -1.85e-16
    sd(nova_df$stay)
    [1] 1

  8. Analise as variáveis age e beds do ponto de vista da simetria da distribuição amostral e da possível existência de observações atípicas, com recurso a estatísticas sumárias e a gráficos.

    library(ggplot2)
library(patchwork)

hosp |>
  select(age, beds) |> 
  summary()
          age            beds    
 Min.   :38.8   Min.   : 29  
 1st Qu.:50.8   1st Qu.:106  
 Median :52.9   Median :186  
 Mean   :53.2   Mean   :252  
 3rd Qu.:56.1   3rd Qu.:312  
 Max.   :65.9   Max.   :835  
 NA's   :3                   
    g1 <- ggplot(hosp) +
        geom_histogram(aes(age), bins = ncla, fill = "cornflowerblue")

g2 <- ggplot(hosp) +
        geom_boxplot(aes(age))

g1 + g2
    Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_bin()`).
    Warning: Removed 3 rows containing non-finite outside the scale range
(`stat_boxplot()`).

    g3 <- ggplot(hosp) +
        geom_histogram(aes(beds), bins = ncla, fill = "cornflowerblue",
                       boundary = 0)

g4 <- ggplot(hosp) +
        geom_boxplot(aes(beds))

g3 + g4

  9. Construa alguns gráficos que permitam analisar possíveis relações entre o risco de infeção e as restantes variáveis.

    Código 
    ggplot(hosp) +
  geom_point(aes(beds, risk, size = stay, color = region), alpha = 0.5)

    Código 
    ggplot(hosp) +
  geom_boxplot(aes(region, risk))

    Código 
    ggplot(hosp) +
  geom_boxplot(aes(medschool, risk))

2 Cotações na BVL

Este novo conjunto de dados contém as cotações diárias das empresas EDP, GALP, MOTA_ENGIL, NOS e NOVABASE na Bolsa de Valores de Lisboa entre 11/3/2019 e 8/3/2021.

Os dados encontram-se em https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/, em 5 ficheiros de texto com as mesmas variáveis separadas por tabs, cada um com o mesmo nome de cada uma das empresas e extensão txt, por exemplo:

https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/EDP.txt (37.7KB).

As variáveis são as seguintes:

Variável Descrição
Date Data da sessão
Open Cotação na abertura da sessão
High Cotacão máxima na sessão
Low Cotacão mínima na sessão
Close Cotação no fecho da sessão
Number of Shares Número de ações transacionadas
Number of Trades Número de transações
Turnover Volume de negócios
vwap Volume-weighted average price

2.1 Leitura e arrumação

  1. Importe os ficheiros de dados para o R e coloque os dados numa única data frame.

    localizacao <- 'https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/'
empresas <- c("EDP", "GALP", "MOTA_ENGIL", "NOS", "NOVABASE")
bolsa <- data.frame()

for (empresa in empresas) {
  url <- paste0(localizacao, empresa, ".txt")
  tmp <- read.delim(url, dec = ",")
  tmp$Company <- empresa                  
  bolsa <- rbind(bolsa, tmp)
}
bolsa$Company <- as.factor(bolsa$Company)

  2. Execute o bloco de código abaixo e confirme a obtenção de um gráfico falhado. Investigue a causa do problema e resolva-o.

    library(ggplot2)

bolsa |>
  filter(Company == "EDP") |> 
  ggplot() + 
    geom_line(aes(x = Date, y = Close))
    Código 
    str(bolsa)
    'data.frame':   2544 obs. of  10 variables:
 $ Date            : chr  "11/03/2019" "12/03/2019" "13/03/2019" "14/03/2019" ...
 $ Open            : num  3.16 3.14 3.19 3.23 3.29 ...
 $ High            : num  3.2 3.24 3.23 3.33 3.3 ...
 $ Low             : num  3.16 3.14 3.14 3.23 3.22 ...
 $ Close           : num  3.19 3.17 3.23 3.27 3.24 ...
 $ Number.of.Shares: int  7198947 8439647 9818184 8582235 45278410 8181860 5943336 8596715 8651093 10727283 ...
 $ Number.of.Trades: int  2923 3457 4019 3893 5668 3471 2296 3061 2325 3635 ...
 $ Turnover        : int  22969557 26933090 31397659 28156811 146695679 26724598 19459171 27967858 28799375 35990661 ...
 $ vwap            : num  3.19 3.19 3.2 3.28 3.24 ...
 $ Company         : Factor w/ 5 levels "EDP","GALP","MOTA_ENGIL",..: 1 1 1 1 1 1 1 1 1 1 ...
    Código 
    library(lubridate)
    
Attaching package: 'lubridate'
    The following objects are masked from 'package:base':

    date, intersect, setdiff, union
    Código 
    bolsa$Date <- dmy(bolsa$Date)
# Alternativa
# bolsa$Date <- as.Date(bolsa$Date, format = "%d/%m/%Y")
    Código 
    bolsa |>
  filter(Company == "EDP") |> 
  ggplot() + 
    geom_line(aes(x = Date, y = Close))

2.2 Análise exploratória

  1. Construa um gráfico semelhante ao anterior mas que inclua todas as 5 empresas.

    Código 
    bolsa |>
  ggplot() + 
    geom_line(aes(x = Date, y = Close, color = Company))

  2. Crie duas novas variáveis com:

    1. o volume de negócios médio por transação (vnpt);
    2. a diferença da cotação no fecho das sessões em dias consecutivos (volatilidade).
    Código 
    bolsa |> 
  mutate(vnpt = Turnover / Number.of.Shares) |> 
  mutate(volatilidade = c(NA, diff(Close)), .by = Company) ->
  bolsa

str(bolsa)
    'data.frame':   2544 obs. of  12 variables:
 $ Date            : Date, format: "2019-03-11" "2019-03-12" ...
 $ Open            : num  3.16 3.14 3.19 3.23 3.29 ...
 $ High            : num  3.2 3.24 3.23 3.33 3.3 ...
 $ Low             : num  3.16 3.14 3.14 3.23 3.22 ...
 $ Close           : num  3.19 3.17 3.23 3.27 3.24 ...
 $ Number.of.Shares: int  7198947 8439647 9818184 8582235 45278410 8181860 5943336 8596715 8651093 10727283 ...
 $ Number.of.Trades: int  2923 3457 4019 3893 5668 3471 2296 3061 2325 3635 ...
 $ Turnover        : int  22969557 26933090 31397659 28156811 146695679 26724598 19459171 27967858 28799375 35990661 ...
 $ vwap            : num  3.19 3.19 3.2 3.28 3.24 ...
 $ Company         : Factor w/ 5 levels "EDP","GALP","MOTA_ENGIL",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ vnpt            : num  3.19 3.19 3.2 3.28 3.24 ...
 $ volatilidade    : num  NA -0.0215 0.0557 0.044 -0.0371 ...

  3. Represente a evolução da nova variável volatilidade para as 5 empresas. Comente o gráfico obtido e tente alguma alternativa.

    Código 
    ggplot(bolsa) +
  geom_line(aes(x = Date, y = volatilidade, color = Company))
    Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_line()`).

    Código 
    ggplot(bolsa) +
  geom_line(aes(x = Date, y = volatilidade), color = "darkorange") +
  facet_wrap(~ Company)
    Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).

  4. Selecione os dados relativos a 2020, calcule as médias mensais da variável vwap e represente-as num gráfico de pontos.

    Código 
    bolsa |>
  mutate(Month = month(Date), Year = year(Date)) |>
  filter(Year == 2020) |>
  summarise(Monthly_vwap_mean = mean(vwap),
            .by = c(Company, Month)) |> 
  ggplot() +
    geom_line(aes(Month, Monthly_vwap_mean, color = Company)) + 
    geom_point(aes(Month, Monthly_vwap_mean, color = Company))