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

Projeto computacional

  1. Os pinguins de Palmer

    Ilustração de @allison_horst
    Ilustração de @allison_horst

     

    O ficheiro https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/penguins.csv (14.9KB) contém observações feitas entre 2007 e 2009 em pinguins adultos de três espécies que se encontram em três ilhas do Arquipélago de Palmer na Antártida.

    As variáveis incluídas no conjunto de dados são as seguintes:

    Variável Descrição
    species espécie (Adelie, Chinstrap ou Gentoo)
    island ilha (Biscoe, Dream ou Torgersen)
    bill_len comprimento do bico (mm)
    bill_dep altura do bico (mm)
    flipper_len comprimento das asas (usadas como barbatanas, claro!) (mm)
    body_mass massa corporal (g)
    sex sexo (female ou male)
    year ano do estudo (2007, 2008 ou 2009)

     

    Ilustração de @allison_horst
    Ilustração de @allison_horst

     

    Referência: Gorman, K. B., Williams, T. D. and Fraser, W. R. (2014) Ecological Sexual Dimorphism and Environmental Variability within a Community of Antarctic Penguins (Genus Pygoscelis). PLoS ONE 9, 3, e90081; doi:10.1371/journal.pone.0090081.

    Leia o ficheiro de dados no R e transforme todas as variáveis de tipo character em fatores.

    Em quantas observações há dados em falta?


    Resolução

    library(dplyr)

url <- "https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/penguins.csv"
pinguins <- read.csv(url)
str(pinguins)
    ## 'data.frame':    344 obs. of  8 variables:
##  $ species    : chr  "Adelie" "Adelie" "Adelie" "Adelie" ...
##  $ island     : chr  "Torgersen" "Torgersen" "Torgersen" "Torgersen" ...
##  $ bill_len   : num  39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
##  $ bill_dep   : num  18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
##  $ flipper_len: int  181 186 195 NA 193 190 181 195 193 190 ...
##  $ body_mass  : int  3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
##  $ sex        : chr  "male" "female" "female" NA ...
##  $ year       : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
    pinguins |> 
  mutate(across(where(is.character), as.factor)) ->
  pinguins
str(pinguins)
    ## 'data.frame':    344 obs. of  8 variables:
##  $ species    : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ island     : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ bill_len   : num  39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
##  $ bill_dep   : num  18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
##  $ flipper_len: int  181 186 195 NA 193 190 181 195 193 190 ...
##  $ body_mass  : int  3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
##  $ sex        : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...
##  $ year       : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
    pinguins$num_vars_falta <- rowSums(is.na(pinguins))
solution <- sum(pinguins$num_vars_falta > 0)

    O valor pretendido é 11.


  2. Elimine as observações em que faltem valores em mais do que uma variável. Em seguida, confirme que apenas faltam valores na variável sex e substitua esses valores em falta por “female”.

    Indique o rácio de pinguins fêmeas por macho no conjunto de dados final.


    Resolução

    pinguins |>
  filter(num_vars_falta < 2) ->
  pinguins
pinguins$num_vars_falta <- NULL

summary(pinguins)
    ##       species          island       bill_len       bill_dep   
##  Adelie   :151   Biscoe   :167   Min.   :32.1   Min.   :13.1  
##  Chinstrap: 68   Dream    :124   1st Qu.:39.2   1st Qu.:15.6  
##  Gentoo   :123   Torgersen: 51   Median :44.5   Median :17.3  
##                                  Mean   :43.9   Mean   :17.2  
##                                  3rd Qu.:48.5   3rd Qu.:18.7  
##                                  Max.   :59.6   Max.   :21.5  
##   flipper_len    body_mass        sex           year     
##  Min.   :172   Min.   :2700   female:165   Min.   :2007  
##  1st Qu.:190   1st Qu.:3550   male  :168   1st Qu.:2007  
##  Median :197   Median :4050   NA's  :  9   Median :2008  
##  Mean   :201   Mean   :4202                Mean   :2008  
##  3rd Qu.:213   3rd Qu.:4750                3rd Qu.:2009  
##  Max.   :231   Max.   :6300                Max.   :2009
    # # Por exemplo:
# sexo_em_falta <- "male"

pinguins |>
  tidyr::replace_na(list(sex = sexo_em_falta)) ->
  pinguins

