Trabalho III - Análise de Dados Categóricos

Vitor Nickhorn Kroeff

GRR20195798

Questão 1

a)

Independência sugere que os sintomas não estão ligados à gravidade da DAVF, tornando o sistema de classificação menos útil. Dependência, por outro lado, indica que os sintomas estão associados à classificação, validando o sistema, o que é a preferência dos pesquisadores para justificar o sistema de classificação.

b)

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

dados <- 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)
tabela <- data.frame(t(matrix(dados, 7,7)), row.names = lin)
tabela <- `colnames<-`(tabela, col)
kable(tabela)

	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

Teremos o teste como

chisq.test(tabela, correct = FALSE)

## Warning in chisq.test(tabela, correct = FALSE): Aproximação do qui-quadrado
## pode estar incorreta

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

Com base no resultado, podemos interpretar que há uma diferença significativa.

c)

Para este caso, o teste qui-quadrado de Pearson pode ser inadequado, pois pressupõe independência entre as categorias.

Questão 2

a)

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

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)

Box Plot:

par(mfrow = c(1,3))
boxplot ( formula = sugar ~ Shelf , data = cereal2 , ylab = "Sugar", xlab = "Prateleira", 
          pars = list ( outpch = NA), col ="steelblue4") 
stripchart (x = cereal2$sugar ~ cereal2$Shelf , lwd = 2, col = "red4", 
            method = "jitter", vertical = TRUE , pch = 1, add = TRUE)


boxplot ( formula = fat ~ Shelf , data = cereal2 , ylab = "Fat", xlab = "Prateleira", 
          pars = list ( outpch = NA), col ="steelblue4") 
stripchart (x = cereal2$fat ~ cereal2$Shelf , lwd = 2, col = "red4", 
            method = "jitter", vertical = TRUE , pch = 1, add = TRUE)

boxplot ( formula = sodium ~ Shelf , data = cereal2 , ylab = "Sodium", xlab = "Prateleira", 
          pars = list ( outpch = NA), col ="steelblue4") 
stripchart (x = cereal2$sodium ~ cereal2$Shelf , lwd = 2, col = "red4", 
            method = "jitter", vertical = TRUE , pch = 1, add = TRUE)

Coordenadas Paralelas:

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

c)

d)

Construindo o ajuste teremos que

require(nnet)
require(MASS)
require(car)
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

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 inclusão da variável fat não é significativa para o modelo.

e)

Nosso modelo com interação será dado por:

require(nnet)
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

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

Podemos observar que o a interação é não significatica de acordo com o resultado do teste.

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

Podemos observar que a prateleira 3 é a que apresenta maior probabilidade.

g)

require(ggplot2)
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" = "red3", "2" = "royalblue", "3" = "orange3", "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

Carregando os 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")))

kable(dados)

	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

Teste:

chisq.test(dados, correct = FALSE)

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

Existe evidência contra a hipótese de independência.

Transformando os dados:

require(dplyr)

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")
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 podemos observar que, ao nível de 5%, a variável Guilda é significativa para o modelo.