Trabalho III - Análise de Dados Categóricos

Helen Eduarda dos Santos Lourenço

Questão 1

a)

Independência significaria que os sintomas não estão relacionados à gravidade da DAVF, tornando o sistema de classificação menos útil. Dependência indicaria que os sintomas estão associados à classificação, validando o sistema. Os pesquisadores prefeririam dependência para justificar o sistema de classificação.

b)

linhas <- c("Hemorragia", "Hipertensão intracraniana", "Déficit neurológico focal",
            "Convulsões", "Deficiência cardíaca", "Mielopatia", "Sintomas não agressivos")
colunas <- c("1", "2a", "2b", "2a e 2b", "3", "4", "5")

dados <- matrix(c(
    0, 0, 2, 1, 10, 19, 5,
    1, 8, 1, 2, 0, 4, 1,
    0, 0, 0, 6, 8, 2, 0,
    0, 1, 0, 2, 1, 3, 0,
    0, 1, 0, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 6,
    83, 17, 7, 6, 6, 1, 0
), nrow = length(linhas), ncol = length(colunas), byrow = TRUE,
dimnames = list(linhas, colunas))

# Tabela
kable(dados)

	1	2a	2b	2a e 2b	3	4	5
Hemorragia	0	0	2	1	10	19	5
Hipertensão intracraniana	1	8	1	2	0	4	1
Déficit neurológico focal	0	0	0	6	8	2	0
Convulsões	0	1	0	2	1	3	0
Deficiência cardíaca	0	1	0	1	0	0	0
Mielopatia	0	0	0	0	0	0	6
Sintomas não agressivos	83	17	7	6	6	1	0

# Teste
chisq.test(dados, correct = FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  dados
## X-squared = 303.2, df = 36, p-value < 2.2e-16

Há evidência contra a independência (p-valor muito pequeno).

c)

No cenário em que os pacientes podem apresentar mais de um sintoma, o uso do teste qui-quadrado de Pearson pode ser inadequado, pois o mesmo pressupõe a independência entre as categorias das variáveis. Se um paciente pode contribuir para mais de uma célula na tabela, a independência entre as categorias não é satisfeita.

Questão 2

Upload dos dados:

cereal = read.csv("http://leg.ufpr.br/~lucambio/ADC/cereal_dillons.csv")

a)

stand01 <- function (x){(x - min(x)) / ( max(x) - min(x))}
cereal2 <- data.frame (Shelf = cereal$Shelf, sugar = stand01 (x = cereal$sugar_g / cereal$size_g ), 
                       fat = stand01 (x = cereal$fat_g / cereal$size_g ), 
                       sodium = stand01 (x = cereal$sodium_mg / cereal$size_g ))
kable(head(cereal2))

Shelf	sugar	fat	sodium
1	0.6428571	0.000	0.5666667
1	0.1285714	0.000	0.9000000
1	0.1285714	0.000	1.0000000
1	0.1125000	0.675	0.8166667
1	0.7800000	0.360	0.6533333
1	0.6387097	0.000	0.5419355

b)

Boxplots:

Coordenadas paralelas:

parcoord(cereal2, col = cereal2$Shelf, lty = cereal2$Shelf, lwd = 1.8)

c)

Acredito que a variável resposta prateleira poderia ser considerada ordinal se alguns outros fatores fossem levados em consideração. Por exemplo: A prateleira 1 fica mais perto da porta de entrada, seguida pelas prateleiras 2, 3 e 4. Apesar disso, o enunciado sugere uma intenção, por parte do supermercado, de colocar as caixas em prateleiras com maior destaque, indicando que, por algum motivo, existe sim uma escala ordinal entre elas (prateleira 1 tem maior destaque, seguida pelas prateleiras 2, 3 e 4).

d)

# Ajuste
ajuste <- multinom(Shelf ~ sugar + fat + sodium, data = cereal2)

## # weights:  20 (12 variable)
## initial  value 55.451774 
## iter  10 value 37.329384
## iter  20 value 33.775257
## iter  30 value 33.608495
## iter  40 value 33.596631
## iter  50 value 33.595909
## iter  60 value 33.595564
## iter  70 value 33.595277
## iter  80 value 33.595147
## final  value 33.595139 
## converged