n_femeas <- sum(pinguins$sex == "female")
n_machos <- nrow(pinguins) - n_femeas
solution <- n_femeas / n_machos

    O valor pretendido é 1.036.


  3. Importante

    Use as seguintes linhas de código para ler o ficheiro de dados que deve usar nesta pergunta e em todas as seguintes.

    url <- ‘https://web.tecnico.ulisboa.pt/paulo.soares/iadme/dados/penguins_alt_female.csv’

    pinguins <- read.csv(url, stringsAsFactors = TRUE)

    Nota: atenção às aspas!

    De entre as variáveis numéricas, a que tem o maior valor do coeficiente de variação é:


    1. body_mass
    2. flipper_len
    3. bill_len
    4. bill_dep
    5. year

    Resolução

    if (medida == "coeficiente de variação") {
  pinguins |>
    summarise(across(where(is.numeric), \(x) sd(x) / mean(x))) ->
    res
} else {
  pinguins |>
    summarise(across(where(is.numeric), mad)) ->
    res
}

solution <- ifelse(extremo == "maior",  which.max(res), which.min(res))
    ##   bill_len bill_dep flipper_len body_mass       year
## 1 0.124302  0.11514   0.0699883  0.190862 0.00040695

    A variável é body_mass.


  4. As cinco variáveis numéricas ordenadas por ordem crescente do coeficiente de correlação linear com a variável body_mass são:


    1. year, body_mass, bill_len, bill_dep, flipper_len
    2. flipper_len, year, bill_len, bill_dep, body_mass
    3. flipper_len, year, bill_dep, bill_len, body_mass
    4. bill_len, flipper_len, year, body_mass, bill_dep
    5. bill_dep, year, bill_len, flipper_len, body_mass

    Resolução

    pinguins |> 
  select(where(is.numeric)) ->
  tmp

if(medida == "da covariância"){
  medidas <- var(tmp)
} else {
  medidas <- cor(tmp)
}

medidas[,"body_mass"] |>
  sort(decreasing = (sentido == "decrescente")) ->
  medidas

medidas |>
  names() ->
  solution
    ##    bill_dep        year    bill_len flipper_len   body_mass 
##  -0.4719156   0.0422094   0.5951098   0.8712018   1.0000000

    A ordem é bill_dep, year, bill_len, flipper_len, body_mass.


  5. Construa um único gráfico que permita analisar visualmente a associação das variáveis categóricas sex e species com a variável body_mass.


    Resolução

    library(ggplot2)
theme_set(theme_light())

cores <- c("darkorange", "mediumorchid1", "cyan4")

ggplot(pinguins) +
  geom_boxplot(aes(sex, body_mass, fill = species)) +
  scale_fill_manual(values = cores)


  6. Construa uma única figura com os gráficos de dispersão que permita analisar visualmente a associação de cada uma das variáveis numéricas (excetuando a variável year) com a variável body_mass. Em todos os gráficos, deve ainda ser possível distinguir as três espécies de pinguins.


    Resolução

    library(ggplot2)
theme_set(theme_light())

cores <- c("darkorange", "mediumorchid1", "cyan4")

pinguins |>
  select(where(is.numeric), -year, species) |> 
  tidyr::pivot_longer(-c(body_mass, species), names_to = "Covariável",
                                              values_to = "Valor") |>
  ggplot() +
    geom_point(aes(Valor, body_mass, color = species), alpha = 0.6) + 
    scale_color_manual(values = cores) +
    facet_wrap(~ Covariável, scales = "free_x") +
    theme(legend.position = "top")


  7. Para avaliar o uso das variáveis species e island como covariáveis num modelo de regressão linear múltipla para explicar a variação da variável resposta body_mass, foi construído o seguinte gráfico:

    O gráfico permite afirmar:


    1. que não há associação entre as variáveis e elas deverão ambas ser incluídas no modelo.
    2. que a associação observada entre as variáveis não impede que ambas sejam incluídas no modelo, deixando-se qualquer possível remoção para uma fase posterior.
    3. que a associação entre as variáveis é perfeita e, por isso, elas não poderão ser usadas conjuntamente.
    4. nada. Quantas mais covariáveis, melhor!
    5. que há uma clara associação entre as variáveis e que uma delas poderá revelar-se redundante no modelo.

    Resolução

    Observa-se a existência de associação entre as variáveis uma vez que alguns valores de uma das variáveis determinam univocamente o valor da outra. No entanto, isso não acontece para todos os valores, ou seja, a associação não é perfeita. Nestas condições, ambas as variáveis podem ser incluídas no modelo. Caso a contribuição de uma delas seja reduzida na presença da outra, poderá então ser removida.


  8. No gráfico da esquerda da figura abaixo representa-se a reta de regressão estimada do modelo de regressão simples flipper_len ~ bill_dep. À direita representam-se os resultados do ajustamento de outro modelo de regressão com as mesmas variáveis mas a que se juntou a covariável categórica species.

    O modelo de regressão linear múltipla que foi considerado pode ser definido por várias fórmulas:


    1. flipper_len ~ (bill_dep + species)^2
    2. flipper_len ~ bill_dep * species - bill_dep:species
    3. flipper_len ~ bill_dep * species
    4. flipper_len ~ bill_dep + species
    5. flipper_len ~ bill_dep + species + bill_dep:species

    Resolução

    Uma vez que as retas ajustadas para cada nível da covariável species não são paralelas, o modelo considerado inclui a interação entre as duas covariáveis. Esse modelo pode ser definido pelas fórmulas equivalentes:


  9. Os gráficos seguintes sugerem que a variável year deverá ter pouca influência na variação da variável resposta body_mass.

    Ajuste um modelo de regressão linear múltipla de 1ª ordem sem essa covariável (Modelo 1) e indique a estimativa da diferença da massa corporal média entre um pinguim macho da ilha Dream e um pinguim fêmea da ilha Biscoe com valores comuns nas restantes covariáveis.


    Resolução

    # # Por exemplo:
# ilha <- "Dream"

modelo_1 <- lm(body_mass ~ . - year, data = pinguins)
summary(modelo_1)
    ## 
## Call:
## lm(formula = body_mass ~ . - year, data = pinguins)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -782.2 -167.5   -8.2  183.7  922.3 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1539.86     574.37   -2.68  0.00771 ** 
## speciesChinstrap  -270.18      88.44   -3.06  0.00243 ** 
## speciesGentoo      984.50     136.71    7.20  4.0e-12 ***
## islandDream        -18.35      58.65   -0.31  0.75452    
## islandTorgersen    -26.32      60.18   -0.44  0.66215    
## bill_len            19.98       7.10    2.81  0.00518 ** 
## bill_dep            70.02      19.50    3.59  0.00038 ***
## flipper_len         15.86       2.93    5.42  1.2e-07 ***
## sexmale            381.16      47.47    8.03  1.7e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 290 on 333 degrees of freedom
## Multiple R-squared:  0.872,  Adjusted R-squared:  0.869 
## F-statistic:  284 on 8 and 333 DF,  p-value: <2e-16
    solution <- coef(modelo_1)[["sexmale"]] + coef(modelo_1)[[paste0("island", ilha)]]

    O valor pretendido é 362.803 g.


  10. Ajuste um novo modelo (Modelo 2) removendo a variável island do Modelo 1. Compare os dois modelos recorrendo a um teste-F adequado e indique o valor observado da estatística de teste.


    Resolução

    modelo_2 <- update(modelo_1, . ~ . - island)
anv <- anova(modelo_2, modelo_1)
anv
    ## Analysis of Variance Table
## 
## Model 1: body_mass ~ species + bill_len + bill_dep + flipper_len + sex
## Model 2: body_mass ~ (species + island + bill_len + bill_dep + flipper_len + 
##     sex + year) - year
##   Res.Df      RSS Df Sum of Sq     F Pr(>F)
## 1    335 28012161                          
## 2    333 27995432  2     16729 0.099  0.905
    solution <- anv$F[2]

    O valor pretendido é 0.099.


  11. Considere um terceiro modelo (Modelo 3) adicionando um termo de interação entre as variáveis species e sex ao Modelo 2. Indique o valor da variação observada no coeficiente de determinação ajustado entre este último modelo e o Modelo 1 com 6 casas decimais.


    Resolução

    modelo_3 <- update(modelo_2, . ~ . + species:sex)
summ_1 <- summary(modelo_1)
summ_3 <- summary(modelo_3)
summ_3
    ## 
## Call:
## lm(formula = body_mass ~ species + bill_len + bill_dep + flipper_len + 
##     sex + species:sex, data = pinguins)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -766.0 -168.8   -7.6  167.8  867.9 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -1713.69     559.55   -3.06  0.00237 ** 
## speciesChinstrap          -118.61      84.93   -1.40  0.16345    
## speciesGentoo              961.29     126.07    7.63  2.6e-13 ***
## bill_len                    22.35       6.85    3.26  0.00121 ** 
## bill_dep                    73.53      18.75    3.92  0.00011 ***
## flipper_len                 15.78       2.83    5.58  4.9e-08 ***
## sexmale                    416.07      54.58    7.62  2.6e-13 ***
## speciesChinstrap:sexmale  -356.75      82.40   -4.33  2.0e-05 ***
## speciesGentoo:sexmale       64.04      68.71    0.93  0.35199    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 279 on 333 degrees of freedom
## Multiple R-squared:  0.882,  Adjusted R-squared:  0.879 
## F-statistic:  310 on 8 and 333 DF,  p-value: <2e-16
    if(tipo == "ajustado") {
    solution <- summ_3$adj.r.squared - summ_1$adj.r.squared
} else {
    solution <- summ_3$r.squared - summ_1$r.squared
}

    O valor pretendido é 0.009626.


  12. Para avaliar o Modelo 3 começou-se por produzir os seguintes três gráficos dos resíduos.

    Com base nestes gráficos é razoável afirmar que:


    1. o pressuposto da normalidade não é posto em causa.
    2. o gráfico dos resíduos vs leverage mostra que há uma observação que se pode considerar excessivamente influente no modelo.
    3. a hipótese de variância constante parece ser invalidada.
    4. existe uma tendência que não é corretamente modelada por um preditor linear.
    5. não há qualquer evidência clara que contrarie os pressupostos do modelo.

    Resolução

    Em conjunto, os gráficos sugerem que:


  13. Com base no Modelo 3 calcule um intervalo de predição a um nível \(\gamma = 0.93\) para o valor da massa corporal de um pinguim com as seguintes características:

    ##   species island    sex bill_len bill_dep flipper_len