summary(ajuste)

## Call:
## multinom(formula = Shelf ~ sugar + fat + sodium, data = cereal2)
## 
## Coefficients:
##   (Intercept)      sugar        fat    sodium
## 2    6.900708   2.693071  4.0647092 -17.49373
## 3   21.680680 -12.216442 -0.5571273 -24.97850
## 4   21.288343 -11.393710 -0.8701180 -24.67385
## 
## Std. Errors:
##   (Intercept)    sugar      fat   sodium
## 2    6.487408 5.051689 2.307250 7.097098
## 3    7.450885 4.887954 2.414963 8.080261
## 4    7.435125 4.871338 2.405710 8.062295
## 
## Residual Deviance: 67.19028 
## AIC: 91.19028

# TRV
Anova(ajuste)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Shelf
##        LR Chisq Df Pr(>Chisq)    
## sugar   22.7648  3  4.521e-05 ***
## fat      5.2836  3     0.1522    
## sodium  26.6197  3  7.073e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A adição da variável fat não é significativa.

e)

# Ajuste com interação
ajuste2 <- multinom(Shelf ~ sugar*fat*sodium, data = cereal2)

## # weights:  36 (24 variable)
## initial  value 55.451774 
## iter  10 value 36.170336
## iter  20 value 31.166546
## iter  30 value 29.963705
## iter  40 value 28.414027
## iter  50 value 27.891712
## iter  60 value 27.763967
## iter  70 value 27.622579
## iter  80 value 27.438263
## iter  90 value 27.015534
## iter 100 value 26.772481
## final  value 26.772481 
## stopped after 100 iterations

summary(ajuste2)

## Call:
## multinom(formula = Shelf ~ sugar * fat * sodium, data = cereal2)
## 
## Coefficients:
##   (Intercept)      sugar       fat     sodium sugar:fat sugar:sodium fat:sodium
## 2   -4.563627   8.944868 22.063003   1.030077  35.60873   -12.250084  -23.75955
## 3   24.498320 -22.248456 35.981865 -27.899087 -17.12487    13.253103  -59.54150
## 4   27.246742 -21.852777  7.298799 -29.106797  41.08251     2.887805  -30.85250
##   sugar:fat:sodium
## 2        -55.88455
## 3         37.71571
## 4        -22.59552
## 
## Std. Errors:
##   (Intercept)    sugar       fat   sodium sugar:fat sugar:sodium fat:sodium
## 2    25.21113 29.72894  96.57821 27.29915  135.1117     31.98647   116.0776
## 3    22.83750 25.81043 101.17670 24.61166  150.1228     26.89827   138.0237
## 4    22.80359 26.00692 100.83444 24.51538  150.6750     28.86631   138.5448
##   sugar:fat:sodium
## 2         158.8091
## 3         212.2222
## 4         217.3953
## 
## Residual Deviance: 53.54496 
## AIC: 101.545

# TRV
Anova(ajuste2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Shelf
##                  LR Chisq Df Pr(>Chisq)    
## sugar             19.2525  3  0.0002424 ***
## fat                6.1167  3  0.1060686    
## sodium            30.8407  3  9.183e-07 ***
## sugar:fat          3.2309  3  0.3573733    
## sugar:sodium       3.0185  3  0.3887844    
## fat:sodium         3.1586  3  0.3678151    
## sugar:fat:sodium   2.5884  3  0.4595299    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

f)

pred <- predict(ajuste, newdata = data.frame(sugar = 0.12, fat = 0.5, sodium = 0.13), type = "probs")
round(pred, 4)

##      1      2      3      4 
## 0.0000 0.0000 0.6012 0.3988

g)

mean_fat <- mean(cereal2$fat)
mean_sodium <- mean(cereal2$sodium)
sugar_values <- seq(0, 1, by = 0.01)
pred <- predict(ajuste, newdata = data.frame(sugar = sugar_values, fat = mean_fat, sodium = mean_sodium), type = "probs")