## 1  Adelie  Dream female     33.9     19.2       172.4

    Indique o comprimento do intervalo obtido.


    Resolução

    # # Por exemplo:
# gama <- 0.95
# novo_pinguim <- data.frame(species = "Chinstrap", island = "Dream",
#                            sex = "female", bill_len = 51.5, bill_dep = 18.8,
#                            flipper_len = 196.4)

tipo <- ifelse(tipo == "confiança", "confidence", "prediction")
intervalo <- predict(modelo_3, novo_pinguim, interval = tipo, level = gama)
intervalo
    ##       fit     lwr     upr
## 1 3175.98 2654.23 3697.73
    solution <- intervalo[3] - intervalo[2]

    O valor pretendido é 1043.5034.


  14. Comece por fixar a semente do gerador de números pseudo-aleatórios em 2015. De seguida, use a função sample para separar o conjunto de dados em dois subconjuntos: treino e teste. Este último conjunto deverá conter 20% das observações.

    Use o conjunto de treino para ajustar um modelo de regressão logística para a variável resposta sex em função de todas as outras variáveis. Remova todas as covariáveis não significativas ao nível de significância de 0.05 e ajuste um novo modelo reduzido.

    Compare os dois modelos recorrendo ao teste de razão de verosimilhanças apropriado e indique o valor-p obtido.


    Resolução

    # # Por exemplo:
# semente <- 2000

set.seed(semente)

n <- nrow(pinguins)
teste.index <- sample(n, size = 0.2 * n)

teste <- pinguins[teste.index, ]
treino <- pinguins[-teste.index, ]

logmodelo <- glm(sex ~ ., family = "binomial", data = treino)
summary(logmodelo)
    ## 
## Call:
## glm(formula = sex ~ ., family = "binomial", data = treino)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.78e+02   6.84e+02   -0.26  0.79449    
## speciesChinstrap -7.43e+00   1.73e+00   -4.28  1.8e-05 ***
## speciesGentoo    -1.02e+01   3.01e+00   -3.39  0.00071 ***
## islandDream      -1.06e-01   8.54e-01   -0.12  0.90118    
## islandTorgersen  -1.62e+00   1.08e+00   -1.50  0.13327    
## bill_len          5.92e-01   1.35e-01    4.37  1.2e-05 ***
## bill_dep          1.46e+00   3.69e-01    3.95  7.7e-05 ***
## flipper_len       6.90e-02   5.89e-02    1.17  0.24189    
## body_mass         5.46e-03   1.17e-03    4.68  2.9e-06 ***
## year              4.78e-02   3.41e-01    0.14  0.88879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 379.71  on 273  degrees of freedom
## Residual deviance: 108.24  on 264  degrees of freedom
## AIC: 128.2
## 
## Number of Fisher Scoring iterations: 7
    logmodelo_R <- update(logmodelo, . ~ . - year - island - flipper_len)

anv <- anova(logmodelo_R, logmodelo, test = "LRT")
anv
    ## Analysis of Deviance Table
## 
## Model 1: sex ~ species + bill_len + bill_dep + body_mass
## Model 2: sex ~ species + island + bill_len + bill_dep + flipper_len + 
##     body_mass + year
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       268      112.3                     
## 2       264      108.2  4    4.086    0.395
    solution <- anv$`Pr(>Chi)`[2]

    O valor pretendido é 0.395.


  15. Recorra ao modelo reduzido anterior para predizer os valores da variável sex das observações no conjunto de teste. Use o valor de corte 0.5 para classificar cada pinguim de acordo com a probabilidade estimada.

    Indique a proporção de pinguins macho corretamente classificados.


    Resolução

    teste$pred_sex <- (predict(logmodelo_R, teste, type = "response") > 0.5)

teste$pred_sex <- as.factor(teste$pred_sex)
levels(teste$pred_sex) <- c("female", "male")

tab <- with(teste, table(sex, pred_sex))
prop <- diag(prop.table(tab, margin = 1))
prop
    ##   female     male 
## 0.941176 0.852941
    solution <- ifelse(sexo == "fêmea", prop[1], prop[2])

    O valor pretendido é 0.8529.