plot_data <- data.frame(sugar = sugar_values, 
                        prob_shelf1 = pred[,1], prob_shelf2 = pred[,2], 
                        prob_shelf3 = pred[,3], prob_shelf4 = pred[,4])

ggplot(plot_data, aes(x = sugar)) +
  geom_line(aes(y = prob_shelf1, color = "1")) +
  geom_line(aes(y = prob_shelf2, color = "2")) +
  geom_line(aes(y = prob_shelf3, color = "3")) +
  geom_line(aes(y = prob_shelf4, color = "4")) +
  labs(x = "Teor de Açúcar", y = "pi.hat", color = "Prateleira") +
  scale_color_manual(values = c("1" = "#ff5454", "2" = "#0090ff", "3" = "peru", "4" = "limegreen")) +
  theme_minimal()

h)

round(exp(coefficients(ajuste)),4)[,-1]

##    sugar     fat sodium
## 2 14.777 58.2480      0
## 3  0.000  0.5729      0
## 4  0.000  0.4189      0

round(exp(confint(ajuste)),2)

## , , 2
## 
##             2.5 %       97.5 %
## (Intercept)  0.00 330392859.41
## sugar        0.00    294843.41
## fat          0.63      5360.63
## sodium       0.00         0.03
## 
## , , 3
## 
##               2.5 %       97.5 %
## (Intercept) 1184.66 5.728017e+15
## sugar          0.00 7.000000e-02
## fat            0.01 6.511000e+01
## sodium         0.00 0.000000e+00
## 
## , , 4
## 
##              2.5 %       97.5 %
## (Intercept) 825.32 3.751459e+15
## sugar         0.00 1.600000e-01
## fat           0.00 4.676000e+01
## sodium        0.00 0.000000e+00

Questão 3

Dados:

dados <- array(data = c(8, 69, 2, 25, 20, 7, 100,
                        39, 139, 39, 57, 48, 20, 177,
                        48, 131, 50, 115, 16, 19, 190), 
               dim = c(7,3),
               dimnames = list(Guilda = c("Air-Arth", "Leaf-Arth", "Wood-Arth", "Ground-Arth", "Fruit", "Seeds", "Nectar"),
                               Degradação = c("Alta", "Moderada", "Baixa")))

dados

##              Degradação
## Guilda        Alta Moderada Baixa
##   Air-Arth       8       39    48
##   Leaf-Arth     69      139   131
##   Wood-Arth      2       39    50
##   Ground-Arth   25       57   115
##   Fruit         20       48    16
##   Seeds          7       20    19
##   Nectar       100      177   190

Análise:

chisq.test(dados, correct = FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  dados
## X-squared = 68.496, df = 12, p-value = 6.115e-10

Forte evidência contra a independência (existe uma relação significativa entre as guildas e os níveis de degradação do habitat). - Transformação dos dados:

df <- melt(dados, varnames = c("Guilda", "Degradação"))

colnames(df) <- c("Guilda", "Degradação", "Dados")
df$Degradação <- factor(df$Degradação, levels = c("Baixa", "Moderada", "Alta"))
df <- df %>% mutate(Valor = Dados) %>% select(-Dados)

Regressão ordinal:

ajuste.ord  <- polr(Degradação ~ Guilda + Valor, data = df, method = "logistic")

TRV:

O summary não está funcionando na minha máquina (não entendi o motivo), então decidi realizar um TRV.

Anova(ajuste.ord)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Degradação
##        LR Chisq Df Pr(>Chisq)    
## Guilda   15.937  6     0.0141 *  
## Valor    18.647  1  1.573e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aqui, fica claro que a inclusão da variável guildas no modelo é significativa a um nível de 5%.

Regressão orientada pelo professor:

df$Guilda.Geral <- with(df, ifelse(grepl("Arth", Guilda), "Arth", "NotArth"))
ajuste.guildas <- polr(Degradação ~ Guilda.Geral + Valor, data = df,method = "logistic")

Summary:

Aqui podemos ver que a guilda NotArth possui um \(\hat\beta\) bem pequeno e um erro padrão elevado, o que nos remete a ideia de que \(\beta = 0\) (variável Guilda.Geral não tem significância estatística no modelo).