stats::glm
y tabla de devianza con car::Anova
summary
emmeans
(ambas escalas)prop.table
doBy
nAGQ
library(tidyverse)
library(magrittr)
library(glue)
library(patchwork)
library(broom)
library(ggforce)
library(ggrepel)
library(paletteer)
library(latex2exp)
library(gt)
library(naniar)
library(ggpointdensity)
library(emmeans)
library(sjPlot)
library(glmmTMB)
library(GGally)
library(lme4)
Factores cruzados vs anidados:
GLS = Generalized least squares
Cuando la misma unidad experimental es sometida a mediciones sucesivas: las observaciones ya no son independeintes, están correlacionadas para cada UE (por ejemplo, para el mismo paciente):
library(nlme)
##
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
##
## lmList
## The following object is masked from 'package:glue':
##
## collapse
## The following object is masked from 'package:dplyr':
##
## collapse
asma <- read_csv("clases.R/asma.csv")
## Parsed with column specification:
## cols(
## paciente = col_double(),
## droga = col_character(),
## tiempo = col_character(),
## vef = col_double()
## )
asma_wide <- pivot_wider(asma, names_from = tiempo, values_from = vef)
GGally::ggpairs(asma_wide, progress = F)
## `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`.
## `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`.
tiempos <- asma_wide[3:11]
round(cov(tiempos), 2)
## basal 1h 2h 3h 4h 5h 6h 7h 8h
## basal 0.28 0.25 0.27 0.26 0.23 0.29 0.26 0.25 0.24
## 1h 0.25 0.58 0.56 0.53 0.48 0.50 0.43 0.38 0.42
## 2h 0.27 0.56 0.60 0.55 0.52 0.54 0.45 0.41 0.45
## 3h 0.26 0.53 0.55 0.56 0.52 0.53 0.45 0.41 0.45
## 4h 0.23 0.48 0.52 0.52 0.54 0.52 0.43 0.41 0.44
## 5h 0.29 0.50 0.54 0.53 0.52 0.60 0.50 0.47 0.50
## 6h 0.26 0.43 0.45 0.45 0.43 0.50 0.49 0.44 0.46
## 7h 0.25 0.38 0.41 0.41 0.41 0.47 0.44 0.49 0.45
## 8h 0.24 0.42 0.45 0.45 0.44 0.50 0.46 0.45 0.50
round(cor(tiempos), 2)
## basal 1h 2h 3h 4h 5h 6h 7h 8h
## basal 1.00 0.63 0.65 0.65 0.60 0.69 0.69 0.68 0.64
## 1h 0.63 1.00 0.95 0.93 0.86 0.84 0.80 0.70 0.77
## 2h 0.65 0.95 1.00 0.95 0.92 0.90 0.83 0.76 0.82
## 3h 0.65 0.93 0.95 1.00 0.95 0.91 0.85 0.79 0.84
## 4h 0.60 0.86 0.92 0.95 1.00 0.91 0.84 0.80 0.85
## 5h 0.69 0.84 0.90 0.91 0.91 1.00 0.93 0.86 0.92
## 6h 0.69 0.80 0.83 0.85 0.84 0.93 1.00 0.90 0.93
## 7h 0.68 0.70 0.76 0.79 0.80 0.86 0.90 1.00 0.90
## 8h 0.64 0.77 0.82 0.84 0.85 0.92 0.93 0.90 1.00
GGally::ggcorr(tiempos, label = T, label_size = 3, limits = F, midpoint = 0, label_round = 2)
asma_wide_2 <-
asma_wide %>%
pivot_longer(cols = matches("\\dh$"), names_to = "tiempo", values_to = "vef")
formulita <- as.formula("vef ~ droga * tiempo + basal")
# Modelo de simetría compuesta
asma_m1 <- gls(formulita, correlation = corCompSymm(form = ~ 1 | paciente), asma_wide_2)
# Modelo de simetria compuesta con varianzas distintas
asma_m2 <- gls(formulita, correlation = corCompSymm(form = ~ 1|paciente), asma_wide_2,
weights = varIdent(form = ~ 1|tiempo))
# Modelo autorregresivo con varianzas iguales
asma_m3 <- gls(formulita, correlation = corAR1(form = ~ 1|paciente), asma_wide_2)
# Modelo autorregresivo con varianzas distintas
asma_m4 <- gls(formulita, correlation = corAR1(form = ~ 1|paciente), asma_wide_2,
weights = varIdent(form = ~ 1|tiempo))
# Matriz desestructurada
asma_m5 <- gls(formulita, correlation = corSymm(form = ~ 1|paciente), asma_wide_2)
El modelo con varianzas iguales estima rho, el de varianzas distintas estima rho y las varianzas, etc.
asma_m1
## Generalized least squares fit by REML
## Model: formulita
## Data: asma_wide_2
## Log-restricted-likelihood: -173.6451
##
## Coefficients:
## (Intercept) drogaB drogaP tiempo2h
## 1.07964094 0.21844510 -0.64440732 -0.07666667
## tiempo3h tiempo4h tiempo5h tiempo6h
## -0.28958333 -0.42708333 -0.42000000 -0.49416667
## tiempo7h tiempo8h basal drogaB:tiempo2h
## -0.60500000 -0.61541667 0.90285162 0.01208333
## drogaP:tiempo2h drogaB:tiempo3h drogaP:tiempo3h drogaB:tiempo4h
## 0.14250000 0.17583333 0.36083333 0.17958333
## drogaP:tiempo4h drogaB:tiempo5h drogaP:tiempo5h drogaB:tiempo6h
## 0.47166667 -0.02166667 0.36166667 -0.11166667
## drogaP:tiempo6h drogaB:tiempo7h drogaP:tiempo7h drogaB:tiempo8h
## 0.48291667 -0.10916667 0.56291667 -0.06541667
## drogaP:tiempo8h
## 0.52083333
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | paciente
## Parameter estimate(s):
## Rho
## 0.7656626
## Degrees of freedom: 576 total; 551 residual
## Residual standard error: 0.5190225
getVarCov(asma_m1)
## Marginal variance covariance matrix
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0.26938 0.20626 0.20626 0.20626 0.20626 0.20626 0.20626 0.20626
## [2,] 0.20626 0.26938 0.20626 0.20626 0.20626 0.20626 0.20626 0.20626
## [3,] 0.20626 0.20626 0.26938 0.20626 0.20626 0.20626 0.20626 0.20626
## [4,] 0.20626 0.20626 0.20626 0.26938 0.20626 0.20626 0.20626 0.20626
## [5,] 0.20626 0.20626 0.20626 0.20626 0.26938 0.20626 0.20626 0.20626
## [6,] 0.20626 0.20626 0.20626 0.20626 0.20626 0.26938 0.20626 0.20626
## [7,] 0.20626 0.20626 0.20626 0.20626 0.20626 0.20626 0.26938 0.20626
## [8,] 0.20626 0.20626 0.20626 0.20626 0.20626 0.20626 0.20626 0.26938
## Standard Deviations: 0.51902 0.51902 0.51902 0.51902 0.51902 0.51902 0.51902 0.51902
asma_m2
## Generalized least squares fit by REML
## Model: formulita
## Data: asma_wide_2
## Log-restricted-likelihood: -170.9409
##
## Coefficients:
## (Intercept) drogaB drogaP tiempo2h
## 1.03114955 0.21889185 -0.64381670 -0.07666667
## tiempo3h tiempo4h tiempo5h tiempo6h
## -0.28958333 -0.42708333 -0.42000000 -0.49416667
## tiempo7h tiempo8h basal drogaB:tiempo2h
## -0.60500000 -0.61541667 0.92102453 0.01208333
## drogaP:tiempo2h drogaB:tiempo3h drogaP:tiempo3h drogaB:tiempo4h
## 0.14250000 0.17583333 0.36083333 0.17958333
## drogaP:tiempo4h drogaB:tiempo5h drogaP:tiempo5h drogaB:tiempo6h
## 0.47166667 -0.02166667 0.36166667 -0.11166667
## drogaP:tiempo6h drogaB:tiempo7h drogaP:tiempo7h drogaB:tiempo8h
## 0.48291667 -0.10916667 0.56291667 -0.06541667
## drogaP:tiempo8h
## 0.52083333
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | paciente
## Parameter estimate(s):
## Rho
## 0.7676762
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | tiempo
## Parameter estimates:
## 1h 2h 3h 4h 5h 6h 7h 8h
## 1.000000 1.029786 1.001402 1.096473 1.071530 1.034629 1.128785 1.116511
## Degrees of freedom: 576 total; 551 residual
## Residual standard error: 0.4895412
getVarCov(asma_m2)
## Marginal variance covariance matrix
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0.23965 0.18945 0.18423 0.20172 0.19713 0.19034 0.20767 0.20541
## [2,] 0.18945 0.25414 0.18972 0.20773 0.20301 0.19601 0.21385 0.21153
## [3,] 0.18423 0.18972 0.24032 0.20201 0.19741 0.19061 0.20796 0.20570
## [4,] 0.20172 0.20773 0.20201 0.28812 0.21615 0.20871 0.22770 0.22523
## [5,] 0.19713 0.20301 0.19741 0.21615 0.27516 0.20396 0.22252 0.22010
## [6,] 0.19034 0.19601 0.19061 0.20871 0.20396 0.25654 0.21486 0.21252
## [7,] 0.20767 0.21385 0.20796 0.22770 0.22252 0.21486 0.30535 0.23186
## [8,] 0.20541 0.21153 0.20570 0.22523 0.22010 0.21252 0.23186 0.29875
## Standard Deviations: 0.48954 0.50412 0.49023 0.53677 0.52456 0.50649 0.55259 0.54658
anova(asma_m1, asma_m2, asma_m3, asma_m4, asma_m5)
## Model df AIC BIC logLik Test L.Ratio p-value
## asma_m1 1 27 401.2902 517.7070 -173.64508
## asma_m2 2 34 409.8817 556.4807 -170.94086 1 vs 2 5.40844 0.6102
## asma_m3 3 27 329.0350 445.4519 -137.51752 2 vs 3 66.84668 <.0001
## asma_m4 4 34 324.5664 471.1653 -128.28318 3 vs 4 18.46869 0.0100
## asma_m5 5 54 261.4179 494.2516 -76.70896 4 vs 5 103.14842 <.0001
Conservamos el modelo 5, menor AIC y p value
clase 11, p. 25
asma_ajuste <- as.data.frame(cbind(
"residuos" = residuals(asma_m5), # Adriana pone que son de Pearson estos
"predichos" = predict(asma_m5)
))
asma_ajuste$droga <- asma_wide_2$droga
ggplot(asma_ajuste) +
aes(predichos, residuos) +
geom_hline(yintercept = 0) +
geom_point()
ggplot(asma_ajuste) +
aes(droga, residuos) +
geom_boxplot()
ggplot(asma_ajuste) +
aes(sample = residuos) +
geom_qq(shape = 1) +
geom_qq_line()
# e<-resid(m5) # residuos de pearson
# pre<-predict(m5) #predichos
# par(mfrow = c(2, 2))
# plot(pre, e, xlab="Predichos", ylab="Residuos de pearson",main="Gr?fico de dispersi?n de RE vs PRED",cex.main=.8 )
# plot(bd2$droga, e, xlab="Predichos", ylab="Residuos de pearson",main="Gr?fico de dispersi?n de RE vs Trat",cex.main=.8 )
# abline(0,0)
# qqnorm(e, cex.main=.8)
# qqline(e)
# par(mfrow = c(1, 1))
# shapiro.test(e)
# hist(e)
library(emmeans)
anova(asma_m5)
## Denom. DF: 551
## numDF F-value p-value
## (Intercept) 1 3526.670 <.0001
## droga 2 9.776 1e-04
## tiempo 7 13.251 <.0001
## basal 1 82.045 <.0001
## droga:tiempo 14 3.963 <.0001
# CLD(emmeans(asma_m5, pairwise ~ droga | tiempo))
Clase 11, diapo 23.
Se modela la respuesta de cada individuo a cada tiempo, con cada tiempo k con su propio efecto fijo (!) no correlacionado con el resto de los tiempos. La correlación se incorpora con el efecto aleatorio de paciente anidado en un tratamiento específico. Los tratamientos tienen efectos fijos y los tiempos también, y puede haber interacción entre tratamiento y tiempo.
Implícitamente, se induce la matriz de simetría compuesta (todos los pacientes con la misma varianza, el mismo rho en la no diagonal):
“No es muy aplicable a DMR”, si bien los resultados para los efectos fijos son los mismos.
asma_m6 <- nlme::lme(vef ~ droga * tiempo + basal, random = ~ 1|paciente, asma_wide_2)
asma_m6
## Linear mixed-effects model fit by REML
## Data: asma_wide_2
## Log-restricted-likelihood: -173.6451
## Fixed: vef ~ droga * tiempo + basal
## (Intercept) drogaB drogaP tiempo2h
## 1.07964094 0.21844510 -0.64440732 -0.07666667
## tiempo3h tiempo4h tiempo5h tiempo6h
## -0.28958333 -0.42708333 -0.42000000 -0.49416667
## tiempo7h tiempo8h basal drogaB:tiempo2h
## -0.60500000 -0.61541667 0.90285162 0.01208333
## drogaP:tiempo2h drogaB:tiempo3h drogaP:tiempo3h drogaB:tiempo4h
## 0.14250000 0.17583333 0.36083333 0.17958333
## drogaP:tiempo4h drogaB:tiempo5h drogaP:tiempo5h drogaB:tiempo6h
## 0.47166667 -0.02166667 0.36166667 -0.11166667
## drogaP:tiempo6h drogaB:tiempo7h drogaP:tiempo7h drogaB:tiempo8h
## 0.48291667 -0.10916667 0.56291667 -0.06541667
## drogaP:tiempo8h
## 0.52083333
##
## Random effects:
## Formula: ~1 | paciente
## (Intercept) Residual
## StdDev: 0.4541552 0.2512505
##
## Number of Observations: 576
## Number of Groups: 72
# Sería igual que
asma_m1
## Generalized least squares fit by REML
## Model: formulita
## Data: asma_wide_2
## Log-restricted-likelihood: -173.6451
##
## Coefficients:
## (Intercept) drogaB drogaP tiempo2h
## 1.07964094 0.21844510 -0.64440732 -0.07666667
## tiempo3h tiempo4h tiempo5h tiempo6h
## -0.28958333 -0.42708333 -0.42000000 -0.49416667
## tiempo7h tiempo8h basal drogaB:tiempo2h
## -0.60500000 -0.61541667 0.90285162 0.01208333
## drogaP:tiempo2h drogaB:tiempo3h drogaP:tiempo3h drogaB:tiempo4h
## 0.14250000 0.17583333 0.36083333 0.17958333
## drogaP:tiempo4h drogaB:tiempo5h drogaP:tiempo5h drogaB:tiempo6h
## 0.47166667 -0.02166667 0.36166667 -0.11166667
## drogaP:tiempo6h drogaB:tiempo7h drogaP:tiempo7h drogaB:tiempo8h
## 0.48291667 -0.10916667 0.56291667 -0.06541667
## drogaP:tiempo8h
## 0.52083333
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | paciente
## Parameter estimate(s):
## Rho
## 0.7656626
## Degrees of freedom: 576 total; 551 residual
## Residual standard error: 0.5190225
La regresión Poisson se usa cuando la VR es un conteo de cosas en un continuo.
malezas <- read.csv("Clases 13 y 14 Regresión Poisson/Malezas.csv")
malezas
## Tratamiento malezas
## 1 herb1 0
## 2 herb1 4
## 3 herb1 1
## 4 herb1 4
## 5 herb1 3
## 6 herb1 4
## 7 herb2 7
## 8 herb2 7
## 9 herb2 6
## 10 herb2 5
## 11 herb2 1
## 12 herb2 6
## 13 control 6
## 14 control 11
## 15 control 12
## 16 control 14
## 17 control 16
## 18 control 18
pinos <- read_csv("Clases 13 y 14 Regresión Poisson/pinos.csv")
## Parsed with column specification:
## cols(
## parcela = col_character(),
## n_conos = col_double(),
## altura = col_double(),
## edad = col_double(),
## area_basal = col_double()
## )
summary(malezas)
## Tratamiento malezas
## control:6 Min. : 0.000
## herb1 :6 1st Qu.: 4.000
## herb2 :6 Median : 6.000
## Mean : 6.944
## 3rd Qu.:10.000
## Max. :18.000
El summary te da el dato de missing values:
summary(pinos)
## parcela n_conos altura edad
## Length:42 Min. : 1.00 Min. : 4.000 Min. :13.00
## Class :character 1st Qu.: 56.25 1st Qu.: 5.550 1st Qu.:16.00
## Mode :character Median :172.00 Median : 7.050 Median :18.00
## Mean :272.57 Mean : 7.217 Mean :18.19
## 3rd Qu.:509.00 3rd Qu.: 8.367 3rd Qu.:20.00
## Max. :919.00 Max. :11.150 Max. :24.00
##
## area_basal
## Min. : 3421
## 1st Qu.: 19737
## Median : 47724
## Mean : 59839
## 3rd Qu.: 89350
## Max. :195565
## NA's :6
fig_missing <-
visdat::vis_miss(pinos) +
theme(
axis.text.x = element_text(angle = 0, size = 8, hjust = .5),
legend.position = "top"
)
fig_nconos <- ggplot(pinos, aes(y = area_basal, x = n_conos)) + geom_miss_point() + guides(color = F)
fig_altura <- ggplot(pinos, aes(y = area_basal, x = altura)) + geom_miss_point() + guides(color = F)
fig_edad <- ggplot(pinos, aes(y = area_basal, x = edad)) + geom_miss_point() + guides(color = F)
fig_missing / (fig_nconos + fig_altura + fig_edad + plot_layout(ncol = 3))
En el TP, Adriana descarta la variable area_basal completa.
malezas %$%
psych::describeBy(malezas, Tratamiento, mat = T) %>%
mutate_if(is.double, round, 3)
## item group1 vars n mean sd median trimmed mad min max range
## 1 1 control 1 6 12.833 4.215 13.0 12.833 3.707 6 18 12
## 2 2 herb1 1 6 2.667 1.751 3.5 2.667 0.741 0 4 4
## 3 3 herb2 1 6 5.333 2.251 6.0 5.333 1.483 1 7 6
## skew kurtosis se
## 1 -0.344 -1.412 1.721
## 2 -0.510 -1.799 0.715
## 3 -1.046 -0.608 0.919
Atención a la opción que considera missing values:
cor(pinos[, -1]) %>% round(2)
## n_conos altura edad area_basal
## n_conos 1.00 0.64 0.63 NA
## altura 0.64 1.00 0.77 NA
## edad 0.63 0.77 1.00 NA
## area_basal NA NA NA 1
cor(pinos[, -1], use = "pairwise.complete.obs") %>% round(2)
## n_conos altura edad area_basal
## n_conos 1.00 0.64 0.63 0.70
## altura 0.64 1.00 0.77 0.87
## edad 0.63 0.77 1.00 0.81
## area_basal 0.70 0.87 0.81 1.00
visdat::vis_cor(pinos[, -1]) +
geom_text(aes(label = round(value, 2)), color = "white")
fig_box <-
ggplot(malezas, aes(Tratamiento, malezas)) +
geom_boxplot(
width = .15,
position = position_nudge(x = -.25),
outlier.shape = 4
) +
geom_jitter(width = .05)
fig_dens <-
ggplot(malezas, aes(malezas)) +
geom_density(aes(fill = factor(Tratamiento)), alpha = 0.8) +
xlab("Número de malezas/m2") +
theme(legend.position = "top")
fig_box + fig_dens
malezas_lm <- lm(malezas ~ Tratamiento, data = malezas)
summary(malezas_lm)
##
## Call:
## lm(formula = malezas ~ Tratamiento, data = malezas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.8333 -1.4583 0.6667 1.3333 5.1667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.833 1.200 10.699 2.04e-08 ***
## Tratamientoherb1 -10.167 1.696 -5.993 2.46e-05 ***
## Tratamientoherb2 -7.500 1.696 -4.421 0.000496 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.938 on 15 degrees of freedom
## Multiple R-squared: 0.7203, Adjusted R-squared: 0.683
## F-statistic: 19.31 on 2 and 15 DF, p-value: 7.089e-05
anova(malezas_lm)
## Analysis of Variance Table
##
## Response: malezas
## Df Sum Sq Mean Sq F value Pr(>F)
## Tratamiento 2 333.44 166.722 19.311 7.089e-05 ***
## Residuals 15 129.50 8.633
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Atención a los residuos con forma de embudo: aumenta la varianza cuando aumenta la media de Y, como en una Poisson. Se viola el supuesto de heterocedasticidad. El QQ se ve feo también.
Adriana plotea los residuos no estandarizados en el Y. Ojo con que en el X los predichos estén en la escala de la VR.
library(patchwork)
malezas_fit <- broom::augment(malezas_lm, data = malezas, type.predict = "response", type.residual = "pearson")
resid_fig <-
malezas_fit %>%
ggplot() +
aes(.fitted, .resid, color = Tratamiento) +
geom_hline(yintercept = 0) +
geom_point() +
labs(
title = "Residuos en forma de embudo",
subtitle = "La varianza crece con la media! Poisson o BN",
x = "Predichos (escala VR)"
)
qq_fig <-
malezas_fit %>%
ggplot() +
aes(sample = .std.resid) +
geom_qq_line() +
geom_qq()
resid_fig + qq_fig
Termina de descartar normalidad con un test SW sobre los residuos no estandarizados. Aunque acá… no se rechaza H0.
shapiro.test(malezas_fit$.resid)
##
## Shapiro-Wilk normality test
##
## data: malezas_fit$.resid
## W = 0.94037, p-value = 0.2941
malezas_glm <- stats::glm(malezas ~ Tratamiento, data = malezas, family = poisson)
Como en Poisson la varianza crece con la media (se espera sea igual a la media!), los residuos son naturalmente heterocedásticos. Por esto, no se suelen plotear los residuos “crudos” (Y - Y sombrero), sino los Pearson (que son estandarizados).
malezas %>%
mutate(
predichos = predict(malezas_glm, type = "response"),
residuos_crudos = malezas - predichos
) %>%
ggplot() +
aes(predichos, residuos_crudos) +
geom_hline(yintercept = 0) +
geom_segment(aes(xend = predichos, yend = 0), color = "Red") +
geom_point()
Los residuos Pearson están estandarizados! Compararlos con el umbral 3 es que estén a 3 desvíos estándar de 0.
pinos <- read_csv("Clases 13 y 14 Regresión Poisson/pinos.csv")
## Parsed with column specification:
## cols(
## parcela = col_character(),
## n_conos = col_double(),
## altura = col_double(),
## edad = col_double(),
## area_basal = col_double()
## )
pinos_glm <- MASS::glm.nb(n_conos ~ edad * altura, data = pinos)
pinos %>%
rowid_to_column("obs_id") %>%
mutate(
residuos_pearson = resid(pinos_glm, type = "pearson"),
predichos = predict(pinos_glm, type = "response")
) %>%
ggplot() +
aes(x = predichos, y = residuos_pearson, color = abs(residuos_pearson) > 3) +
geom_hline(yintercept = c(0, -3, 3)) +
geom_point() +
geom_text(aes(label = ifelse(abs(residuos_pearson) > 3, obs_id, "")),
position = position_nudge(x = 30)) +
scale_color_manual(values = c("Black", "Red"))
La devianza es un cociente de verosimilitud del modelo de interés vs. el modelo saturado. El modelo saturado es el modelo que tiene un parámetro por observación, de modo que tiene ajuste perfecto (pero obviamente overfitea). La devianza se usa para comparar modelos, no por sí sola.
\[ D = 2 \log \frac{ \mathcal{L}(\text{Modelo saturado}) }{ \mathcal{L}(\text{Modelo reducido}) } = 2 \left( \log\{ P(y|\hat{\theta}_s) \} - \log\{ P(y|\hat{\theta}_0 \} \right) \]
Eso es la devianza de un modelo. Pero luego, calculamos un porcentaje de devianza explicada, de nuestro modelo respecto del modelo NULO (ojo, no el saturado!), que es el que sólo predice con un intercept (y por ende con la media de los datos).
\[ \text{Devianza explicada} = 1 - \frac{D_0}{D_{\text{Nulo}}} \]
devianza_explicada <- 1 - (malezas_glm$deviance / malezas_glm$null.deviance)
100 * devianza_explicada
## [1] 67.67756
A mano:
modelo_nulo <- glm(malezas ~ 1, data = malezas, family = poisson)
malezas_con_id <- malezas %>% rowid_to_column("obs_id")
modelo_saturado <- glm(malezas ~ as.factor(obs_id), data = malezas_con_id, family = poisson)
logLik(modelo_nulo)
## 'log Lik.' -65.40679 (df=1)
logLik(malezas_glm)
## 'log Lik.' -41.88149 (df=3)
str(logLik(modelo_saturado))
## Class 'logLik' : -30.65 (df=18)
devianza <- 2 * (logLik(modelo_saturado) - logLik(malezas_glm))
malezas_glm$deviance
## [1] 22.47112
devianza # Coincide
## 'log Lik.' 22.47112 (df=18)
stats::glm
y tabla de devianza con car::Anova
El test del cociente de verosimilitud (LR por likelihood ratio). El estadístico del test es el siguiente, y asintóticamente se distribuye como chi cuadrado. Si se rechaza H0, significa que el tratamiento tiene efecto.
Este test sólo es válido para modelos anidados. Nótese que el modelo nulo está anidado en cualquier otro modelo con intercept.
\[ H_0: \alpha_i = 0 \\ \Delta D = D_{\text{Nulo}} - D_{\text{modelo}} \approx \chi^2_{ \text{GL}_{\text{modelo}} - \text{GL}_{\text{Nulo}} } \approx \chi^2_{ \text{parametros}_0 - \text{parametros}_{\text{Nulo}} } \]
car::Anova(malezas_glm, test.staistic = "LR")
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
## Analysis of Deviance Table (Type II tests)
##
## Response: malezas
## LR Chisq Df Pr(>Chisq)
## Tratamiento 47.051 2 6.069e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El mismo test (asintótico) de cociente de verosimilitud lo hace drop1
, contra una distro chi cuadrado:
drop1(malezas_glm, test = "Chisq")
## Single term deletions
##
## Model:
## malezas ~ Tratamiento
## Df Deviance AIC LRT Pr(>Chi)
## <none> 22.471 89.763
## Tratamiento 2 69.522 132.814 47.051 6.069e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El objeto glm tiene muchos atributos para usar:
names(malezas_glm)
## [1] "coefficients" "residuals" "fitted.values"
## [4] "effects" "R" "rank"
## [7] "qr" "family" "linear.predictors"
## [10] "deviance" "aic" "null.deviance"
## [13] "iter" "weights" "prior.weights"
## [16] "df.residual" "df.null" "y"
## [19] "converged" "boundary" "model"
## [22] "call" "formula" "terms"
## [25] "data" "offset" "control"
## [28] "method" "contrasts" "xlevels"
summary
Acá se ven los \(\hat{\beta}_j\), parámetros estimados del efecto de los tratamientos, con el intercept representando el tratamiento “base” de los controles.
La solución se encuentra con un algoritmo iterativo, no hay solución cerrada en GLM. No hubo ajuste por sub/sobre dispersión.
summary(malezas_glm)
##
## Call:
## stats::glm(formula = malezas ~ Tratamiento, family = poisson,
## data = malezas)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3094 -0.4524 0.2829 0.7418 1.3588
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.5520 0.1140 22.394 < 2e-16 ***
## Tratamientoherb1 -1.5712 0.2747 -5.719 1.07e-08 ***
## Tratamientoherb2 -0.8781 0.2103 -4.175 2.98e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 69.522 on 17 degrees of freedom
## Residual deviance: 22.471 on 15 degrees of freedom
## AIC: 89.763
##
## Number of Fisher Scoring iterations: 5
Atenti: la magnitud del efecto (multiplicativo) del tratamiento \(j\)-ésimo es \(e^{\beta_j}\)! Acá están ambas escalas.
Notar que el \(\beta_0\) está asociado al control, o a veces al primer nivel del factor si no hay controles. Luego, los \(e^{\beta_j}\) del resto de los tratamientos son las magnitudes por las cuales se multiplican ese efecto de base.
\[ \mathbb{E}(y_{ij}) = e^{\beta_0} \cdot e^{\beta_i} \]
malezas_glm$coefficients
## (Intercept) Tratamientoherb1 Tratamientoherb2
## 2.5520460 -1.5712167 -0.8780695
exp(malezas_glm$coefficients)
## (Intercept) Tratamientoherb1 Tratamientoherb2
## 12.8333333 0.2077922 0.4155844
confint(malezas_glm)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 2.320048 2.7673841
## Tratamientoherb1 -2.145398 -1.0609226
## Tratamientoherb2 -1.303601 -0.4760456
exp(confint(malezas_glm))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 10.1761630 15.9169419
## Tratamientoherb1 0.1170215 0.3461363
## Tratamientoherb2 0.2715520 0.6212351
sjPlot::plot_model devuelve un ggplot y tiene esta $data, se puede usar:
sjPlot::plot_model(malezas_glm)$data
## Waiting for profiling to be done...
## term estimate std.error conf.low conf.high statistic
## 1 Tratamientoherb1 0.2077922 0.274749 0.1170215 0.3461363 -5.718735
## 2 Tratamientoherb2 0.4155844 0.210326 0.2715520 0.6212351 -4.174803
## df.residual p.value p.stars p.label group xpos xmin xmax
## 1 15 1.073200e-08 *** 0.21 *** neg 2 1.825 2.175
## 2 15 2.982438e-05 *** 0.42 *** neg 1 0.825 1.175
tab_model
da una tabla linda para presentación de los estimados, IC y p valores:
sjPlot::tab_model(malezas_glm)
## Waiting for profiling to be done...
malezas | |||
---|---|---|---|
Predictors | Incidence Rate Ratios | CI | p |
(Intercept) | 12.83 | 10.18 – 15.92 | <0.001 |
Tratamiento: herb1 | 0.21 | 0.12 – 0.35 | <0.001 |
Tratamiento: herb2 | 0.42 | 0.27 – 0.62 | <0.001 |
Observations | 18 | ||
R2 Nagelkerke | 0.947 |
sjPlot::plot_model(malezas_glm) +
geom_hline(yintercept = 1, size = 2) +
geom_text(aes(label = round(estimate, 2)), position = position_nudge(x = .25)) +
geom_text(aes(label = p.stars), position = position_nudge(x = -.25))
## Waiting for profiling to be done...
Los residuos Pearson son los usuales: y menos y sombrero, estandarizado. Los de la devianza…?
\[ e_{i, \text{Pearson}} = \frac{Y_i - \hat{Y}_i}{\sqrt{\hat{\sigma}^2_{Y_i}}} = \frac{Y_i - \hat{\mu}_i}{\sqrt{\hat{\mu}^2_i}} \quad \text{para Poisson} \]
\[ e_{i, \text{devianza}} = \text{signo}(Y_i - \hat{Y}_i) \sqrt{d_i} \quad d_i = \dots \]
malezas_tabla <-
tibble(
Trat = rep(c("herb1", "herb2", "control"), each = 6),
Malezas = malezas$malezas,
Predichos_escala_log = predict(malezas_glm),
Predichos_escala_VR = fitted(malezas_glm),
resid_Pearson = resid(malezas_glm, type = "pearson"),
resid_devianza = resid(malezas_glm),
d_cook = cooks.distance(malezas_glm)
) %>%
mutate_if(is.double, round, 2)
malezas_tabla
## # A tibble: 18 x 7
## Trat Malezas Predichos_escal… Predichos_escal… resid_Pearson
## <chr> <int> <dbl> <dbl> <dbl>
## 1 herb1 0 0.98 2.67 -1.63
## 2 herb1 4 0.98 2.67 0.82
## 3 herb1 1 0.98 2.67 -1.02
## 4 herb1 4 0.98 2.67 0.82
## 5 herb1 3 0.98 2.67 0.2
## 6 herb1 4 0.98 2.67 0.82
## 7 herb2 7 1.67 5.33 0.72
## 8 herb2 7 1.67 5.33 0.72
## 9 herb2 6 1.67 5.33 0.290
## 10 herb2 5 1.67 5.33 -0.14
## 11 herb2 1 1.67 5.33 -1.88
## 12 herb2 6 1.67 5.33 0.290
## 13 cont… 6 2.55 12.8 -1.91
## 14 cont… 11 2.55 12.8 -0.51
## 15 cont… 12 2.55 12.8 -0.23
## 16 cont… 14 2.55 12.8 0.33
## 17 cont… 16 2.55 12.8 0.88
## 18 cont… 18 2.55 12.8 1.44
## # … with 2 more variables: resid_devianza <dbl>, d_cook <dbl>
malezas_fit_poisson <- augment(malezas_glm, data = malezas)
Contrastar este ajuste con el del modelo lineal que asumía normalidad! Acá la magnitud de los residuos no crecen con la media, está bien.
fig_resid <-
malezas_tabla %>%
ggplot() +
aes(x = Predichos_escala_VR, y = resid_Pearson, color = Trat) +
geom_hline(yintercept = 0) +
geom_point() +
guides(color = F)
fig_cook <-
malezas_fit_poisson %>%
rowid_to_column("obs_id") %>%
ggplot() +
aes(x = obs_id, y = .cooksd, color = Tratamiento) +
geom_segment(aes(xend = obs_id, yend = 0), size = .25) +
geom_point()
fig_resid + fig_cook
pinos_glm_edad <- MASS::glm.nb(n_conos ~ edad, data = pinos, link = log)
new_edad <- data.frame(edad = min(pinos$edad) : max(pinos$edad))
eta <- predict(pinos_glm_edad, new_edad)
pinos %>%
ggplot() +
aes(x = edad, y = n_conos) +
geom_point() +
geom_line(data = new_edad, aes(x = edad, y = exp(eta)))
La distribución Poisson implica que la varianza y la media son iguales. Se testea que el ajuste no viole este supuesto con el siguiente estadístico:
\[ \hat{\phi} = \frac{ \sum^{n}_{i=1} e^2_{i, \text{Pearson}} }{ \text{g.l. residuos }} \]
Se toleran valores \(< 1.5\) como aceptables para un ajuste Poisson. Valores \(< 1\) indican subdispersión, pero son poco usuales en datasets biológicos. Valores muy por encima de \(1\) indican sobre dispersión, un problema que debe ser abordado con soluciones varias (ver abajo).
calcular_estadistico_de_dispersion <- function(modelo, gl_menos_1 = F) {
residuos_cuadrados <- resid(modelo, type = "pearson")^2
# grados_libertad <- modelo$df.residual
n <- length(modelo$residuals)
p <- length(modelo$coefficients)
grados_libertad_residuos <- n - p
if (gl_menos_1) { grados_libertad_residuos <- grados_libertad_residuos - 1 }
return (sum(residuos_cuadrados/grados_libertad_residuos))
}
calcular_estadistico_de_dispersion(malezas_glm)
## [1] 1.161472
emmeans
(ambas escalas)El paquete emmeans
(Estimated Marginal Means) permite calcular y plotear contrastes entre medias de los tratamientos, con sus IC, en ambas escalas. Usa IC de Tukey por default.
# ajust = "tukey"
contrastes <- emmeans::emmeans(malezas_glm, pairwise ~ Tratamiento)
contrastes_VR <- emmeans::emmeans(malezas_glm, pairwise ~ Tratamiento, type = "response")
summary(contrastes)
## $emmeans
## Tratamiento emmean SE df asymp.LCL asymp.UCL
## control 2.552 0.114 Inf 2.329 2.78
## herb1 0.981 0.250 Inf 0.491 1.47
## herb2 1.674 0.177 Inf 1.328 2.02
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## control - herb1 1.571 0.275 Inf 5.719 <.0001
## control - herb2 0.878 0.210 Inf 4.175 0.0001
## herb1 - herb2 -0.693 0.306 Inf -2.264 0.0610
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
summary(contrastes_VR)
## $emmeans
## Tratamiento rate SE df asymp.LCL asymp.UCL
## control 12.83 1.462 Inf 10.26 16.05
## herb1 2.67 0.667 Inf 1.63 4.35
## herb2 5.33 0.943 Inf 3.77 7.54
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df z.ratio p.value
## control / herb1 4.81 1.322 Inf 5.719 <.0001
## control / herb2 2.41 0.506 Inf 4.175 0.0001
## herb1 / herb2 0.50 0.153 Inf -2.264 0.0610
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scale
También acepta un parámetro offset:
emmeans::emmeans(malezas_glm, pairwise ~ Tratamiento, offset = 0)$emmeans
## Tratamiento emmean SE df asymp.LCL asymp.UCL
## control 2.552 0.114 Inf 2.329 2.78
## herb1 0.981 0.250 Inf 0.491 1.47
## herb2 1.674 0.177 Inf 1.328 2.02
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
plot(contrastes, comparison = T) + plot(contrastes_VR, comparison = T)
confint(contrastes)
## $emmeans
## Tratamiento emmean SE df asymp.LCL asymp.UCL
## control 2.552 0.114 Inf 2.329 2.78
## herb1 0.981 0.250 Inf 0.491 1.47
## herb2 1.674 0.177 Inf 1.328 2.02
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df asymp.LCL asymp.UCL
## control - herb1 1.571 0.275 Inf 0.927 2.2151
## control - herb2 0.878 0.210 Inf 0.385 1.3710
## herb1 - herb2 -0.693 0.306 Inf -1.411 0.0245
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
confint(contrastes_VR)
## $emmeans
## Tratamiento rate SE df asymp.LCL asymp.UCL
## control 12.83 1.462 Inf 10.26 16.05
## herb1 2.67 0.667 Inf 1.63 4.35
## herb2 5.33 0.943 Inf 3.77 7.54
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df asymp.LCL asymp.UCL
## control / herb1 4.81 1.322 Inf 2.528 9.16
## control / herb2 2.41 0.506 Inf 1.470 3.94
## herb1 / herb2 0.50 0.153 Inf 0.244 1.02
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
## Intervals are back-transformed from the log scale
malezas_medias <- contrastes_VR$emmeans %>% as_tibble
malezas_medias
## # A tibble: 3 x 6
## Tratamiento rate SE df asymp.LCL asymp.UCL
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 control 12.8 1.46 Inf 10.3 16.0
## 2 herb1 2.67 0.667 Inf 1.63 4.35
## 3 herb2 5.33 0.943 Inf 3.77 7.54
malezas_medias %>%
ggplot() +
aes(x = Tratamiento, y = rate) +
geom_errorbar(
aes(ymin = rate - SE, ymax = rate + SE), width = .1
) +
geom_point(colour = "blue", size = 1) +
labs(y = "Malezas/m2") +
geom_text(
label = c("A", "B", "B"),
position = position_nudge(x = .2)
) +
geom_text(
mapping = aes(label = round(rate, 1)),
position = position_nudge(x = -.2),
size = 3.5
)
Cuando hay sobredispersión, en lugar de Poisson conviene la BN.
roadkills_full <- read_tsv("tps/TP 7/roadkills.txt")
## Parsed with column specification:
## cols(
## .default = col_double()
## )
## See spec(...) for full column specifications.
# OPEN.L + MONT.S + SQ.POLIC + D.PARK + SQ.SHRUB + SQ.WAT.RES + L.WAT.C + SQ.L.P.ROAD + SQ.D.WAT.COUR
roadkills <- select(roadkills_full, TOT.N, OPEN.L, MONT.S, POLIC, D.PARK, SHRUB, WAT.RES, L.WAT.C, L.P.ROAD, D.WAT.COUR)
roadkills
## # A tibble: 52 x 10
## TOT.N OPEN.L MONT.S POLIC D.PARK SHRUB WAT.RES L.WAT.C L.P.ROAD
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 22 22.7 0 4.81 250. 0.406 0.043 0.583 1.98
## 2 14 24.7 0 2.22 741. 0.735 0.182 1.42 1.76
## 3 65 30.1 0.258 1.95 1240. 0.474 0.453 2.00 1.25
## 4 55 50.3 1.78 0.625 1740. 0.607 0.026 1.92 0.666
## 5 88 43.6 2.43 0.791 2232. 0.173 0 2.17 0.653
## 6 104 31.4 0 0.054 2724. 0.325 0.039 2.39 1.31
## 7 49 24.8 0 0.022 3216. 0.055 0.114 1.16 0.685
## 8 66 56.2 0 11.3 3709. 0.092 0.224 2.43 0.677
## 9 26 48.7 1.11 1.24 4206. 1.74 0.177 2.42 0.664
## 10 47 15.6 0 0.119 4704. 0 0 0.211 0.654
## # … with 42 more rows, and 1 more variable: D.WAT.COUR <dbl>
ggpairs para AED de muchas covariables
sqrt
o log
Si las covariables tienen valores extremos que “distorsionan” la escala, se pueden modificar con sqrt
o más fuertemente con log
:
fig_polic <- roadkills %>% ggplot(aes(y = log(TOT.N), x = POLIC)) + geom_point()
fig_polic_sq <- roadkills %>% ggplot(aes(y = log(TOT.N), x = sqrt(POLIC))) + geom_point()
fig_polic_log <- roadkills %>% ggplot(aes(y = log(TOT.N), x = log(POLIC))) + geom_point()
fig_polic + fig_polic_sq + fig_polic_log
roadkills$SQ.POLIC <- sqrt(roadkills$POLIC)
roadkills$SQ.SHRUB <- sqrt(roadkills$SHRUB)
roadkills$SQ.WAT.RES <- sqrt(roadkills$WAT.RES)
roadkills$SQ.L.P.ROAD <- sqrt(roadkills$L.P.ROAD)
roadkills$SQ.D.WAT.COUR <- sqrt(roadkills$D.WAT.COUR)
roadkills %<>% select(-POLIC, -SHRUB, -WAT.RES, -L.P.ROAD, -D.WAT.COUR)
road_glm <- glm(TOT.N ~ ., family = poisson, data = roadkills)
calcular_estadistico_de_dispersion(road_glm)
## [1] 5.928003
AER::dispersiontest(road_glm)
## Registered S3 method overwritten by 'AER':
## method from
## nobs.survreg insight
##
## Overdispersion test
##
## data: road_glm
## z = 4.0738, p-value = 2.313e-05
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 4.797236
# AER::dispersiontest(road_glm, trafo = 1) # 2 para NB, es mu^2
# AER::dispersiontest(road_glm, alternative = "less")
Hay sobredispersión! (estadístico > 1.5, criterio de Zuur). Se puede probar con:
Podemos implementar todos los modelos y compararlos con AIC
para elegir uno de ellos. También OLRE para sobredispersión.
Se reescala el error estándar de los betas con la raíz del parámetro de dispersión. Los estimados no cambian, pero sí sus IC.
\[ \sigma^*_{\beta_j} = \sigma_{\beta_j} \cdot \sqrt{\hat{\phi}} \]
road_glm_quasi <- glm(TOT.N ~ ., family = quasipoisson, data = roadkills)
summary(road_glm_quasi)
##
## Call:
## glm(formula = TOT.N ~ ., family = quasipoisson, data = roadkills)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.8398 -1.3965 -0.1409 1.4641 4.3749
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.749e+00 3.814e-01 9.830 1.86e-12 ***
## OPEN.L -3.025e-03 3.847e-03 -0.786 0.43604
## MONT.S 8.697e-02 3.309e-02 2.628 0.01194 *
## D.PARK -1.301e-04 1.445e-05 -9.004 2.33e-11 ***
## L.WAT.C 3.355e-01 1.005e-01 3.338 0.00177 **
## SQ.POLIC -1.787e-01 1.139e-01 -1.570 0.12400
## SQ.SHRUB -6.112e-01 2.863e-01 -2.135 0.03867 *
## SQ.WAT.RES 2.243e-01 1.717e-01 1.306 0.19851
## SQ.L.P.ROAD 4.517e-01 3.282e-01 1.376 0.17597
## SQ.D.WAT.COUR 7.355e-03 1.188e-02 0.619 0.53910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.928003)
##
## Null deviance: 1071.44 on 51 degrees of freedom
## Residual deviance: 270.23 on 42 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
Vamos a agregar de a una variable por vez, hasta que no sean significativas las nuevas variables agregadas.
m0 <- MASS::glm.nb(TOT.N ~ 1, data = roadkills, link = log)
m0
##
## Call: MASS::glm.nb(formula = TOT.N ~ 1, data = roadkills, link = log,
## init.theta = 1.183705161)
##
## Coefficients:
## (Intercept)
## 3.254
##
## Degrees of Freedom: 51 Total (i.e. Null); 51 Residual
## Null Deviance: 57.57
## Residual Deviance: 57.57 AIC: 447.7
f <- as.formula("TOT.N ~ OPEN.L + MONT.S + D.PARK + L.WAT.C + SQ.POLIC + SQ.SHRUB + SQ.WAT.RES + SQ.L.P.ROAD + SQ.D.WAT.COUR")
add1(m0, f, test = "Chisq")
## Single term additions
##
## Model:
## TOT.N ~ 1
## Df Deviance AIC LRT Pr(>Chi)
## <none> 57.569 445.66
## OPEN.L 1 52.060 442.16 5.509 0.01892 *
## MONT.S 1 57.424 447.52 0.145 0.70289
## D.PARK 1 21.310 411.41 36.259 1.727e-09 ***
## L.WAT.C 1 57.538 447.63 0.031 0.86039
## SQ.POLIC 1 55.076 445.17 2.493 0.11432
## SQ.SHRUB 1 54.427 444.52 3.142 0.07629 .
## SQ.WAT.RES 1 57.283 447.38 0.286 0.59267
## SQ.L.P.ROAD 1 56.750 446.85 0.819 0.36555
## SQ.D.WAT.COUR 1 57.562 447.66 0.007 0.93230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- update(m0, . ~ . + D.PARK)
add1(m1, f, test = "Chisq")
## Single term additions
##
## Model:
## TOT.N ~ D.PARK
## Df Deviance AIC LRT Pr(>Chi)
## <none> 54.742 391.09
## OPEN.L 1 47.315 385.67 7.4273 0.006424 **
## MONT.S 1 51.604 389.96 3.1378 0.076499 .
## L.WAT.C 1 52.399 390.75 2.3434 0.125818
## SQ.POLIC 1 54.614 392.97 0.1277 0.720860
## SQ.SHRUB 1 54.712 393.06 0.0300 0.862535
## SQ.WAT.RES 1 54.525 392.88 0.2168 0.641481
## SQ.L.P.ROAD 1 52.056 390.41 2.6858 0.101245
## SQ.D.WAT.COUR 1 54.593 392.94 0.1486 0.699856
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- update(m1, . ~ . + OPEN.L)
add1(m2, f, test = "Chisq")
## Single term additions
##
## Model:
## TOT.N ~ D.PARK + OPEN.L
## Df Deviance AIC LRT Pr(>Chi)
## <none> 51.839 385.43
## MONT.S 1 49.921 385.51 1.9183 0.16604
## L.WAT.C 1 46.954 382.55 4.8848 0.02709 *
## SQ.POLIC 1 51.577 387.17 0.2620 0.60878
## SQ.SHRUB 1 51.514 387.11 0.3248 0.56871
## SQ.WAT.RES 1 51.738 387.33 0.1013 0.75029
## SQ.L.P.ROAD 1 47.961 383.55 3.8784 0.04891 *
## SQ.D.WAT.COUR 1 50.966 386.56 0.8737 0.34993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- update(m2, . ~ . + L.WAT.C)
add1(m3, f, test = "Chisq")
## Single term additions
##
## Model:
## TOT.N ~ D.PARK + OPEN.L + L.WAT.C
## Df Deviance AIC LRT Pr(>Chi)
## <none> 52.003 382.25
## MONT.S 1 50.448 382.70 1.55461 0.2125
## SQ.POLIC 1 51.999 384.25 0.00388 0.9504
## SQ.SHRUB 1 51.954 384.20 0.04914 0.8246
## SQ.WAT.RES 1 51.988 384.24 0.01474 0.9034
## SQ.L.P.ROAD 1 49.336 381.59 2.66655 0.1025
## SQ.D.WAT.COUR 1 51.999 384.25 0.00382 0.9507
Cuando las variables que quedan no son significativas, dejo de agrandar el modelo.
AIC(m0, m1, m2, m3)
## df AIC
## m0 2 447.6643
## m1 3 393.0929
## m2 4 387.4317
## m3 5 384.2536
anova(m0, m1, m2, m3)
## Likelihood ratio tests of Negative Binomial Models
##
## Response: TOT.N
## Model theta Resid. df 2 x log-lik. Test
## 1 1 1.183705 51 -443.6643
## 2 D.PARK 3.681040 50 -387.0929 1 vs 2
## 3 D.PARK + OPEN.L 4.132845 49 -379.4317 2 vs 3
## 4 D.PARK + OPEN.L + L.WAT.C 4.725770 48 -374.2536 3 vs 4
## df LR stat. Pr(Chi)
## 1
## 2 1 56.571344 5.417888e-14
## 3 1 7.661245 5.641956e-03
## 4 1 5.178081 2.287358e-02
Lo mismo puede hacerse con stepAIC
, step
lmerTest::step(m3)
## Start: AIC=382.25
## TOT.N ~ D.PARK + OPEN.L + L.WAT.C
##
## Df Deviance AIC
## <none> 52.003 382.25
## - L.WAT.C 1 57.501 385.75
## - OPEN.L 1 63.077 391.33
## - D.PARK 1 168.593 496.84
##
## Call: MASS::glm.nb(formula = TOT.N ~ D.PARK + OPEN.L + L.WAT.C, data = roadkills,
## init.theta = 4.725770367, link = log)
##
## Coefficients:
## (Intercept) D.PARK OPEN.L L.WAT.C
## 4.4666276 -0.0001158 -0.0107139 0.1867415
##
## Degrees of Freedom: 51 Total (i.e. Null); 48 Residual
## Null Deviance: 189.7
## Residual Deviance: 52 AIC: 384.3
MASS::stepAIC(m3)
## Start: AIC=382.25
## TOT.N ~ D.PARK + OPEN.L + L.WAT.C
##
## Df AIC
## <none> 382.25
## - L.WAT.C 1 385.43
## - OPEN.L 1 390.70
## - D.PARK 1 442.70
##
## Call: MASS::glm.nb(formula = TOT.N ~ D.PARK + OPEN.L + L.WAT.C, data = roadkills,
## init.theta = 4.725770367, link = log)
##
## Coefficients:
## (Intercept) D.PARK OPEN.L L.WAT.C
## 4.4666276 -0.0001158 -0.0107139 0.1867415
##
## Degrees of Freedom: 51 Total (i.e. Null); 48 Residual
## Null Deviance: 189.7
## Residual Deviance: 52 AIC: 384.3
Se resta 1 a los grados de libertad porque “se estimó k” (el parámetro extra de la BN).
road_glm_bn <- m3
calcular_estadistico_de_dispersion(road_glm_bn, gl_menos_1 = T)
## [1] 1.04612
\[ VIF_j = \frac{1}{1 - R^2_j} \]
donde R^2_j se obtiene al hacer una regresión lineal múltiple con Xj como VR y las demás covariables como VE. El VIF va de 1 a infinito. Tomamos valores superiores a 5 como indicación de colinealidad importante. Ver diapo de clase 8, p. 12. Se resuelve eliminando variables o combinando VE asociadas en una nueva variable.
car::vif(road_glm_bn)
## D.PARK OPEN.L L.WAT.C
## 1.105490 1.127361 1.130730
No hay demasiada colinealidad.
Hacemos lo usual para un ajuste ahora: residuos vs. predichos, summary, confint, exp(confint()).
summary(road_glm_bn)
##
## Call:
## MASS::glm.nb(formula = TOT.N ~ D.PARK + OPEN.L + L.WAT.C, data = roadkills,
## init.theta = 4.725770367, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7799 -0.7976 -0.1921 0.4513 2.0395
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.467e+00 1.716e-01 26.022 < 2e-16 ***
## D.PARK -1.158e-04 1.082e-05 -10.708 < 2e-16 ***
## OPEN.L -1.071e-02 3.146e-03 -3.405 0.000661 ***
## L.WAT.C 1.867e-01 7.911e-02 2.361 0.018246 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.7258) family taken to be 1)
##
## Null deviance: 189.709 on 51 degrees of freedom
## Residual deviance: 52.003 on 48 degrees of freedom
## AIC: 384.25
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.73
## Std. Err.: 1.16
##
## 2 x log-likelihood: -374.254
exp(confint(road_glm_bn))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 60.9673187 125.8304307
## D.PARK 0.9998625 0.9999055
## OPEN.L 0.9830986 0.9956006
## L.WAT.C 1.0310932 1.4103815
car::Anova(road_glm_bn)
## Analysis of Deviance Table (Type II tests)
##
## Response: TOT.N
## LR Chisq Df Pr(>Chisq)
## D.PARK 116.590 1 < 2.2e-16 ***
## OPEN.L 11.074 1 0.0008754 ***
## L.WAT.C 5.498 1 0.0190389 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nuevos_datos <-
roadkills %$%
tibble(
D.PARK = seq(min(D.PARK), max(D.PARK), length.out = nrow(roadkills)),
OPEN.L = mean(OPEN.L),
L.WAT.C = mean(L.WAT.C)
)
nuevos_datos %<>%
mutate(
TOT.N = predict(road_glm_bn, newdata = nuevos_datos, type = "response")
)
nuevos_datos %>%
ggplot() +
aes(x = D.PARK, y = TOT.N) +
geom_point(data = roadkills, mapping = aes(color = "Datos")) +
geom_line(aes(color = "Predichos"))
aranias <- read_tsv("tps/TP 7/arania.txt") %>% mutate(trat = factor(trat))
## Parsed with column specification:
## cols(
## trat = col_character(),
## n_huevos = col_double()
## )
aranias
## # A tibble: 30 x 2
## trat n_huevos
## <fct> <dbl>
## 1 control 277
## 2 control 250
## 3 control 231
## 4 control 253
## 5 control 258
## 6 control 246
## 7 control 243
## 8 control 266
## 9 control 248
## 10 control 234
## # … with 20 more rows
aranias_pois <- stats::glm(n_huevos ~ trat, data = aranias, family = poisson)
calcular_estadistico_de_dispersion(aranias_pois)
## [1] 0.6218766
Hay subdispersión, no está bien el supuesto de Poisson. Usamos Conwell Poisson:
aranias_compois <-
glmmTMB::glmmTMB(n_huevos ~ trat, data = aranias, family = glmmTMB::compois)
calcular_estadistico_de_dispersion(aranias_compois)
## [1] Inf
aranias_compois
## Formula: n_huevos ~ trat
## Data: aranias
## AIC BIC logLik df.resid
## 235.0349 240.6397 -113.5174 26
##
## Number of obs: 30
##
## Overdispersion parameter for compois family (): 0.554
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) tratglifo alta tratglifo baja
## 5.5239 -0.3928 -0.2206
$$ _i = ( ) = _0 + _1 Ancho_i \
odds_i = e^{_i} = e^{_0} e^{_1 Ancho_i} $$
abejas <- read_csv("Clases 15 y 16 Regresión Bernoulli y binomial/abejas.csv")
## Parsed with column specification:
## cols(
## parasitos = col_double(),
## ancho = col_double()
## )
abejas
## # A tibble: 917 x 2
## parasitos ancho
## <dbl> <dbl>
## 1 1 5.16
## 2 0 5.23
## 3 0 5
## 4 1 5.72
## 5 1 5.5
## 6 0 5.63
## 7 1 5.04
## 8 1 5.89
## 9 0 5.12
## 10 1 5.15
## # … with 907 more rows
table(abejas$parasitos)
##
## 0 1
## 407 510
summary(abejas)
## parasitos ancho
## Min. :0.0000 Min. :3.52
## 1st Qu.:0.0000 1st Qu.:4.92
## Median :1.0000 Median :5.18
## Mean :0.5562 Mean :5.19
## 3rd Qu.:1.0000 3rd Qu.:5.48
## Max. :1.0000 Max. :6.64
ggplot(abejas) +
aes(x = ancho) +
geom_histogram(binwidth = .25)
base_fig <-
ggplot(abejas) +
aes(ancho, parasitos) +
geom_jitter(height = .05, alpha = .5)
base_fig +
geom_smooth(method = "lm", se = F) # Para que esto fucnione, la VR NO tiene que se factor!
base_fig +
geom_smooth(method = "glm", method.args = list(family = "binomial"), se = F)
Para la regresión logística o Bernoulli, usamos family = binomial con glm
:
abejas_m1 <- glm(parasitos ~ ancho, data = abejas, family = binomial)
# En el R de la clase lo hizo sin Chisq:
drop1(abejas_m1, test = "Chisq")
## Single term deletions
##
## Model:
## parasitos ~ ancho
## Df Deviance AIC LRT Pr(>Chi)
## <none> 1104.9 1108.9
## ancho 1 1259.6 1261.6 154.73 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coefficients(abejas_m1)
## (Intercept) ancho
## -11.244759 2.217536
exp(coefficients(abejas_m1))
## (Intercept) ancho
## 1.307565e-05 9.184671e+00
exp(confint(abejas_m1))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 1.577729e-06 9.780426e-05
## ancho 6.225532e+00 1.382548e+01
predicciones_logit <- predict(abejas_m1, type = "link")
predicciones_odds <- exp(predicciones_logit)
predicciones_prob <- predict(abejas_m1, type = "response")
abejas_ajuste <- as.data.frame(cbind(
"ancho" = abejas$ancho,
"parasitos" = abejas$parasitos,
"pred_logit" = round(predicciones_logit, 2),
"pred_odds" = round(predicciones_odds, 2),
"pred_prob" = round(predicciones_prob, 2),
"residuos_pearson" = residuals(abejas_m1, type = "pearson")
))
abejas_ajuste
## ancho parasitos pred_logit pred_odds pred_prob residuos_pearson
## 1 5.16 1 0.20 1.22 0.55 0.9058664
## 2 5.23 0 0.35 1.42 0.59 -1.1930072
## 3 5.00 0 -0.16 0.85 0.46 -0.9244656
## 4 5.72 1 1.44 4.22 0.81 0.4868625
## 5 5.50 1 0.95 2.59 0.72 0.6213601
## 6 5.63 0 1.24 3.46 0.78 -1.8588991
## 7 5.04 1 -0.07 0.93 0.48 1.0347799
## 8 5.89 1 1.82 6.15 0.86 0.4032236
## 9 5.12 0 0.11 1.12 0.53 -1.0560259
## 10 5.15 1 0.18 1.19 0.54 0.9159662
## 11 4.57 0 -1.11 0.33 0.25 -0.5738945
## 12 4.71 0 -0.80 0.45 0.31 -0.6702650
## 13 5.40 0 0.73 2.07 0.67 -1.4404676
## 14 5.35 1 0.62 1.86 0.65 0.7337922
## 15 5.00 0 -0.16 0.85 0.46 -0.9244656
## 16 5.35 1 0.62 1.86 0.65 0.7337922
## 17 5.49 1 0.93 2.53 0.72 0.6282878
## 18 5.35 1 0.62 1.86 0.65 0.7337922
## 19 4.96 0 -0.25 0.78 0.44 -0.8843608
## 20 4.90 0 -0.38 0.68 0.41 -0.8274420
## 21 6.08 1 2.24 9.37 0.90 0.3266291
## 22 4.77 1 -0.67 0.51 0.34 1.3959233
## 23 5.61 1 1.20 3.31 0.77 0.5500154
## 24 5.05 0 -0.05 0.95 0.49 -0.9771638
## 25 5.07 0 0.00 1.00 0.50 -0.9990748
## 26 5.39 1 0.71 2.03 0.67 0.7019591
## 27 4.50 1 -1.27 0.28 0.22 1.8831076
## 28 5.40 1 0.73 2.07 0.67 0.6942190
## 29 5.35 1 0.62 1.86 0.65 0.7337922
## 30 5.00 0 -0.16 0.85 0.46 -0.9244656
## 31 5.19 1 0.26 1.30 0.57 0.8762301
## 32 4.66 0 -0.91 0.40 0.29 -0.6341178
## 33 5.40 1 0.73 2.07 0.67 0.6942190
## 34 5.23 1 0.35 1.42 0.59 0.8382179
## 35 4.76 1 -0.69 0.50 0.33 1.4114870
## 36 5.06 1 -0.02 0.98 0.49 1.0120858
## 37 5.14 1 0.15 1.17 0.54 0.9261787
## 38 5.39 1 0.71 2.03 0.67 0.7019591
## 39 5.58 0 1.13 3.09 0.76 -1.7586492
## 40 5.21 0 0.31 1.36 0.58 -1.1668431
## 41 5.02 0 -0.11 0.89 0.47 -0.9451950
## 42 5.15 1 0.18 1.19 0.54 0.9159662
## 43 5.35 0 0.62 1.86 0.65 -1.3627836
## 44 5.57 1 1.11 3.03 0.75 0.5749580
## 45 4.56 0 -1.13 0.32 0.24 -0.5675665
## 46 5.18 0 0.24 1.27 0.56 -1.1286687
## 47 5.15 0 0.18 1.19 0.54 -1.0917433
## 48 5.37 1 0.66 1.94 0.66 0.7176992
## 49 5.50 1 0.95 2.59 0.72 0.6213601
## 50 4.98 1 -0.20 0.82 0.45 1.1059612
## 51 6.64 1 3.48 32.45 0.97 0.1755485
## 52 5.11 0 0.09 1.09 0.52 -1.0443817
## 53 4.98 1 -0.20 0.82 0.45 1.1059612
## 54 4.16 0 -2.02 0.13 0.12 -0.3642538
## 55 5.63 0 1.24 3.46 0.78 -1.8588991
## 56 5.55 1 1.06 2.89 0.74 0.5878503
## 57 4.93 0 -0.31 0.73 0.42 -0.8554281
## 58 5.71 1 1.42 4.13 0.80 0.4922907
## 59 6.37 1 2.88 17.83 0.95 0.2368158
## 60 5.24 0 0.38 1.46 0.59 -1.2063085
## 61 5.58 1 1.13 3.09 0.76 0.5686182
## 62 5.05 0 -0.05 0.95 0.49 -0.9771638
## 63 5.08 0 0.02 1.02 0.51 -1.0102138
## 64 5.39 1 0.71 2.03 0.67 0.7019591
## 65 5.36 1 0.64 1.90 0.66 0.7257011
## 66 6.39 0 2.93 18.64 0.95 -4.3173771
## 67 5.36 0 0.64 1.90 0.66 -1.3779778
## 68 5.91 1 1.86 6.43 0.87 0.3943804
## 69 4.24 0 -1.84 0.16 0.14 -0.3980399
## 70 4.72 0 -0.78 0.46 0.31 -0.6777380
## 71 4.90 0 -0.38 0.68 0.41 -0.8274420
## 72 5.41 0 0.75 2.12 0.68 -1.4565279
## 73 5.00 0 -0.16 0.85 0.46 -0.9244656
## 74 5.67 1 1.33 3.78 0.79 0.5146156
## 75 4.90 0 -0.38 0.68 0.41 -0.8274420
## 76 5.54 1 1.04 2.83 0.74 0.5944045
## 77 5.52 1 1.00 2.71 0.73 0.6077328
## 78 5.55 1 1.06 2.89 0.74 0.5878503
## 79 5.96 1 1.97 7.18 0.88 0.3731116
## 80 5.55 1 1.06 2.89 0.74 0.5878503
## 81 5.02 0 -0.11 0.89 0.47 -0.9451950
## 82 5.64 1 1.26 3.53 0.78 0.5320211
## 83 5.07 1 0.00 1.00 0.50 1.0009261
## 84 4.82 1 -0.56 0.57 0.36 1.3206416
## 85 5.18 1 0.24 1.27 0.56 0.8859996
## 86 5.07 1 0.00 1.00 0.50 1.0009261
## 87 5.64 1 1.26 3.53 0.78 0.5320211
## 88 5.47 1 0.89 2.42 0.71 0.6423760
## 89 5.57 1 1.11 3.03 0.75 0.5749580
## 90 4.07 0 -2.22 0.11 0.10 -0.3296599
## 91 5.13 1 0.13 1.14 0.53 0.9365050
## 92 4.69 0 -0.84 0.43 0.30 -0.6555652
## 93 4.71 1 -0.80 0.45 0.31 1.4919472
## 94 4.96 0 -0.25 0.78 0.44 -0.8843608
## 95 5.25 0 0.40 1.49 0.60 -1.2197581
## 96 4.69 0 -0.84 0.43 0.30 -0.6555652
## 97 5.05 0 -0.05 0.95 0.49 -0.9771638
## 98 5.42 1 0.77 2.17 0.68 0.6789939
## 99 5.40 1 0.73 2.07 0.67 0.6942190
## 100 4.86 0 -0.47 0.63 0.39 -0.7915463
## 101 5.40 0 0.73 2.07 0.67 -1.4404676
## 102 5.40 1 0.73 2.07 0.67 0.6942190
## 103 5.39 1 0.71 2.03 0.67 0.7019591
## 104 4.66 0 -0.91 0.40 0.29 -0.6341178
## 105 5.12 0 0.11 1.12 0.53 -1.0560259
## 106 4.55 0 -1.15 0.32 0.24 -0.5613083
## 107 5.09 1 0.04 1.04 0.51 0.9789745
## 108 4.96 1 -0.25 0.78 0.44 1.1307602
## 109 5.08 0 0.02 1.02 0.51 -1.0102138
## 110 5.33 1 0.57 1.78 0.64 0.7502461
## 111 4.93 1 -0.31 0.73 0.42 1.1690053
## 112 5.61 1 1.20 3.31 0.77 0.5500154
## 113 5.10 1 0.06 1.07 0.52 0.9681799
## 114 4.85 0 -0.49 0.61 0.38 -0.7828183
## 115 5.05 0 -0.05 0.95 0.49 -0.9771638
## 116 4.94 1 -0.29 0.75 0.43 1.1561153
## 117 4.46 0 -1.35 0.26 0.21 -0.5079999
## 118 5.04 1 -0.07 0.93 0.48 1.0347799
## 119 5.21 1 0.31 1.36 0.58 0.8570133
## 120 5.15 1 0.18 1.19 0.54 0.9159662
## 121 5.14 1 0.15 1.17 0.54 0.9261787
## 122 5.96 1 1.97 7.18 0.88 0.3731116
## 123 5.77 1 1.55 4.71 0.82 0.4606062
## 124 4.79 1 -0.62 0.54 0.35 1.3653089
## 125 4.96 0 -0.25 0.78 0.44 -0.8843608
## 126 4.61 0 -1.02 0.36 0.26 -0.5999200
## 127 4.82 0 -0.56 0.57 0.36 -0.7572077
## 128 5.25 1 0.40 1.49 0.60 0.8198347
## 129 5.07 1 0.00 1.00 0.50 1.0009261
## 130 5.18 1 0.24 1.27 0.56 0.8859996
## 131 4.82 0 -0.56 0.57 0.36 -0.7572077
## 132 5.26 0 0.42 1.52 0.60 -1.2333577
## 133 5.61 1 1.20 3.31 0.77 0.5500154
## 134 5.43 0 0.80 2.22 0.69 -1.4891877
## 135 5.69 0 1.37 3.95 0.80 -1.9867706
## 136 5.50 0 0.95 2.59 0.72 -1.6093728
## 137 5.26 0 0.42 1.52 0.60 -1.2333577
## 138 5.18 1 0.24 1.27 0.56 0.8859996
## 139 5.55 0 1.06 2.89 0.74 -1.7011134
## 140 5.26 1 0.42 1.52 0.60 0.8107948
## 141 5.57 1 1.11 3.03 0.75 0.5749580
## 142 5.36 1 0.64 1.90 0.66 0.7257011
## 143 5.35 1 0.62 1.86 0.65 0.7337922
## 144 5.43 0 0.80 2.22 0.69 -1.4891877
## 145 5.93 1 1.91 6.72 0.87 0.3857311
## 146 5.11 1 0.09 1.09 0.52 0.9575043
## 147 4.77 0 -0.67 0.51 0.34 -0.7163717
## 148 5.72 1 1.44 4.22 0.81 0.4868625
## 149 4.69 0 -0.84 0.43 0.30 -0.6555652
## 150 4.69 0 -0.84 0.43 0.30 -0.6555652
## 151 5.47 0 0.89 2.42 0.71 -1.5567208
## 152 5.20 1 0.29 1.33 0.57 0.8665685
## 153 5.57 1 1.11 3.03 0.75 0.5749580
## 154 5.20 1 0.29 1.33 0.57 0.8665685
## 155 5.05 0 -0.05 0.95 0.49 -0.9771638
## 156 5.38 1 0.69 1.98 0.66 0.7097855
## 157 5.91 1 1.86 6.43 0.87 0.3943804
## 158 4.82 1 -0.56 0.57 0.36 1.3206416
## 159 6.00 1 2.06 7.85 0.89 0.3569254
## 160 5.43 0 0.80 2.22 0.69 -1.4891877
## 161 5.75 0 1.51 4.51 0.82 -2.1234382
## 162 5.46 1 0.86 2.37 0.70 0.6495381
## 163 5.66 1 1.31 3.69 0.79 0.5203532
## 164 5.11 0 0.09 1.09 0.52 -1.0443817
## 165 5.40 0 0.73 2.07 0.67 -1.4404676
## 166 4.94 1 -0.29 0.75 0.43 1.1561153
## 167 5.36 1 0.64 1.90 0.66 0.7257011
## 168 5.50 0 0.95 2.59 0.72 -1.6093728
## 169 5.26 1 0.42 1.52 0.60 0.8107948
## 170 4.94 0 -0.29 0.75 0.43 -0.8649656
## 171 5.55 0 1.06 2.89 0.74 -1.7011134
## 172 5.40 0 0.73 2.07 0.67 -1.4404676
## 173 5.58 1 1.13 3.09 0.76 0.5686182
## 174 6.06 1 2.19 8.97 0.90 0.3339531
## 175 5.36 0 0.64 1.90 0.66 -1.3779778
## 176 5.11 0 0.09 1.09 0.52 -1.0443817
## 177 5.05 0 -0.05 0.95 0.49 -0.9771638
## 178 5.42 1 0.77 2.17 0.68 0.6789939
## 179 5.29 1 0.49 1.63 0.62 0.7842689
## 180 6.11 0 2.30 10.02 0.91 -3.1651265
## 181 5.22 0 0.33 1.39 0.58 -1.1798526
## 182 5.33 0 0.57 1.78 0.64 -1.3328960
## 183 5.40 1 0.73 2.07 0.67 0.6942190
## 184 5.46 1 0.86 2.37 0.70 0.6495381
## 185 4.74 0 -0.73 0.48 0.32 -0.6929350
## 186 4.90 1 -0.38 0.68 0.41 1.2085439
## 187 5.00 0 -0.16 0.85 0.46 -0.9244656
## 188 5.69 0 1.37 3.95 0.80 -1.9867706
## 189 4.96 0 -0.25 0.78 0.44 -0.8843608
## 190 5.17 1 0.22 1.25 0.55 0.8958779
## 191 5.16 1 0.20 1.22 0.55 0.9058664
## 192 5.00 0 -0.16 0.85 0.46 -0.9244656
## 193 5.46 1 0.86 2.37 0.70 0.6495381
## 194 4.92 0 -0.33 0.72 0.42 -0.8459958
## 195 4.57 0 -1.11 0.33 0.25 -0.5738945
## 196 4.85 0 -0.49 0.61 0.38 -0.7828183
## 197 5.20 0 0.29 1.33 0.57 -1.1539769
## 198 5.43 1 0.80 2.22 0.69 0.6715070
## 199 5.65 0 1.28 3.61 0.78 -1.9005813
## 200 4.99 1 -0.18 0.84 0.46 1.0937663
## 201 5.38 1 0.69 1.98 0.66 0.7097855
## 202 4.75 0 -0.71 0.49 0.33 -0.7006608
## 203 5.67 1 1.33 3.78 0.79 0.5146156
## 204 4.92 1 -0.33 0.72 0.42 1.1820389
## 205 5.18 0 0.24 1.27 0.56 -1.1286687
## 206 5.36 1 0.64 1.90 0.66 0.7257011
## 207 4.86 0 -0.47 0.63 0.39 -0.7915463
## 208 5.59 1 1.15 3.16 0.76 0.5623484
## 209 5.10 1 0.06 1.07 0.52 0.9681799
## 210 4.46 1 -1.35 0.26 0.21 1.9685045
## 211 5.52 1 1.00 2.71 0.73 0.6077328
## 212 4.18 0 -1.98 0.14 0.12 -0.3724214
## 213 4.46 1 -1.35 0.26 0.21 1.9685045
## 214 5.16 1 0.20 1.22 0.55 0.9058664
## 215 4.96 0 -0.25 0.78 0.44 -0.8843608
## 216 5.40 0 0.73 2.07 0.67 -1.4404676
## 217 5.76 1 1.53 4.61 0.82 0.4657416
## 218 4.60 1 -1.04 0.35 0.26 1.6854736
## 219 5.70 0 1.40 4.04 0.80 -2.0089218
## 220 4.86 0 -0.47 0.63 0.39 -0.7915463
## 221 5.94 1 1.93 6.87 0.87 0.3814779
## 222 5.25 0 0.40 1.49 0.60 -1.2197581
## 223 4.98 1 -0.20 0.82 0.45 1.1059612
## 224 5.16 0 0.20 1.22 0.55 -1.1039156
## 225 5.05 1 -0.05 0.95 0.49 1.0233699
## 226 6.46 1 3.08 21.77 0.96 0.2143250
## 227 5.67 1 1.33 3.78 0.79 0.5146156
## 228 5.00 0 -0.16 0.85 0.46 -0.9244656
## 229 5.04 1 -0.07 0.93 0.48 1.0347799
## 230 5.90 0 1.84 6.29 0.86 -2.5076643
## 231 5.70 1 1.40 4.04 0.80 0.4977795
## 232 4.71 0 -0.80 0.45 0.31 -0.6702650
## 233 5.34 1 0.60 1.82 0.64 0.7419735
## 234 5.45 0 0.84 2.32 0.70 -1.5225799
## 235 5.48 1 0.91 2.48 0.71 0.6352928
## 236 5.74 1 1.48 4.41 0.82 0.4761850
## 237 6.14 1 2.37 10.71 0.91 0.3056068
## 238 6.14 1 2.37 10.71 0.91 0.3056068
## 239 4.86 0 -0.47 0.63 0.39 -0.7915463
## 240 4.81 0 -0.58 0.56 0.36 -0.7488584
## 241 5.85 1 1.73 5.63 0.85 0.4215093
## 242 5.55 0 1.06 2.89 0.74 -1.7011134
## 243 4.79 0 -0.62 0.54 0.35 -0.7324350
## 244 5.20 1 0.29 1.33 0.57 0.8665685
## 245 4.95 1 -0.27 0.76 0.43 1.1433675
## 246 5.11 0 0.09 1.09 0.52 -1.0443817
## 247 5.20 1 0.29 1.33 0.57 0.8665685
## 248 4.86 0 -0.47 0.63 0.39 -0.7915463
## 249 4.90 0 -0.38 0.68 0.41 -0.8274420
## 250 5.15 1 0.18 1.19 0.54 0.9159662
## 251 4.77 1 -0.67 0.51 0.34 1.3959233
## 252 4.87 0 -0.45 0.64 0.39 -0.8003715
## 253 4.82 0 -0.56 0.57 0.36 -0.7572077
## 254 5.61 1 1.20 3.31 0.77 0.5500154
## 255 5.26 1 0.42 1.52 0.60 0.8107948
## 256 5.15 0 0.18 1.19 0.54 -1.0917433
## 257 4.57 0 -1.11 0.33 0.25 -0.5738945
## 258 5.48 1 0.91 2.48 0.71 0.6352928
## 259 5.03 1 -0.09 0.91 0.48 1.0463170
## 260 5.07 1 0.00 1.00 0.50 1.0009261
## 261 5.08 1 0.02 1.02 0.51 0.9898894
## 262 5.52 1 1.00 2.71 0.73 0.6077328
## 263 5.22 1 0.33 1.39 0.58 0.8475635
## 264 5.23 1 0.35 1.42 0.59 0.8382179
## 265 5.66 1 1.31 3.69 0.79 0.5203532
## 266 5.77 1 1.55 4.71 0.82 0.4606062
## 267 4.86 1 -0.47 0.63 0.39 1.2633500
## 268 4.92 1 -0.33 0.72 0.42 1.1820389
## 269 5.70 1 1.40 4.04 0.80 0.4977795
## 270 4.90 1 -0.38 0.68 0.41 1.2085439
## 271 5.50 1 0.95 2.59 0.72 0.6213601
## 272 5.96 1 1.97 7.18 0.88 0.3731116
## 273 5.20 1 0.29 1.33 0.57 0.8665685
## 274 5.43 1 0.80 2.22 0.69 0.6715070
## 275 5.50 1 0.95 2.59 0.72 0.6213601
## 276 5.33 1 0.57 1.78 0.64 0.7502461
## 277 4.86 0 -0.47 0.63 0.39 -0.7915463
## 278 4.07 0 -2.22 0.11 0.10 -0.3296599
## 279 5.72 1 1.44 4.22 0.81 0.4868625
## 280 5.17 1 0.22 1.25 0.55 0.8958779
## 281 5.46 1 0.86 2.37 0.70 0.6495381
## 282 5.10 1 0.06 1.07 0.52 0.9681799
## 283 5.35 1 0.62 1.86 0.65 0.7337922
## 284 4.46 0 -1.35 0.26 0.21 -0.5079999
## 285 5.54 1 1.04 2.83 0.74 0.5944045
## 286 4.96 0 -0.25 0.78 0.44 -0.8843608
## 287 5.35 1 0.62 1.86 0.65 0.7337922
## 288 5.16 1 0.20 1.22 0.55 0.9058664
## 289 5.42 1 0.77 2.17 0.68 0.6789939
## 290 5.49 1 0.93 2.53 0.72 0.6282878
## 291 4.41 0 -1.47 0.23 0.19 -0.4806036
## 292 4.99 1 -0.18 0.84 0.46 1.0937663
## 293 4.92 0 -0.33 0.72 0.42 -0.8459958
## 294 5.52 1 1.00 2.71 0.73 0.6077328
## 295 5.04 1 -0.07 0.93 0.48 1.0347799
## 296 5.52 1 1.00 2.71 0.73 0.6077328
## 297 4.90 0 -0.38 0.68 0.41 -0.8274420
## 298 5.00 1 -0.16 0.85 0.46 1.0817060
## 299 4.96 0 -0.25 0.78 0.44 -0.8843608
## 300 5.30 1 0.51 1.66 0.62 0.7756212
## 301 4.61 0 -1.02 0.36 0.26 -0.5999200
## 302 4.56 0 -1.13 0.32 0.24 -0.5675665
## 303 5.23 1 0.35 1.42 0.59 0.8382179
## 304 4.96 0 -0.25 0.78 0.44 -0.8843608
## 305 4.66 0 -0.91 0.40 0.29 -0.6341178
## 306 5.36 0 0.64 1.90 0.66 -1.3779778
## 307 4.46 0 -1.35 0.26 0.21 -0.5079999
## 308 4.07 0 -2.22 0.11 0.10 -0.3296599
## 309 4.46 0 -1.35 0.26 0.21 -0.5079999
## 310 5.76 1 1.53 4.61 0.82 0.4657416
## 311 4.77 0 -0.67 0.51 0.34 -0.7163717
## 312 5.18 0 0.24 1.27 0.56 -1.1286687
## 313 5.18 1 0.24 1.27 0.56 0.8859996
## 314 4.90 1 -0.38 0.68 0.41 1.2085439
## 315 4.96 0 -0.25 0.78 0.44 -0.8843608
## 316 5.07 1 0.00 1.00 0.50 1.0009261
## 317 4.75 0 -0.71 0.49 0.33 -0.7006608
## 318 5.57 1 1.11 3.03 0.75 0.5749580
## 319 5.48 1 0.91 2.48 0.71 0.6352928
## 320 4.39 1 -1.51 0.22 0.18 2.1273730
## 321 4.85 1 -0.49 0.61 0.38 1.2774356
## 322 4.32 0 -1.67 0.19 0.16 -0.4349598
## 323 5.31 1 0.53 1.70 0.63 0.7670689
## 324 4.87 0 -0.45 0.64 0.39 -0.8003715
## 325 5.50 1 0.95 2.59 0.72 0.6213601
## 326 5.45 0 0.84 2.32 0.70 -1.5225799
## 327 5.59 1 1.15 3.16 0.76 0.5623484
## 328 5.83 1 1.68 5.38 0.84 0.4309609
## 329 4.87 0 -0.45 0.64 0.39 -0.8003715
## 330 5.12 1 0.11 1.12 0.53 0.9469464
## 331 5.50 1 0.95 2.59 0.72 0.6213601
## 332 5.39 1 0.71 2.03 0.67 0.7019591
## 333 5.83 1 1.68 5.38 0.84 0.4309609
## 334 5.45 0 0.84 2.32 0.70 -1.5225799
## 335 5.06 0 -0.02 0.98 0.49 -0.9880585
## 336 5.00 0 -0.16 0.85 0.46 -0.9244656
## 337 4.96 0 -0.25 0.78 0.44 -0.8843608
## 338 5.18 0 0.24 1.27 0.56 -1.1286687
## 339 5.16 0 0.20 1.22 0.55 -1.1039156
## 340 4.31 0 -1.69 0.19 0.16 -0.4301638
## 341 4.74 1 -0.73 0.48 0.32 1.4431368
## 342 4.94 1 -0.29 0.75 0.43 1.1561153
## 343 5.18 0 0.24 1.27 0.56 -1.1286687
## 344 5.36 0 0.64 1.90 0.66 -1.3779778
## 345 5.40 1 0.73 2.07 0.67 0.6942190
## 346 4.89 1 -0.40 0.67 0.40 1.2220184
## 347 4.83 0 -0.53 0.59 0.37 -0.7656501
## 348 5.91 1 1.86 6.43 0.87 0.3943804
## 349 5.15 1 0.18 1.19 0.54 0.9159662
## 350 5.54 1 1.04 2.83 0.74 0.5944045
## 351 5.00 1 -0.16 0.85 0.46 1.0817060
## 352 5.19 0 0.26 1.30 0.57 -1.1412527
## 353 5.08 1 0.02 1.02 0.51 0.9898894
## 354 5.64 1 1.26 3.53 0.78 0.5320211
## 355 5.05 1 -0.05 0.95 0.49 1.0233699
## 356 5.39 1 0.71 2.03 0.67 0.7019591
## 357 5.37 1 0.66 1.94 0.66 0.7176992
## 358 4.80 1 -0.60 0.55 0.35 1.3502544
## 359 4.51 0 -1.24 0.29 0.22 -0.5369578
## 360 5.00 0 -0.16 0.85 0.46 -0.9244656
## 361 4.99 0 -0.18 0.84 0.46 -0.9142721
## 362 5.64 1 1.26 3.53 0.78 0.5320211
## 363 4.79 0 -0.62 0.54 0.35 -0.7324350
## 364 4.80 1 -0.60 0.55 0.35 1.3502544
## 365 4.73 1 -0.76 0.47 0.32 1.4592269
## 366 5.05 0 -0.05 0.95 0.49 -0.9771638
## 367 5.00 1 -0.16 0.85 0.46 1.0817060
## 368 5.44 0 0.82 2.27 0.69 -1.5057912
## 369 5.29 1 0.49 1.63 0.62 0.7842689
## 370 5.55 0 1.06 2.89 0.74 -1.7011134
## 371 5.24 1 0.38 1.46 0.59 0.8289753
## 372 4.94 0 -0.29 0.75 0.43 -0.8649656
## 373 5.20 1 0.29 1.33 0.57 0.8665685
## 374 5.91 1 1.86 6.43 0.87 0.3943804
## 375 5.26 1 0.42 1.52 0.60 0.8107948
## 376 4.98 1 -0.20 0.82 0.45 1.1059612
## 377 5.07 0 0.00 1.00 0.50 -0.9990748
## 378 6.30 1 2.73 15.27 0.94 0.2559280
## 379 5.58 1 1.13 3.09 0.76 0.5686182
## 380 5.33 1 0.57 1.78 0.64 0.7502461
## 381 5.23 1 0.35 1.42 0.59 0.8382179
## 382 5.42 1 0.77 2.17 0.68 0.6789939
## 383 4.90 1 -0.38 0.68 0.41 1.2085439
## 384 5.48 1 0.91 2.48 0.71 0.6352928
## 385 5.16 1 0.20 1.22 0.55 0.9058664
## 386 4.94 0 -0.29 0.75 0.43 -0.8649656
## 387 5.22 1 0.33 1.39 0.58 0.8475635
## 388 5.40 0 0.73 2.07 0.67 -1.4404676
## 389 4.87 1 -0.45 0.64 0.39 1.2494198
## 390 5.22 0 0.33 1.39 0.58 -1.1798526
## 391 5.15 0 0.18 1.19 0.54 -1.0917433
## 392 5.10 1 0.06 1.07 0.52 0.9681799
## 393 5.06 0 -0.02 0.98 0.49 -0.9880585
## 394 5.01 1 -0.13 0.87 0.47 1.0697786
## 395 4.93 1 -0.31 0.73 0.42 1.1690053
## 396 4.61 0 -1.02 0.36 0.26 -0.5999200
## 397 5.96 1 1.97 7.18 0.88 0.3731116
## 398 4.68 0 -0.87 0.42 0.30 -0.6483367
## 399 5.14 1 0.15 1.17 0.54 0.9261787
## 400 5.22 0 0.33 1.39 0.58 -1.1798526
## 401 5.37 1 0.66 1.94 0.66 0.7176992
## 402 5.00 1 -0.16 0.85 0.46 1.0817060
## 403 4.35 0 -1.60 0.20 0.17 -0.4496712
## 404 6.11 1 2.30 10.02 0.91 0.3159431
## 405 5.59 0 1.15 3.16 0.76 -1.7782570
## 406 4.87 0 -0.45 0.64 0.39 -0.8003715
## 407 4.98 1 -0.20 0.82 0.45 1.1059612
## 408 5.72 1 1.44 4.22 0.81 0.4868625
## 409 5.07 1 0.00 1.00 0.50 1.0009261
## 410 3.93 0 -2.53 0.08 0.07 -0.2822616
## 411 6.07 1 2.22 9.17 0.90 0.3302708
## 412 5.96 1 1.97 7.18 0.88 0.3731116
## 413 5.36 1 0.64 1.90 0.66 0.7257011
## 414 5.55 0 1.06 2.89 0.74 -1.7011134
## 415 4.90 0 -0.38 0.68 0.41 -0.8274420
## 416 5.64 0 1.26 3.53 0.78 -1.8796246
## 417 5.19 0 0.26 1.30 0.57 -1.1412527
## 418 4.17 0 -2.00 0.14 0.12 -0.3683150
## 419 5.55 0 1.06 2.89 0.74 -1.7011134
## 420 5.88 1 1.79 6.02 0.86 0.4077193
## 421 5.27 0 0.44 1.56 0.61 -1.2471088
## 422 4.80 1 -0.60 0.55 0.35 1.3502544
## 423 5.81 0 1.64 5.15 0.84 -2.2695070
## 424 4.51 0 -1.24 0.29 0.22 -0.5369578
## 425 5.57 1 1.11 3.03 0.75 0.5749580
## 426 5.36 0 0.64 1.90 0.66 -1.3779778
## 427 5.52 0 1.00 2.71 0.73 -1.6454599
## 428 5.26 1 0.42 1.52 0.60 0.8107948
## 429 5.33 1 0.57 1.78 0.64 0.7502461
## 430 5.13 1 0.13 1.14 0.53 0.9365050
## 431 5.20 1 0.29 1.33 0.57 0.8665685
## 432 5.49 1 0.93 2.53 0.72 0.6282878
## 433 5.30 0 0.51 1.66 0.62 -1.2892891
## 434 4.95 1 -0.27 0.76 0.43 1.1433675
## 435 5.42 1 0.77 2.17 0.68 0.6789939
## 436 5.08 1 0.02 1.02 0.51 0.9898894
## 437 5.11 0 0.09 1.09 0.52 -1.0443817
## 438 5.73 1 1.46 4.31 0.81 0.4814941
## 439 6.32 1 2.77 15.96 0.94 0.2503152
## 440 5.23 1 0.35 1.42 0.59 0.8382179
## 441 5.21 0 0.31 1.36 0.58 -1.1668431
## 442 4.71 0 -0.80 0.45 0.31 -0.6702650
## 443 4.91 0 -0.36 0.70 0.41 -0.8366675
## 444 5.33 0 0.57 1.78 0.64 -1.3328960
## 445 5.33 1 0.57 1.78 0.64 0.7502461
## 446 5.18 1 0.24 1.27 0.56 0.8859996
## 447 5.11 0 0.09 1.09 0.52 -1.0443817
## 448 5.00 0 -0.16 0.85 0.46 -0.9244656
## 449 6.06 1 2.19 8.97 0.90 0.3339531
## 450 4.97 1 -0.22 0.80 0.44 1.1182919
## 451 5.40 1 0.73 2.07 0.67 0.6942190
## 452 5.74 0 1.48 4.41 0.82 -2.1000242
## 453 4.82 0 -0.56 0.57 0.36 -0.7572077
## 454 5.07 1 0.00 1.00 0.50 1.0009261
## 455 4.46 0 -1.35 0.26 0.21 -0.5079999
## 456 4.83 0 -0.53 0.59 0.37 -0.7656501
## 457 4.62 1 -1.00 0.37 0.27 1.6485090
## 458 5.49 1 0.93 2.53 0.72 0.6282878
## 459 4.89 0 -0.40 0.67 0.40 -0.8183183
## 460 4.26 0 -1.80 0.17 0.14 -0.4069652
## 461 5.47 0 0.89 2.42 0.71 -1.5567208
## 462 4.70 0 -0.82 0.44 0.31 -0.6628744
## 463 5.58 0 1.13 3.09 0.76 -1.7586492
## 464 5.54 1 1.04 2.83 0.74 0.5944045
## 465 5.71 1 1.42 4.13 0.80 0.4922907
## 466 5.15 1 0.18 1.19 0.54 0.9159662
## 467 5.18 0 0.24 1.27 0.56 -1.1286687
## 468 5.25 0 0.40 1.49 0.60 -1.2197581
## 469 5.83 0 1.68 5.38 0.84 -2.3203963
## 470 6.37 1 2.88 17.83 0.95 0.2368158
## 471 5.11 0 0.09 1.09 0.52 -1.0443817
## 472 4.68 0 -0.87 0.42 0.30 -0.6483367
## 473 5.48 0 0.91 2.48 0.71 -1.5740772
## 474 5.14 1 0.15 1.17 0.54 0.9261787
## 475 5.14 1 0.15 1.17 0.54 0.9261787
## 476 5.01 1 -0.13 0.87 0.47 1.0697786
## 477 5.00 1 -0.16 0.85 0.46 1.0817060
## 478 4.92 1 -0.33 0.72 0.42 1.1820389
## 479 5.14 1 0.15 1.17 0.54 0.9261787
## 480 5.15 1 0.18 1.19 0.54 0.9159662
## 481 5.03 1 -0.09 0.91 0.48 1.0463170
## 482 5.01 0 -0.13 0.87 0.47 -0.9347728
## 483 5.40 0 0.73 2.07 0.67 -1.4404676
## 484 5.64 1 1.26 3.53 0.78 0.5320211
## 485 3.93 0 -2.53 0.08 0.07 -0.2822616
## 486 5.20 0 0.29 1.33 0.57 -1.1539769
## 487 5.73 1 1.46 4.31 0.81 0.4814941
## 488 4.63 0 -0.98 0.38 0.27 -0.6133721
## 489 5.32 1 0.55 1.74 0.63 0.7586109
## 490 5.72 0 1.44 4.22 0.81 -2.0539680
## 491 4.59 0 -1.07 0.34 0.26 -0.5867630
## 492 5.04 0 -0.07 0.93 0.48 -0.9663891
## 493 4.69 0 -0.84 0.43 0.30 -0.6555652
## 494 5.50 0 0.95 2.59 0.72 -1.6093728
## 495 5.72 1 1.44 4.22 0.81 0.4868625
## 496 5.55 0 1.06 2.89 0.74 -1.7011134
## 497 5.25 1 0.40 1.49 0.60 0.8198347
## 498 5.21 0 0.31 1.36 0.58 -1.1668431
## 499 5.14 1 0.15 1.17 0.54 0.9261787
## 500 5.40 0 0.73 2.07 0.67 -1.4404676
## 501 4.61 0 -1.02 0.36 0.26 -0.5999200
## 502 5.11 0 0.09 1.09 0.52 -1.0443817
## 503 5.42 1 0.77 2.17 0.68 0.6789939
## 504 5.37 1 0.66 1.94 0.66 0.7176992
## 505 5.12 0 0.11 1.12 0.53 -1.0560259
## 506 5.37 1 0.66 1.94 0.66 0.7176992
## 507 5.18 1 0.24 1.27 0.56 0.8859996
## 508 5.59 1 1.15 3.16 0.76 0.5623484
## 509 5.00 1 -0.16 0.85 0.46 1.0817060
## 510 6.32 0 2.77 15.96 0.94 -3.9949630
## 511 5.46 1 0.86 2.37 0.70 0.6495381
## 512 4.98 1 -0.20 0.82 0.45 1.1059612
## 513 4.96 0 -0.25 0.78 0.44 -0.8843608
## 514 5.27 1 0.44 1.56 0.61 0.8018546
## 515 5.72 0 1.44 4.22 0.81 -2.0539680
## 516 5.49 0 0.93 2.53 0.72 -1.5916272
## 517 5.07 0 0.00 1.00 0.50 -0.9990748
## 518 4.41 0 -1.47 0.23 0.19 -0.4806036
## 519 5.54 1 1.04 2.83 0.74 0.5944045
## 520 5.30 1 0.51 1.66 0.62 0.7756212
## 521 5.40 1 0.73 2.07 0.67 0.6942190
## 522 5.09 1 0.04 1.04 0.51 0.9789745
## 523 6.06 1 2.19 8.97 0.90 0.3339531
## 524 5.44 0 0.82 2.27 0.69 -1.5057912
## 525 4.61 0 -1.02 0.36 0.26 -0.5999200
## 526 5.72 1 1.44 4.22 0.81 0.4868625
## 527 4.94 1 -0.29 0.75 0.43 1.1561153
## 528 4.51 0 -1.24 0.29 0.22 -0.5369578
## 529 5.37 1 0.66 1.94 0.66 0.7176992
## 530 5.11 0 0.09 1.09 0.52 -1.0443817
## 531 5.19 0 0.26 1.30 0.57 -1.1412527
## 532 5.25 0 0.40 1.49 0.60 -1.2197581
## 533 3.96 0 -2.46 0.09 0.08 -0.2918083
## 534 5.29 1 0.49 1.63 0.62 0.7842689
## 535 4.82 0 -0.56 0.57 0.36 -0.7572077
## 536 4.61 1 -1.02 0.36 0.26 1.6668889
## 537 5.87 1 1.77 5.88 0.85 0.4122651
## 538 5.67 1 1.33 3.78 0.79 0.5146156
## 539 5.63 1 1.24 3.46 0.78 0.5379528
## 540 4.41 1 -1.47 0.23 0.19 2.0807170
## 541 5.01 1 -0.13 0.87 0.47 1.0697786
## 542 4.79 0 -0.62 0.54 0.35 -0.7324350
## 543 4.70 1 -0.82 0.44 0.31 1.5085815
## 544 5.36 1 0.64 1.90 0.66 0.7257011
## 545 5.43 1 0.80 2.22 0.69 0.6715070
## 546 5.83 1 1.68 5.38 0.84 0.4309609
## 547 5.20 1 0.29 1.33 0.57 0.8665685
## 548 5.11 0 0.09 1.09 0.52 -1.0443817
## 549 5.52 1 1.00 2.71 0.73 0.6077328
## 550 6.26 1 2.64 13.97 0.93 0.2675341
## 551 4.67 0 -0.89 0.41 0.29 -0.6411878
## 552 5.48 1 0.91 2.48 0.71 0.6352928
## 553 4.68 1 -0.87 0.42 0.30 1.5424085
## 554 4.32 0 -1.67 0.19 0.16 -0.4349598
## 555 6.06 1 2.19 8.97 0.90 0.3339531
## 556 5.22 1 0.33 1.39 0.58 0.8475635
## 557 4.71 0 -0.80 0.45 0.31 -0.6702650
## 558 4.94 0 -0.29 0.75 0.43 -0.8649656
## 559 5.88 1 1.79 6.02 0.86 0.4077193
## 560 5.50 0 0.95 2.59 0.72 -1.6093728
## 561 5.33 0 0.57 1.78 0.64 -1.3328960
## 562 4.90 1 -0.38 0.68 0.41 1.2085439
## 563 4.21 0 -1.91 0.15 0.13 -0.3850176
## 564 5.35 1 0.62 1.86 0.65 0.7337922
## 565 5.25 1 0.40 1.49 0.60 0.8198347
## 566 4.80 1 -0.60 0.55 0.35 1.3502544
## 567 4.89 1 -0.40 0.67 0.40 1.2220184
## 568 5.00 0 -0.16 0.85 0.46 -0.9244656
## 569 6.26 1 2.64 13.97 0.93 0.2675341
## 570 4.30 0 -1.71 0.18 0.15 -0.4254206
## 571 5.25 0 0.40 1.49 0.60 -1.2197581
## 572 5.67 1 1.33 3.78 0.79 0.5146156
## 573 5.29 1 0.49 1.63 0.62 0.7842689
## 574 5.10 1 0.06 1.07 0.52 0.9681799
## 575 4.46 0 -1.35 0.26 0.21 -0.5079999
## 576 5.01 0 -0.13 0.87 0.47 -0.9347728
## 577 5.98 0 2.02 7.51 0.88 -2.7402611
## 578 4.61 0 -1.02 0.36 0.26 -0.5999200
## 579 5.39 1 0.71 2.03 0.67 0.7019591
## 580 5.26 1 0.42 1.52 0.60 0.8107948
## 581 5.15 0 0.18 1.19 0.54 -1.0917433
## 582 4.81 0 -0.58 0.56 0.36 -0.7488584
## 583 4.85 0 -0.49 0.61 0.38 -0.7828183
## 584 4.82 0 -0.56 0.57 0.36 -0.7572077
## 585 5.46 1 0.86 2.37 0.70 0.6495381
## 586 5.50 0 0.95 2.59 0.72 -1.6093728
## 587 5.11 0 0.09 1.09 0.52 -1.0443817
## 588 4.90 0 -0.38 0.68 0.41 -0.8274420
## 589 6.37 1 2.88 17.83 0.95 0.2368158
## 590 4.86 0 -0.47 0.63 0.39 -0.7915463
## 591 5.05 0 -0.05 0.95 0.49 -0.9771638
## 592 5.28 0 0.46 1.59 0.61 -1.2610133
## 593 5.52 1 1.00 2.71 0.73 0.6077328
## 594 4.83 0 -0.53 0.59 0.37 -0.7656501
## 595 5.30 1 0.51 1.66 0.62 0.7756212
## 596 5.16 0 0.20 1.22 0.55 -1.1039156
## 597 4.93 1 -0.31 0.73 0.42 1.1690053
## 598 5.55 1 1.06 2.89 0.74 0.5878503
## 599 4.80 0 -0.60 0.55 0.35 -0.7406012
## 600 4.96 0 -0.25 0.78 0.44 -0.8843608
## 601 5.58 1 1.13 3.09 0.76 0.5686182
## 602 5.22 0 0.33 1.39 0.58 -1.1798526
## 603 4.18 0 -1.98 0.14 0.12 -0.3724214
## 604 5.62 1 1.22 3.38 0.77 0.5439507
## 605 5.11 0 0.09 1.09 0.52 -1.0443817
## 606 5.07 1 0.00 1.00 0.50 1.0009261
## 607 4.44 0 -1.40 0.25 0.20 -0.4968588
## 608 5.47 1 0.89 2.42 0.71 0.6423760
## 609 5.41 1 0.75 2.12 0.68 0.6865643
## 610 5.34 1 0.60 1.82 0.64 0.7419735
## 611 4.80 0 -0.60 0.55 0.35 -0.7406012
## 612 4.07 0 -2.22 0.11 0.10 -0.3296599
## 613 5.22 0 0.33 1.39 0.58 -1.1798526
## 614 5.64 0 1.26 3.53 0.78 -1.8796246
## 615 5.80 1 1.62 5.04 0.83 0.4455370
## 616 4.67 0 -0.89 0.41 0.29 -0.6411878
## 617 5.23 1 0.35 1.42 0.59 0.8382179
## 618 5.03 1 -0.09 0.91 0.48 1.0463170
## 619 5.17 1 0.22 1.25 0.55 0.8958779
## 620 4.85 1 -0.49 0.61 0.38 1.2774356
## 621 5.98 1 2.02 7.51 0.88 0.3649287
## 622 5.39 1 0.71 2.03 0.67 0.7019591
## 623 5.70 1 1.40 4.04 0.80 0.4977795
## 624 4.93 1 -0.31 0.73 0.42 1.1690053
## 625 5.98 1 2.02 7.51 0.88 0.3649287
## 626 5.15 0 0.18 1.19 0.54 -1.0917433
## 627 5.47 1 0.89 2.42 0.71 0.6423760
## 628 5.50 1 0.95 2.59 0.72 0.6213601
## 629 5.36 0 0.64 1.90 0.66 -1.3779778
## 630 5.66 1 1.31 3.69 0.79 0.5203532
## 631 4.76 0 -0.69 0.50 0.33 -0.7084727
## 632 6.06 1 2.19 8.97 0.90 0.3339531
## 633 4.21 0 -1.91 0.15 0.13 -0.3850176
## 634 4.26 0 -1.80 0.17 0.14 -0.4069652
## 635 5.93 1 1.91 6.72 0.87 0.3857311
## 636 5.96 1 1.97 7.18 0.88 0.3731116
## 637 5.55 1 1.06 2.89 0.74 0.5878503
## 638 5.71 1 1.42 4.13 0.80 0.4922907
## 639 5.15 1 0.18 1.19 0.54 0.9159662
## 640 5.46 1 0.86 2.37 0.70 0.6495381
## 641 5.15 1 0.18 1.19 0.54 0.9159662
## 642 5.77 1 1.55 4.71 0.82 0.4606062
## 643 4.42 1 -1.44 0.24 0.19 2.0577741
## 644 5.70 1 1.40 4.04 0.80 0.4977795
## 645 5.50 1 0.95 2.59 0.72 0.6213601
## 646 5.55 1 1.06 2.89 0.74 0.5878503
## 647 5.24 1 0.38 1.46 0.59 0.8289753
## 648 5.05 0 -0.05 0.95 0.49 -0.9771638
## 649 4.72 0 -0.78 0.46 0.31 -0.6777380
## 650 4.94 0 -0.29 0.75 0.43 -0.8649656
## 651 4.72 1 -0.78 0.46 0.31 1.4754964
## 652 5.28 0 0.46 1.59 0.61 -1.2610133
## 653 5.08 1 0.02 1.02 0.51 0.9898894
## 654 5.04 0 -0.07 0.93 0.48 -0.9663891
## 655 4.85 0 -0.49 0.61 0.38 -0.7828183
## 656 4.89 0 -0.40 0.67 0.40 -0.8183183
## 657 5.81 0 1.64 5.15 0.84 -2.2695070
## 658 4.52 1 -1.22 0.29 0.23 1.8418086
## 659 4.48 0 -1.31 0.27 0.21 -0.5193908
## 660 4.90 1 -0.38 0.68 0.41 1.2085439
## 661 4.74 1 -0.73 0.48 0.32 1.4431368
## 662 4.71 0 -0.80 0.45 0.31 -0.6702650
## 663 4.48 1 -1.31 0.27 0.21 1.9253327
## 664 5.35 1 0.62 1.86 0.65 0.7337922
## 665 4.74 0 -0.73 0.48 0.32 -0.6929350
## 666 5.40 1 0.73 2.07 0.67 0.6942190
## 667 4.61 0 -1.02 0.36 0.26 -0.5999200
## 668 4.77 0 -0.67 0.51 0.34 -0.7163717
## 669 5.50 0 0.95 2.59 0.72 -1.6093728
## 670 3.96 0 -2.46 0.09 0.08 -0.2918083
## 671 5.02 1 -0.11 0.89 0.47 1.0579828
## 672 4.86 0 -0.47 0.63 0.39 -0.7915463
## 673 4.85 0 -0.49 0.61 0.38 -0.7828183
## 674 5.49 1 0.93 2.53 0.72 0.6282878
## 675 5.23 1 0.35 1.42 0.59 0.8382179
## 676 4.66 0 -0.91 0.40 0.29 -0.6341178
## 677 4.85 1 -0.49 0.61 0.38 1.2774356
## 678 5.13 0 0.13 1.14 0.53 -1.0678000
## 679 4.61 0 -1.02 0.36 0.26 -0.5999200
## 680 5.54 1 1.04 2.83 0.74 0.5944045
## 681 5.40 0 0.73 2.07 0.67 -1.4404676
## 682 5.65 1 1.28 3.61 0.78 0.5261548
## 683 4.86 0 -0.47 0.63 0.39 -0.7915463
## 684 5.14 1 0.15 1.17 0.54 0.9261787
## 685 5.25 0 0.40 1.49 0.60 -1.2197581
## 686 5.31 0 0.53 1.70 0.63 -1.3036639
## 687 5.36 1 0.64 1.90 0.66 0.7257011
## 688 4.85 1 -0.49 0.61 0.38 1.2774356
## 689 5.43 1 0.80 2.22 0.69 0.6715070
## 690 5.13 0 0.13 1.14 0.53 -1.0678000
## 691 5.67 1 1.33 3.78 0.79 0.5146156
## 692 5.30 0 0.51 1.66 0.62 -1.2892891
## 693 4.51 0 -1.24 0.29 0.22 -0.5369578
## 694 5.48 1 0.91 2.48 0.71 0.6352928
## 695 5.00 0 -0.16 0.85 0.46 -0.9244656
## 696 5.48 1 0.91 2.48 0.71 0.6352928
## 697 3.93 0 -2.53 0.08 0.07 -0.2822616
## 698 5.58 1 1.13 3.09 0.76 0.5686182
## 699 5.31 0 0.53 1.70 0.63 -1.3036639
## 700 4.42 0 -1.44 0.24 0.19 -0.4859620
## 701 5.20 0 0.29 1.33 0.57 -1.1539769
## 702 5.26 1 0.42 1.52 0.60 0.8107948
## 703 5.40 0 0.73 2.07 0.67 -1.4404676
## 704 5.35 1 0.62 1.86 0.65 0.7337922
## 705 5.14 1 0.15 1.17 0.54 0.9261787
## 706 5.16 0 0.20 1.22 0.55 -1.1039156
## 707 4.97 1 -0.22 0.80 0.44 1.1182919
## 708 5.11 1 0.09 1.09 0.52 0.9575043
## 709 5.50 1 0.95 2.59 0.72 0.6213601
## 710 4.18 0 -1.98 0.14 0.12 -0.3724214
## 711 5.64 1 1.26 3.53 0.78 0.5320211
## 712 5.32 0 0.55 1.74 0.63 -1.3181989
## 713 5.00 0 -0.16 0.85 0.46 -0.9244656
## 714 5.13 0 0.13 1.14 0.53 -1.0678000
## 715 5.72 1 1.44 4.22 0.81 0.4868625
## 716 5.14 1 0.15 1.17 0.54 0.9261787
## 717 5.42 1 0.77 2.17 0.68 0.6789939
## 718 5.09 0 0.04 1.04 0.51 -1.0214771
## 719 4.85 0 -0.49 0.61 0.38 -0.7828183
## 720 5.64 1 1.26 3.53 0.78 0.5320211
## 721 5.50 1 0.95 2.59 0.72 0.6213601
## 722 5.01 0 -0.13 0.87 0.47 -0.9347728
## 723 4.59 1 -1.07 0.34 0.26 1.7042656
## 724 4.96 1 -0.25 0.78 0.44 1.1307602
## 725 5.46 1 0.86 2.37 0.70 0.6495381
## 726 5.70 1 1.40 4.04 0.80 0.4977795
## 727 5.98 0 2.02 7.51 0.88 -2.7402611
## 728 5.85 0 1.73 5.63 0.85 -2.3724267
## 729 5.97 0 1.99 7.34 0.88 -2.7100457
## 730 5.74 1 1.48 4.41 0.82 0.4761850
## 731 5.25 0 0.40 1.49 0.60 -1.2197581
## 732 5.22 1 0.33 1.39 0.58 0.8475635
## 733 5.23 1 0.35 1.42 0.59 0.8382179
## 734 4.46 0 -1.35 0.26 0.21 -0.5079999
## 735 5.42 1 0.77 2.17 0.68 0.6789939
## 736 5.70 1 1.40 4.04 0.80 0.4977795
## 737 4.16 0 -2.02 0.13 0.12 -0.3642538
## 738 4.60 0 -1.04 0.35 0.26 -0.5933050
## 739 5.85 1 1.73 5.63 0.85 0.4215093
## 740 4.57 0 -1.11 0.33 0.25 -0.5738945
## 741 4.68 1 -0.87 0.42 0.30 1.5424085
## 742 5.50 1 0.95 2.59 0.72 0.6213601
## 743 4.26 0 -1.80 0.17 0.14 -0.4069652
## 744 4.83 0 -0.53 0.59 0.37 -0.7656501
## 745 5.00 1 -0.16 0.85 0.46 1.0817060
## 746 5.11 0 0.09 1.09 0.52 -1.0443817
## 747 5.22 1 0.33 1.39 0.58 0.8475635
## 748 5.20 1 0.29 1.33 0.57 0.8665685
## 749 5.09 1 0.04 1.04 0.51 0.9789745
## 750 4.96 0 -0.25 0.78 0.44 -0.8843608
## 751 5.40 0 0.73 2.07 0.67 -1.4404676
## 752 5.40 1 0.73 2.07 0.67 0.6942190
## 753 5.15 0 0.18 1.19 0.54 -1.0917433
## 754 5.30 0 0.51 1.66 0.62 -1.2892891
## 755 5.21 1 0.31 1.36 0.58 0.8570133
## 756 5.28 1 0.46 1.59 0.61 0.7930130
## 757 4.26 0 -1.80 0.17 0.14 -0.4069652
## 758 4.86 0 -0.47 0.63 0.39 -0.7915463
## 759 5.18 0 0.24 1.27 0.56 -1.1286687
## 760 5.70 0 1.40 4.04 0.80 -2.0089218
## 761 5.69 0 1.37 3.95 0.80 -1.9867706
## 762 5.56 0 1.08 2.96 0.75 -1.7200797
## 763 5.65 1 1.28 3.61 0.78 0.5261548
## 764 4.79 0 -0.62 0.54 0.35 -0.7324350
## 765 4.80 1 -0.60 0.55 0.35 1.3502544
## 766 4.78 1 -0.64 0.52 0.34 1.3805312
## 767 5.18 1 0.24 1.27 0.56 0.8859996
## 768 6.37 1 2.88 17.83 0.95 0.2368158
## 769 5.89 1 1.82 6.15 0.86 0.4032236
## 770 5.11 0 0.09 1.09 0.52 -1.0443817
## 771 4.82 0 -0.56 0.57 0.36 -0.7572077
## 772 4.82 0 -0.56 0.57 0.36 -0.7572077
## 773 4.80 1 -0.60 0.55 0.35 1.3502544
## 774 5.15 1 0.18 1.19 0.54 0.9159662
## 775 5.36 1 0.64 1.90 0.66 0.7257011
## 776 4.71 0 -0.80 0.45 0.31 -0.6702650
## 777 4.90 0 -0.38 0.68 0.41 -0.8274420
## 778 5.37 1 0.66 1.94 0.66 0.7176992
## 779 4.55 0 -1.15 0.32 0.24 -0.5613083
## 780 4.98 0 -0.20 0.82 0.45 -0.9041909
## 781 4.42 0 -1.44 0.24 0.19 -0.4859620
## 782 5.42 1 0.77 2.17 0.68 0.6789939
## 783 4.96 0 -0.25 0.78 0.44 -0.8843608
## 784 5.87 1 1.77 5.88 0.85 0.4122651
## 785 4.71 0 -0.80 0.45 0.31 -0.6702650
## 786 5.40 0 0.73 2.07 0.67 -1.4404676
## 787 5.55 1 1.06 2.89 0.74 0.5878503
## 788 5.15 1 0.18 1.19 0.54 0.9159662
## 789 6.14 1 2.37 10.71 0.91 0.3056068
## 790 4.89 0 -0.40 0.67 0.40 -0.8183183
## 791 5.40 0 0.73 2.07 0.67 -1.4404676
## 792 4.07 0 -2.22 0.11 0.10 -0.3296599
## 793 4.90 0 -0.38 0.68 0.41 -0.8274420
## 794 4.66 0 -0.91 0.40 0.29 -0.6341178
## 795 5.23 1 0.35 1.42 0.59 0.8382179
## 796 5.18 1 0.24 1.27 0.56 0.8859996
## 797 5.72 1 1.44 4.22 0.81 0.4868625
## 798 5.50 0 0.95 2.59 0.72 -1.6093728
## 799 4.59 0 -1.07 0.34 0.26 -0.5867630
## 800 5.29 1 0.49 1.63 0.62 0.7842689
## 801 4.92 1 -0.33 0.72 0.42 1.1820389
## 802 5.30 0 0.51 1.66 0.62 -1.2892891
## 803 5.26 1 0.42 1.52 0.60 0.8107948
## 804 5.15 0 0.18 1.19 0.54 -1.0917433
## 805 5.35 1 0.62 1.86 0.65 0.7337922
## 806 4.93 0 -0.31 0.73 0.42 -0.8554281
## 807 5.61 1 1.20 3.31 0.77 0.5500154
## 808 5.70 1 1.40 4.04 0.80 0.4977795
## 809 4.96 0 -0.25 0.78 0.44 -0.8843608
## 810 5.26 0 0.42 1.52 0.60 -1.2333577
## 811 4.86 0 -0.47 0.63 0.39 -0.7915463
## 812 5.64 1 1.26 3.53 0.78 0.5320211
## 813 5.18 0 0.24 1.27 0.56 -1.1286687
## 814 5.20 1 0.29 1.33 0.57 0.8665685
## 815 5.24 1 0.38 1.46 0.59 0.8289753
## 816 5.48 1 0.91 2.48 0.71 0.6352928
## 817 4.87 0 -0.45 0.64 0.39 -0.8003715
## 818 5.04 0 -0.07 0.93 0.48 -0.9663891
## 819 5.11 0 0.09 1.09 0.52 -1.0443817
## 820 6.06 1 2.19 8.97 0.90 0.3339531
## 821 5.29 1 0.49 1.63 0.62 0.7842689
## 822 5.03 1 -0.09 0.91 0.48 1.0463170
## 823 5.25 1 0.40 1.49 0.60 0.8198347
## 824 4.70 1 -0.82 0.44 0.31 1.5085815
## 825 5.15 1 0.18 1.19 0.54 0.9159662
## 826 4.29 0 -1.73 0.18 0.15 -0.4207297
## 827 4.93 1 -0.31 0.73 0.42 1.1690053
## 828 5.05 0 -0.05 0.95 0.49 -0.9771638
## 829 4.68 0 -0.87 0.42 0.30 -0.6483367
## 830 5.05 0 -0.05 0.95 0.49 -0.9771638
## 831 5.75 1 1.51 4.51 0.82 0.4709344
## 832 5.54 1 1.04 2.83 0.74 0.5944045
## 833 5.73 1 1.46 4.31 0.81 0.4814941
## 834 5.30 0 0.51 1.66 0.62 -1.2892891
## 835 4.66 0 -0.91 0.40 0.29 -0.6341178
## 836 5.58 1 1.13 3.09 0.76 0.5686182
## 837 5.22 0 0.33 1.39 0.58 -1.1798526
## 838 5.08 1 0.02 1.02 0.51 0.9898894
## 839 5.67 1 1.33 3.78 0.79 0.5146156
## 840 5.28 1 0.46 1.59 0.61 0.7930130
## 841 4.82 0 -0.56 0.57 0.36 -0.7572077
## 842 5.12 1 0.11 1.12 0.53 0.9469464
## 843 5.64 1 1.26 3.53 0.78 0.5320211
## 844 4.86 0 -0.47 0.63 0.39 -0.7915463
## 845 4.86 0 -0.47 0.63 0.39 -0.7915463
## 846 5.22 1 0.33 1.39 0.58 0.8475635
## 847 5.55 0 1.06 2.89 0.74 -1.7011134
## 848 5.07 1 0.00 1.00 0.50 1.0009261
## 849 5.18 1 0.24 1.27 0.56 0.8859996
## 850 5.70 1 1.40 4.04 0.80 0.4977795
## 851 4.90 0 -0.38 0.68 0.41 -0.8274420
## 852 5.72 1 1.44 4.22 0.81 0.4868625
## 853 5.79 1 1.59 4.93 0.83 0.4505045
## 854 5.13 1 0.13 1.14 0.53 0.9365050
## 855 5.30 0 0.51 1.66 0.62 -1.2892891
## 856 5.08 0 0.02 1.02 0.51 -1.0102138
## 857 4.80 1 -0.60 0.55 0.35 1.3502544
## 858 4.93 0 -0.31 0.73 0.42 -0.8554281
## 859 4.96 1 -0.25 0.78 0.44 1.1307602
## 860 5.14 1 0.15 1.17 0.54 0.9261787
## 861 4.54 0 -1.18 0.31 0.24 -0.5551190
## 862 4.76 1 -0.69 0.50 0.33 1.4114870
## 863 5.22 0 0.33 1.39 0.58 -1.1798526
## 864 6.04 1 2.15 8.58 0.90 0.3414414
## 865 4.35 1 -1.60 0.20 0.17 2.2238471
## 866 4.90 0 -0.38 0.68 0.41 -0.8274420
## 867 4.74 1 -0.73 0.48 0.32 1.4431368
## 868 4.52 1 -1.22 0.29 0.23 1.8418086
## 869 4.74 1 -0.73 0.48 0.32 1.4431368
## 870 5.01 0 -0.13 0.87 0.47 -0.9347728
## 871 3.52 1 -3.44 0.03 0.03 5.5818263
## 872 5.55 0 1.06 2.89 0.74 -1.7011134
## 873 5.34 1 0.60 1.82 0.64 0.7419735
## 874 4.75 0 -0.71 0.49 0.33 -0.7006608
## 875 5.71 1 1.42 4.13 0.80 0.4922907
## 876 4.77 1 -0.67 0.51 0.34 1.3959233
## 877 5.26 1 0.42 1.52 0.60 0.8107948
## 878 5.55 0 1.06 2.89 0.74 -1.7011134
## 879 5.11 0 0.09 1.09 0.52 -1.0443817
## 880 4.57 0 -1.11 0.33 0.25 -0.5738945
## 881 5.16 0 0.20 1.22 0.55 -1.1039156
## 882 5.40 1 0.73 2.07 0.67 0.6942190
## 883 5.82 1 1.66 5.27 0.84 0.4357658
## 884 4.66 0 -0.91 0.40 0.29 -0.6341178
## 885 5.22 1 0.33 1.39 0.58 0.8475635
## 886 5.30 0 0.51 1.66 0.62 -1.2892891
## 887 4.73 0 -0.76 0.47 0.32 -0.6852944
## 888 5.40 0 0.73 2.07 0.67 -1.4404676
## 889 5.40 1 0.73 2.07 0.67 0.6942190
## 890 5.59 1 1.15 3.16 0.76 0.5623484
## 891 5.08 1 0.02 1.02 0.51 0.9898894
## 892 4.96 0 -0.25 0.78 0.44 -0.8843608
## 893 6.06 1 2.19 8.97 0.90 0.3339531
## 894 4.61 0 -1.02 0.36 0.26 -0.5999200
## 895 5.76 1 1.53 4.61 0.82 0.4657416
## 896 5.64 1 1.26 3.53 0.78 0.5320211
## 897 5.64 1 1.26 3.53 0.78 0.5320211
## 898 4.67 0 -0.89 0.41 0.29 -0.6411878
## 899 5.07 1 0.00 1.00 0.50 1.0009261
## 900 5.41 1 0.75 2.12 0.68 0.6865643
## 901 5.14 1 0.15 1.17 0.54 0.9261787
## 902 4.90 0 -0.38 0.68 0.41 -0.8274420
## 903 4.98 1 -0.20 0.82 0.45 1.1059612
## 904 5.47 0 0.89 2.42 0.71 -1.5567208
## 905 5.49 1 0.93 2.53 0.72 0.6282878
## 906 5.50 1 0.95 2.59 0.72 0.6213601
## 907 5.58 1 1.13 3.09 0.76 0.5686182
## 908 6.18 1 2.46 11.70 0.92 0.2923491
## 909 5.37 1 0.66 1.94 0.66 0.7176992
## 910 5.73 1 1.46 4.31 0.81 0.4814941
## 911 5.30 1 0.51 1.66 0.62 0.7756212
## 912 4.85 1 -0.49 0.61 0.38 1.2774356
## 913 4.69 0 -0.84 0.43 0.30 -0.6555652
## 914 4.95 1 -0.27 0.76 0.43 1.1433675
## 915 5.15 1 0.18 1.19 0.54 0.9159662
## 916 5.43 1 0.80 2.22 0.69 0.6715070
## 917 5.08 1 0.02 1.02 0.51 0.9898894
ggplot(abejas_ajuste) + aes(ancho, pred_logit) + geom_line() + geom_point()
ggplot(abejas_ajuste) + aes(ancho, pred_odds) + geom_line() + geom_point()
ggplot(abejas_ajuste) + aes(ancho, pred_prob) + geom_line() + geom_point()
tab_model te da también el IC del OR:
sjPlot::plot_model(abejas_m1, show.values = T)
## Waiting for profiling to be done...
# tab_model(m1)
sjPlot::plot_model(abejas_m1, type = "pred")
## Data were 'prettified'. Consider using `terms="ancho [all]"` to get smooth plots.
## $ancho
Predecimos las probabilidades de que ciertos anchos observados tengan parásitos:
newdata <- data.frame(ancho = c(5, 6))
predict(abejas_m1, newdata, type = "response")
## 1 2
## 0.4608108 0.8870000
((abejas_m1$null.deviance-abejas_m1$deviance)/abejas_m1$null.deviance)*100
## [1] 12.28365
ATENCIÓN!!! Clase 15, página 15: los gráficos de residuos vs. predichos no son informativos porque la variable toma sólo dos valores! Lo mismo observados vs. predichos.
Hay que hacer el test de Hosmer-Lemeshow (más abajo).
ggplot(abejas_ajuste) +
aes(pred_prob, residuos_pearson) +
geom_point()
ggplot(abejas_ajuste) +
aes(pred_prob, parasitos) +
geom_point()
ggplot(abejas_ajuste) +
aes(ancho, parasitos) +
geom_jitter(height = .1, alpha = .5) +
geom_smooth(
method = "glm",
method.args = list(family = "binomial"),
se = F
) +
labs(
y = "Probabilidad estimada de presencia de parásitos",
x = "Ancho de celda (mm)"
)
En el test de Hosmer Lemeshow, la hipótesis nula es de mal ajuste, de modo que un p valor alto es bueno!
observados <- abejas$parasitos
esperados <- abejas_ajuste$pred_prob
ResourceSelection::hoslem.test(observados, esperados, g = 10) # Cuantiles, así lo usa Adriana en la clase
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: observados, esperados
## X-squared = 3.7642, df = 8, p-value = 0.8777
“Como la variables es dicotómica, no hay posibilidad de sobre o subdispersión”, pero esto para Bernoulli. Para Binomial sí va la prueba, parece (clase 16, p.30).
acaricidas <-
read_csv2("Clases 15 y 16 Regresión Bernoulli y binomial/acaricidas.csv") %>%
mutate(
total = as.integer(total),
muertos = as.integer(muertos),
acaricida = factor(acaricida),
vivos = total - muertos,
prop = muertos/total, # NO redondear! glm espera prop exacta
logit = log(prop/(1-prop))
)
## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
## Parsed with column specification:
## cols(
## ID = col_double(),
## concentracion = col_double(),
## acaricida = col_double(),
## muertos = col_double(),
## total = col_double()
## )
acaricidas
## # A tibble: 97 x 8
## ID concentracion acaricida muertos total vivos prop logit
## <dbl> <dbl> <fct> <int> <int> <int> <dbl> <dbl>
## 1 1 0 1 0 4 4 0 -Inf
## 2 2 0 1 0 4 4 0 -Inf
## 3 3 0 1 0 4 4 0 -Inf
## 4 4 0 1 0 4 4 0 -Inf
## 5 5 0 1 0 4 4 0 -Inf
## 6 6 0.25 1 2 4 2 0.5 0
## 7 7 0.25 1 1 4 3 0.25 -1.10
## 8 8 0.25 1 3 4 1 0.75 1.10
## 9 9 0.25 1 2 4 2 0.5 0
## 10 10 0.25 1 1 4 3 0.25 -1.10
## # … with 87 more rows
acaricidas %$%
psych::describeBy(prop, list(concentracion, acaricida), mat = T) %>%
mutate_if(is.double, round, 2)
## item group1 group2 vars n mean sd median trimmed mad min max
## 1 1 0 1 1 5 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 2 2 0.25 1 1 5 0.45 0.21 0.50 0.45 0.37 0.25 0.75
## 3 3 0.5 1 1 5 0.80 0.21 0.75 0.80 0.37 0.50 1.00
## 4 4 1 1 1 5 0.80 0.21 0.75 0.80 0.37 0.50 1.00
## 5 5 2 1 1 5 1.00 0.00 1.00 1.00 0.00 1.00 1.00
## 6 6 0 2 1 5 0.04 0.09 0.00 0.04 0.00 0.00 0.20
## 7 7 0.25 2 1 5 0.72 0.18 0.60 0.72 0.00 0.60 1.00
## 8 8 0.5 2 1 5 0.92 0.11 1.00 0.92 0.00 0.80 1.00
## 9 9 1 2 1 5 1.00 0.00 1.00 1.00 0.00 1.00 1.00
## 10 10 2 2 NA NA NA NA NA NA NA NA NA
## 11 11 0 3 1 8 0.12 0.17 0.00 0.12 0.00 0.00 0.33
## 12 12 0.25 3 1 8 0.54 0.25 0.50 0.54 0.25 0.33 1.00
## 13 13 0.5 3 1 8 0.88 0.17 1.00 0.88 0.00 0.67 1.00
## 14 14 1 3 1 8 0.96 0.12 1.00 0.96 0.00 0.67 1.00
## 15 15 2 3 NA NA NA NA NA NA NA NA NA
## 16 16 0 4 1 5 0.04 0.09 0.00 0.04 0.00 0.00 0.20
## 17 17 0.25 4 NA NA NA NA NA NA NA NA NA
## 18 18 0.5 4 1 5 0.48 0.23 0.40 0.48 0.30 0.20 0.80
## 19 19 1 4 1 5 0.80 0.14 0.80 0.80 0.00 0.60 1.00
## 20 20 2 4 1 5 0.92 0.11 1.00 0.92 0.00 0.80 1.00
## range skew kurtosis se
## 1 0.00 NaN NaN 0.00
## 2 0.50 0.25 -1.82 0.09
## 3 0.50 -0.25 -1.82 0.09
## 4 0.50 -0.25 -1.82 0.09
## 5 0.00 NaN NaN 0.00
## 6 0.20 1.07 -0.92 0.04
## 7 0.40 0.60 -1.67 0.08
## 8 0.20 -0.29 -2.25 0.05
## 9 0.00 NaN NaN 0.00
## 10 NA NA NA NA
## 11 0.33 0.42 -2.03 0.06
## 12 0.67 0.54 -1.27 0.09
## 13 0.33 -0.42 -2.03 0.06
## 14 0.33 -1.86 1.70 0.04
## 15 NA NA NA NA
## 16 0.20 1.07 -0.92 0.04
## 17 NA NA NA NA
## 18 0.60 0.19 -1.75 0.10
## 19 0.40 0.00 -1.40 0.06
## 20 0.20 -0.29 -2.25 0.05
Se va a ver lineal en la escala del p.l., es decir, en log odds o logit. Igual acá lo forzó a mano con el “lm”.
fig1 <-
acaricidas %>%
ggplot() +
aes(concentracion, prop, color = acaricida) +
geom_point() +
geom_smooth(se = F, method = "glm", method.args = list(family = "binomial"))
fig2 <-
ggplot(acaricidas) +
aes(concentracion, logit, color = acaricida) +
geom_point() +
geom_smooth(se = F, method = "lm")
fig1 + fig2
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning: Removed 53 rows containing non-finite values (stat_smooth).
Hay dos maneras de hacer la regresión binomial:
Atención: estos p-valores no importan! Primero miramos el p-value de anova
o drop1
y chequeamos los supuestos.
# acaricida_glm <- glm(cbind(muertos, vivos) ~ acaricida * concentracion, data = acaricidas, family = binomial)
# summary(acaricida_glm)
acaricida_glm <- glm(prop ~ acaricida * concentracion, family = binomial, weights = total, data = acaricidas)
summary(acaricida_glm)
##
## Call:
## glm(formula = prop ~ acaricida * concentracion, family = binomial,
## data = acaricidas, weights = total)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3707 -0.9532 0.1052 0.5776 2.1640
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6167 0.4368 -3.701 0.000215 ***
## acaricida2 -0.7360 0.7674 -0.959 0.337553
## acaricida3 0.0751 0.6250 0.120 0.904364
## acaricida4 -0.2032 0.6244 -0.325 0.744875
## concentracion 4.0994 0.9272 4.421 9.81e-06 ***
## acaricida2:concentracion 7.3729 2.6309 2.802 0.005072 **
## acaricida3:concentracion 2.1346 1.6356 1.305 0.191855
## acaricida4:concentracion -1.3016 1.0923 -1.192 0.233391
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 313.105 on 96 degrees of freedom
## Residual deviance: 91.451 on 89 degrees of freedom
## AIC: 186.38
##
## Number of Fisher Scoring iterations: 6
Supuestos de la regresión binomial (clase 16. p.30):
acaricidas_fit <-
augment(acaricida_glm, data = acaricidas, type.residual = "pearson") %>%
mutate(
.fitted = predict(acaricida_glm, type = "response")
)
fig_resid <-
acaricidas_fit %>%
ggplot() +
aes(x = .fitted, y = .resid, color = acaricida) +
geom_hline(yintercept = 0) +
geom_point() +
geom_label_repel(
mapping = aes(label = ID),
data = filter(acaricidas_fit, abs(.resid) > 2),
alpha = .6
)
fig_cook <-
ggplot(acaricidas_fit) +
aes(ID, .cooksd, color = acaricida) +
geom_point() +
geom_label_repel(
data = filter(acaricidas_fit, abs(.cooksd) > .05),
mapping = aes(label = ID),
alpha = .6
)
fig_resid + fig_cook
Acá deciden sacar outliers mirando los residuos estandarizados, no pasan por regresión robusta. Sacamos el 72.
acaricidas <- filter(acaricidas, ID != 72)
drop1
primero quita las interacciones. Acá se ve que hay interacción significativa:
acaricida_glm <- glm(prop ~ acaricida * concentracion, family = binomial, weights = total, data = acaricidas)
# car::Anova(acaricida_glm) # igual resultado
drop1(acaricida_glm, test = "LR")
## Single term deletions
##
## Model:
## prop ~ acaricida * concentracion
## Df Deviance AIC LRT Pr(>Chi)
## <none> 84.626 177.94
## acaricida:concentracion 3 110.687 198.00 26.061 9.262e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# calcular_sobredispersion(acaricida_glm)
ggplot(acaricidas_fit) +
aes(x = concentracion, y = prop, color = acaricida) +
geom_point() +
geom_smooth(
method = "glm",
method.args = list(family = binomial),
se = F
)
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!
emmeans::emmeans(acaricida_glm, pairwise ~ acaricida | concentracion) # para la dosis promedio, en escala logit
## $emmeans
## concentracion = 0.604:
## acaricida emmean SE df asymp.LCL asymp.UCL
## 1 0.86 0.348 Inf 0.178 1.542
## 2 4.58 1.015 Inf 2.589 6.568
## 3 2.90 0.678 Inf 1.569 4.225
## 4 -0.13 0.263 Inf -0.645 0.386
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## concentracion = 0.604:
## contrast estimate SE df z.ratio p.value
## 1 - 2 -3.72 1.073 Inf -3.465 0.0030
## 1 - 3 -2.04 0.762 Inf -2.674 0.0376
## 1 - 4 0.99 0.436 Inf 2.269 0.1054
## 2 - 3 1.68 1.220 Inf 1.378 0.5134
## 2 - 4 4.71 1.048 Inf 4.490 <.0001
## 3 - 4 3.03 0.727 Inf 4.165 0.0002
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimates
emmeans::emmeans(acaricida_glm, pairwise ~ acaricida | concentracion, type = "response") # para la dosis promedio, en odds ratio
## $emmeans
## concentracion = 0.604:
## acaricida prob SE df asymp.LCL asymp.UCL
## 1 0.703 0.0727 Inf 0.544 0.824
## 2 0.990 0.0102 Inf 0.930 0.999
## 3 0.948 0.0336 Inf 0.828 0.986
## 4 0.468 0.0655 Inf 0.344 0.595
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## concentracion = 0.604166666666667:
## contrast odds.ratio SE df z.ratio p.value
## 1 / 2 0.0243 0.0260 Inf -3.465 0.0030
## 1 / 3 0.1304 0.0993 Inf -2.674 0.0376
## 1 / 4 2.6901 1.1735 Inf 2.269 0.1054
## 2 / 3 5.3725 6.5563 Inf 1.378 0.5134
## 2 / 4 110.8357 116.2088 Inf 4.490 <.0001
## 3 / 4 20.6304 14.9938 Inf 4.165 0.0002
##
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale
emmeans::emtrends(acaricida_glm, pairwise ~ acaricida, var = "concentracion")
## $emtrends
## acaricida concentracion.trend SE df asymp.LCL asymp.UCL
## 1 4.10 0.927 Inf 2.28 5.92
## 2 11.47 2.462 Inf 6.65 16.30
## 3 7.91 1.693 Inf 4.59 11.23
## 4 2.80 0.577 Inf 1.67 3.93
##
## Trends are based on the logit (transformed) scale
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## 1 - 2 -7.37 2.63 Inf -2.802 0.0261
## 1 - 3 -3.81 1.93 Inf -1.976 0.1972
## 1 - 4 1.30 1.09 Inf 1.192 0.6321
## 2 - 3 3.56 2.99 Inf 1.191 0.6326
## 2 - 4 8.67 2.53 Inf 3.430 0.0034
## 3 - 4 5.12 1.79 Inf 2.859 0.0221
##
## P value adjustment: tukey method for comparing a family of 4 estimates
#calcula una pendiente para cada acaricida y las compara
emmeans::emmip(acaricida_glm, acaricida ~ concentracion, cov.reduce = range)
Se calcula igual que siempre.
m <- acaricida_glm
(m$null.deviance - m$deviance)/m$null.deviance
## [1] 0.7296899
The Hosmer–Lemeshow test is a statistical test for goodness of fit for logistic regression models. It is used frequently in risk prediction models. The test assesses whether or not the observed event rates match expected event rates in subgroups of the model population. The Hosmer–Lemeshow test specifically identifies subgroups as the deciles of fitted risk values. Models for which expected and observed event rates in subgroups are similar are called well calibrated.
acaricidas_fit %$%
ResourceSelection::hoslem.test(prop, .fitted, g = 10)
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: prop, .fitted
## X-squared = 5.9022, df = 8, p-value = 0.6582
nuevos_datos <-tibble(
acaricida = c("1", "1", "2", "2"),
concentracion = c(.5, .99, .5, .99)
)
nuevos_datos %>% predict(acaricida_glm, newdata = ., type = "response")
## 1 2 3 4
## 0.6065893 0.9199512 0.9671834 0.9998772
# sjPlot::plot_model(m1, show.values = TRUE)
# sjPlot::plot_model(m1, type="pred")
fig1 <-
sjPlot::plot_model(acaricida_glm, show.values = T) +
geom_hline(yintercept = 1, color = "Red2")
## Waiting for profiling to be done...
fig2 <-
sjPlot::plot_model(acaricida_glm, type = "pred")
fig1 + fig2
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#model-definition
La VR es un conteo (de polen en este caso), se modela con Poisson, pero hay conteos en diferentes colmenas. Las colmenas individuales no nos interesan, pero sí tenemos que dar cuenta de la variabilidad asociada al cambio de colmenas entre un tratamiento y otro, para no confundir esa variabilidad con un efecto del cambio de nivel del tratamiento.
El modelo es:
$$ _i = ( E[Y_{ij}] ) = _0 + 1 P{ij} + 2 A{ij} + _j \
j NID(0, ^2{colmena}) \
E[Y_i] = e^{_i} \ E[Y_i|colmena=j] = e^{_i + _j} $$
Acá entra GLMM: modelo lineal generalizado (porque usamos Poisson para la VR en lugar de Normal) mixto (porque tenemos efectos fijos, asociados a cada nivel del tratamiento, y efectos aleatorios, asociados a las diferentes colmenas).
kiwi <-
read_csv("Clase 17 MLGM/kiwi.csv") %>%
mutate(
# El orden de los niveles del factor es importante porque luego el primer nivel
# es el que está asociado al Intercept en la regresión, me interesa poner
# el control como tratamiento "base":
tratamiento = factor(tratamiento, levels = c("control", "proteina", "almibar")),
colmena = factor(colmena)
)
## Parsed with column specification:
## cols(
## tratamiento = col_character(),
## colmena = col_double(),
## muestreo = col_double(),
## polen = col_double()
## )
La regresión puede hacerse con lme4
o con glmmTBB
, pero ya no con stats::glm
:
Atenti al “Laplace approximation”, es el modo en el que se aproxima la máxima verosimilitud de los datos y parámetros:
kiwi_glmm <- lme4::glmer(polen ~ tratamiento + (1|colmena), data = kiwi, family = poisson,
nAGQ = 3) # nAGQ es el número de puntos de cuadratura para aproximas la integral de la verosimilitud
kiwi_glmm
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 3) [glmerMod]
## Family: poisson ( log )
## Formula: polen ~ tratamiento + (1 | colmena)
## Data: kiwi
## AIC BIC logLik deviance df.resid
## 867.6770 876.0544 -429.8385 859.6770 56
## Random effects:
## Groups Name Std.Dev.
## colmena (Intercept) 0.9127
## Number of obs: 60, groups: colmena, 15
## Fixed Effects:
## (Intercept) tratamientoproteina tratamientoalmibar
## 3.005 1.286 1.741
summary(kiwi_glmm)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 3) [glmerMod]
## Family: poisson ( log )
## Formula: polen ~ tratamiento + (1 | colmena)
## Data: kiwi
##
## AIC BIC logLik deviance df.resid
## 867.7 876.1 -429.8 859.7 56
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.2771 -2.4390 0.2772 1.9828 9.0443
##
## Random effects:
## Groups Name Variance Std.Dev.
## colmena (Intercept) 0.8331 0.9127
## Number of obs: 60, groups: colmena, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0048 0.4136 7.265 3.73e-13 ***
## tratamientoproteina 1.2856 0.5817 2.210 0.02710 *
## tratamientoalmibar 1.7414 0.5818 2.993 0.00276 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmntp
## tratmntprtn -0.711
## tratmntlmbr -0.711 0.505
kiwi_glmm_2 <- glmmTMB::glmmTMB(polen ~ tratamiento + (1|colmena), data = kiwi, family = poisson)
kiwi_glmm_2
## Formula: polen ~ tratamiento + (1 | colmena)
## Data: kiwi
## AIC BIC logLik df.resid
## 1206.8604 1215.2378 -599.4302 56
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## colmena (Intercept) 0.9127
##
## Number of obs: 60 / Conditional model: colmena, 15
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) tratamientoproteina tratamientoalmibar
## 3.005 1.286 1.741
summary(kiwi_glmm_2)
## Family: poisson ( log )
## Formula: polen ~ tratamiento + (1 | colmena)
## Data: kiwi
##
## AIC BIC logLik deviance df.resid
## 1206.9 1215.2 -599.4 1198.9 56
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## colmena (Intercept) 0.833 0.9127
## Number of obs: 60, groups: colmena, 15
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0047 0.4136 7.265 3.73e-13 ***
## tratamientoproteina 1.2857 0.5817 2.210 0.02709 *
## tratamientoalmibar 1.7414 0.5817 2.993 0.00276 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(glmmTMB)
fixef(kiwi_glmm)
## (Intercept) tratamientoproteina tratamientoalmibar
## 3.004780 1.285631 1.741392
ranef(kiwi_glmm)
## $colmena
## (Intercept)
## 1 0.61776739
## 2 0.92189120
## 3 -1.89439666
## 4 0.20569137
## 5 0.15874980
## 6 0.49166956
## 7 -0.14263967
## 8 -0.15053311
## 9 -0.07425801
## 10 -0.11549101
## 11 0.40891518
## 12 -1.31914603
## 13 1.82714070
## 14 0.40891518
## 15 -1.27614389
##
## with conditional variances for "colmena"
Las mediciones de una misma colmena no son repeticiones independientes del tratamiento, pero sí son repeticiones independientes de la variabilidad dentro de la colmena. CON UN CAVEAT: ojo con muestreos en orden, secuencia temporal que tenga una correlación. Hay que analizar correlación entre muestreos de la misma colmena.
Para considerarlas como repeticiones aleatorias, las correlaciones deberían permenecer constantes entre todos los pares de muestreos (simetría compuesta).
for (i in 1:4) {
kiwi[[glue("muestreo_{i}")]] <- filter(kiwi, muestreo == i) %>% pull(polen)
}
kiwi %>%
select(matches("muestreo_")) %>%
GGally::ggpairs(progress = F)
Modelos condicionales o mixtos: se incluyen VE de efecots aleatorios que inducen la correlación entre observaciones. - parámetros sujeto-específicos: proporcionan un modelo para cada individuo, condicional a un nivel dado de efecto aleatorio (i.e. | colmena = 1)
Modelos marginales: no incluyen VE de efectos aleatorios, pero se explicita una estructura de correlación entre las observaciones. - parámetros promedios poblacionales: ajustan un modelo general para la estructura promedio de la población de individuos.
Con errores de distribución Gaussiana, ambos enfoques son equivalentes, pero con otras distribuciones los parámetros son intrínsecamente diferentes, se interpretan como el sub-bullet.
Atenti a los grados de libertad! Ahora hay un gl menos por estimar la varianza de los efectos aleatorios.
e1 <- resid(kiwi_glmm, type = "pearson")
n <- nrow(kiwi)
k <- length(fixef(kiwi_glmm)) + 1 # parametros estimados; se suma varianza(colmenas)
dispersion <- sum(e1^2) / (n - k)
# Se pueden pedir así los GL también:
df.residual(kiwi_glmm) == (n - k)
## [1] TRUE
dispersion
## [1] 12.76808
calcular_sobredispersion <- function(modelo) {
grados_lib_residuos <- df.residual(modelo)
residuos_P <- residuals(modelo, type = "pearson")
chi_obs <- sum(residuos_P^2)
ratio <- chi_obs/grados_lib_residuos
# La distribución chi cuadrado para esto es aproximada, y sólo válida cuando
# lambda > 5 (cuál lambda??) en Poisson y min(np, n(1-p)) > 5 para
# binomial
p_valor <- pchisq(chi_obs, df = grados_lib_residuos, lower.tail = F)
return(list(
ratio = ratio,
chi_obs = chi_obs,
p_valor = p_valor
))
}
calcular_sobredispersion(kiwi_glmm)
## $ratio
## [1] 12.76808
##
## $chi_obs
## [1] 715.0125
##
## $p_valor
## [1] 4.706282e-115
Predichos vs. pearson, pearson por tratamiento, pearson por nivel del factor aleatorio
kiwi_fit <- augment(
kiwi_glmm, data = kiwi,
type.residual = "pearson",
# type.predict = "response" # No hizo lo que yo quería
) %>%
mutate(
.fitted = predict(kiwi_glmm, type = "response")
)
## Warning in hatvalues.merMod(model): the hat matrix may not make sense for
## GLMMs
kiwi_fit
## tratamiento colmena muestreo polen muestreo_1 muestreo_2 muestreo_3
## 1 almibar 1 1 236 236 167 227
## 2 almibar 1 2 167 306 308 241
## 3 almibar 1 3 227 36 10 2
## 4 almibar 1 4 225 249 208 66
## 5 almibar 2 1 306 139 145 146
## 6 almibar 2 2 308 147 140 153
## 7 almibar 2 3 241 129 43 42
## 8 almibar 2 4 304 35 89 100
## 9 almibar 3 1 36 111 87 48
## 10 almibar 3 2 10 98 76 27
## 11 almibar 3 3 2 32 50 9
## 12 almibar 3 4 19 13 1 0
## 13 almibar 4 1 249 130 122 127
## 14 almibar 4 2 208 31 51 11
## 15 almibar 4 3 66 11 2 1
## 16 almibar 4 4 43 236 167 227
## 17 almibar 5 1 139 306 308 241
## 18 almibar 5 2 145 36 10 2
## 19 almibar 5 3 146 249 208 66
## 20 almibar 5 4 110 139 145 146
## 21 proteina 6 1 147 147 140 153
## 22 proteina 6 2 140 129 43 42
## 23 proteina 6 3 153 35 89 100
## 24 proteina 6 4 38 111 87 48
## 25 proteina 7 1 129 98 76 27
## 26 proteina 7 2 43 32 50 9
## 27 proteina 7 3 42 13 1 0
## 28 proteina 7 4 39 130 122 127
## 29 proteina 8 1 35 31 51 11
## 30 proteina 8 2 89 11 2 1
## 31 proteina 8 3 100 236 167 227
## 32 proteina 8 4 27 306 308 241
## 33 proteina 9 1 111 36 10 2
## 34 proteina 9 2 87 249 208 66
## 35 proteina 9 3 48 139 145 146
## 36 proteina 9 4 25 147 140 153
## 37 proteina 10 1 98 129 43 42
## 38 proteina 10 2 76 35 89 100
## 39 proteina 10 3 27 111 87 48
## 40 proteina 10 4 59 98 76 27
## 41 control 11 1 32 32 50 9
## 42 control 11 2 50 13 1 0
## 43 control 11 3 9 130 122 127
## 44 control 11 4 31 31 51 11
## 45 control 12 1 13 11 2 1
## 46 control 12 2 1 236 167 227
## 47 control 12 3 0 306 308 241
## 48 control 12 4 6 36 10 2
## 49 control 13 1 130 249 208 66
## 50 control 13 2 122 139 145 146
## 51 control 13 3 127 147 140 153
## 52 control 13 4 125 129 43 42
## 53 control 14 1 31 35 89 100
## 54 control 14 2 51 111 87 48
## 55 control 14 3 11 98 76 27
## 56 control 14 4 29 32 50 9
## 57 control 15 1 11 13 1 0
## 58 control 15 2 2 130 122 127
## 59 control 15 3 1 31 51 11
## 60 control 15 4 7 11 2 1
## muestreo_4 .fitted .resid .cooksd .fixed .mu
## 1 225 213.564610 1.50944707 9.207808e-04 4.746172 213.564610
## 2 304 213.564610 -3.31411989 3.966441e-03 4.746172 213.564610
## 3 19 213.564610 0.90996473 3.302095e-04 4.746172 213.564610
## 4 43 213.564610 0.77567166 2.392164e-04 4.746172 213.564610
## 5 110 289.473344 0.96233140 2.718703e-04 4.746172 289.473344
## 6 38 289.473344 1.07759583 3.416535e-04 4.746172 289.473344
## 7 39 289.473344 -2.93467877 2.338825e-03 4.746172 289.473344
## 8 27 289.473344 0.84681441 2.100502e-04 4.746172 289.473344
## 9 25 17.318501 3.91440162 9.842389e-02 4.746172 17.318501
## 10 59 17.318501 -1.91133380 1.510502e-02 4.746172 17.318501
## 11 31 17.318501 -4.69068113 6.617739e-02 4.746172 17.318501
## 12 6 17.318501 0.39776855 7.973891e-04 4.746172 17.318501
## 13 125 141.438273 8.15721283 4.828349e-02 4.746172 141.438273
## 14 29 141.438273 5.22658757 1.848980e-02 4.746172 141.438273
## 15 7 141.438273 -7.08978169 2.375016e-02 4.746172 141.438273
## 16 225 141.438273 -9.72005680 4.044000e-02 4.746172 141.438273
## 17 304 134.952360 0.34670656 7.511008e-05 4.746172 134.952360
## 18 19 134.952360 0.85450303 4.628311e-04 4.746172 134.952360
## 19 43 134.952360 0.93844553 5.595430e-04 4.746172 134.952360
## 20 110 134.952360 -2.21981945 2.854419e-03 4.746172 134.952360
## 21 38 119.352452 2.44137273 4.481414e-03 4.290411 119.352452
## 22 39 119.352452 1.83905714 2.499417e-03 4.290411 119.352452
## 23 27 119.352452 2.94994127 6.637563e-03 4.290411 119.352452
## 24 25 119.352452 -8.70191380 3.880102e-02 4.290411 119.352452
## 25 59 63.292806 7.23132916 9.018154e-02 4.290411 63.292806
## 26 31 63.292806 -2.70932792 8.601536e-03 4.290411 63.292806
## 27 6 63.292806 -2.85254958 9.470166e-03 4.290411 63.292806
## 28 125 63.292806 -3.28895570 1.232671e-02 4.290411 63.292806
## 29 29 62.795174 -3.83056615 1.639460e-02 4.290411 62.795174
## 30 7 62.795174 3.10954944 1.457218e-02 4.290411 62.795174
## 31 225 62.795174 4.31841870 2.937383e-02 4.290411 62.795174
## 32 304 62.795174 -5.10020736 2.719011e-02 4.290411 62.795174
## 33 19 67.772285 4.80356480 3.403324e-02 4.290411 67.772285
## 34 43 67.772285 2.23649335 6.733395e-03 4.290411 67.772285
## 35 110 67.772285 -2.53557618 7.120203e-03 4.290411 67.772285
## 36 38 67.772285 -5.97333189 3.331989e-02 4.290411 67.772285
## 37 39 65.034658 3.79981015 2.149702e-02 4.290411 65.034658
## 38 27 65.034658 1.32396894 2.378524e-03 4.290411 65.034658
## 39 25 65.034658 -5.34778548 2.861687e-02 4.290411 65.034658
## 40 59 65.034658 -0.76035240 7.203914e-04 4.290411 65.034658
## 41 31 30.377286 0.29185552 2.398434e-04 3.004780 30.377286
## 42 6 30.377286 3.25382184 3.507212e-02 3.004780 30.377286
## 43 125 30.377286 -4.56706708 4.162450e-02 3.004780 30.377286
## 44 29 30.377286 0.11260049 3.532002e-05 3.004780 30.377286
## 45 7 5.395871 2.76657579 1.735295e-01 3.004780 5.395871
## 46 225 5.395871 -2.32819099 5.799146e-02 3.004780 5.395871
## 47 304 5.395871 -3.28507857 8.737704e-02 3.004780 5.395871
## 48 19 5.395871 0.25543560 1.095303e-03 3.004780 5.395871
## 49 43 125.451682 0.40366346 1.097670e-04 3.004780 125.451682
## 50 110 125.451682 -0.30960118 6.321667e-05 3.004780 125.451682
## 51 38 125.451682 0.13795333 1.272009e-05 3.004780 125.451682
## 52 39 125.451682 -0.04035112 1.082520e-06 3.004780 125.451682
## 53 27 30.377286 0.11260049 3.532002e-05 3.004780 30.377286
## 54 25 30.377286 3.40644686 3.873785e-02 3.004780 30.377286
## 55 59 30.377286 -4.05055240 3.420029e-02 3.004780 30.377286
## 56 31 30.377286 -0.25181538 1.727795e-04 3.004780 30.377286
## 57 6 5.632966 1.99740672 7.916695e-02 3.004780 5.632966
## 58 125 5.632966 -1.76747730 3.627420e-02 3.004780 5.632966
## 59 29 5.632966 -2.41011607 5.899204e-02 3.004780 5.632966
## 60 7 5.632966 0.55476881 5.136106e-03 3.004780 5.632966
## .offset .sqrtXwt .sqrtrwt .weights .wtres .eta
## 1 0 14.613850 0.06842824 1 1.53521417 5.363939
## 2 0 14.613850 0.06842824 1 -3.18633424 5.363939
## 3 0 14.613850 0.06842824 1 0.91936003 5.363939
## 4 0 14.613850 0.06842824 1 0.78250356 5.363939
## 5 0 17.013916 0.05877542 1 0.97136108 5.668063
## 6 0 17.013916 0.05877542 1 1.08891191 5.668063
## 7 0 17.013916 0.05877542 1 -2.84904096 5.668063
## 8 0 17.013916 0.05877542 1 0.85381025 5.668063
## 9 0 4.161550 0.24029506 1 4.48907198 2.851775
## 10 0 4.161550 0.24029506 1 -1.75859968 2.851775
## 11 0 4.161550 0.24029506 1 -3.68096019 2.851775
## 12 0 4.161550 0.24029506 1 0.40405589 2.851775
## 13 0 11.892782 0.08408461 1 9.04428617 4.951863
## 14 0 11.892782 0.08408461 1 5.59681705 4.951863
## 15 0 11.892782 0.08408461 1 -6.34319795 4.951863
## 16 0 11.892782 0.08408461 1 -8.27714404 4.951863
## 17 0 11.616900 0.08608149 1 0.34842689 4.904922
## 18 0 11.616900 0.08608149 1 0.86491581 4.904922
## 19 0 11.616900 0.08608149 1 0.95099729 4.904922
## 20 0 11.616900 0.08608149 1 -2.14793623 4.904922
## 21 0 10.924855 0.09153440 1 2.53070166 4.782081
## 22 0 10.924855 0.09153440 1 1.88996088 4.782081
## 23 0 10.924855 0.09153440 1 3.07990805 4.782081
## 24 0 10.924855 0.09153440 1 -7.44654767 4.782081
## 25 0 7.955678 0.12569640 1 8.25915754 4.147772
## 26 0 7.955678 0.12569640 1 -2.55073253 4.147772
## 27 0 7.955678 0.12569640 1 -2.67642893 4.147772
## 28 0 7.955678 0.12569640 1 -3.05351812 4.147772
## 29 0 7.924341 0.12619346 1 -3.50756936 4.139878
## 30 0 7.924341 0.12619346 1 3.30687774 4.139878
## 31 0 7.924341 0.12619346 1 4.69500585 4.139878
## 32 0 7.924341 0.12619346 1 -4.51711708 4.139878
## 33 0 8.232392 0.12147137 1 5.25092991 4.216153
## 34 0 8.232392 0.12147137 1 2.33561698 4.216153
## 35 0 8.232392 0.12147137 1 -2.40176653 4.216153
## 36 0 8.232392 0.12147137 1 -5.19560809 4.216153
## 37 0 8.064407 0.12400168 1 4.08775773 4.174920
## 38 0 8.064407 0.12400168 1 1.35972078 4.174920
## 39 0 8.064407 0.12400168 1 -4.71636153 4.174920
## 40 0 8.064407 0.12400168 1 -0.74830778 4.174920
## 41 0 5.511559 0.18143686 1 0.29442010 3.413695
## 42 0 5.511559 0.18143686 1 3.56028353 3.413695
## 43 0 5.511559 0.18143686 1 -3.87862761 3.413695
## 44 0 5.511559 0.18143686 1 0.11298324 3.413695
## 45 0 2.322901 0.43049612 1 3.27354817 1.685634
## 46 0 2.322901 0.43049612 1 -1.89240522 1.685634
## 47 0 2.322901 0.43049612 1 -2.32290133 1.685634
## 48 0 2.322901 0.43049612 1 0.26007536 1.685634
## 49 0 11.200522 0.08928156 1 0.40608090 4.831921
## 50 0 11.200522 0.08928156 1 -0.30817156 4.831921
## 51 0 11.200522 0.08928156 1 0.13823623 4.831921
## 52 0 11.200522 0.08928156 1 -0.04032689 4.831921
## 53 0 5.511559 0.18143686 1 0.11298324 3.413695
## 54 0 5.511559 0.18143686 1 3.74172038 3.413695
## 55 0 5.511559 0.18143686 1 -3.51575389 3.413695
## 56 0 5.511559 0.18143686 1 -0.24989047 3.413695
## 57 0 2.373387 0.42133879 1 2.26133968 1.728636
## 58 0 2.373387 0.42133879 1 -1.53070942 1.728636
## 59 0 2.373387 0.42133879 1 -1.95204821 1.728636
## 60 0 2.373387 0.42133879 1 0.57598452 1.728636
fig1 <- ggplot(kiwi_fit) +
aes(x = .fitted, y = .resid, color = tratamiento) +
geom_hline(yintercept = 0) +
geom_point() +
labs(
y = "Residuos Pearson",
x = "Predichos"
) +
guides(color = F)
fig2 <-
ggplot(kiwi_fit) +
aes(x = tratamiento, y = .resid) +
geom_hline(yintercept = 0) +
geom_boxplot(width = .25)
fig3 <-
ggplot(kiwi_fit) +
aes(x = colmena, y = .resid) +
geom_hline(yintercept = 0) +
geom_boxplot(width = .8)
fig4 <-
ggplot(kiwi_fit) +
aes(x = .fitted, y = polen, color = tratamiento) +
geom_abline() +
geom_point() +
theme(legend.position = "top")
fig4 + fig1 + fig2 + fig3 + plot_layout(nrow = 1)
$$ (0, 1) \
H_0, _j = 0: \
(0, 1) \ $$
Recordar que es asintótico, testea H0 si beta_j = 0, los p-valores están en el summary para cada beta:
# Test de Wald en los p-valores!
summary(kiwi_glmm)$coefficients
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.004780 0.4135917 7.265088 3.727955e-13
## tratamientoproteina 1.285631 0.5816981 2.210135 2.709580e-02
## tratamientoalmibar 1.741392 0.5817504 2.993366 2.759186e-03
drop1(kiwi_glmm, test = "Chisq")
## Single term deletions
##
## Model:
## polen ~ tratamiento + (1 | colmena)
## Df AIC LRT Pr(Chi)
## <none> 867.68
## tratamiento 2 871.16 7.4833 0.02371 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kiwi_glmm_sin_trat <- glmmTMB(polen ~ (1|colmena), data = kiwi, family = poisson)
AIC(kiwi_glmm, kiwi_glmm_sin_trat)
## df AIC
## kiwi_glmm 4 867.677
## kiwi_glmm_sin_trat 2 1210.344
Agregado por mí, no está en la teoríca. Atenti con efectos aleatorios, deviance(modelo)
devianza <- deviance(kiwi_glmm)
kiwi_nulo <- glm(polen ~ 1, data = kiwi, family = poisson)
devianza_nula <- deviance(kiwi_nulo)
devianza_explicada <- 1 - (devianza/devianza_nula)
100 * devianza_explicada
## [1] 83.57083
# nbinom1: variance increases linearly with the mean, i.e. variance=mu*(1+phi)
# nbinom2: variance increases quadratically with the mean, i.e. variance=mu*(1+mu/phi)
kiwi_glmm_nb <- glmmTMB::glmmTMB(polen ~ tratamiento + (1|colmena), data = kiwi, family = "nbinom2")
df.residual(kiwi_glmm)
## [1] 56
df.residual(kiwi_glmm_nb) # Tiene un grado de libertad menos!
## [1] 55
summary(kiwi_glmm_nb)
## Family: nbinom2 ( log )
## Formula: polen ~ tratamiento + (1 | colmena)
## Data: kiwi
##
## AIC BIC logLik deviance df.resid
## 626.6 637.0 -308.3 616.6 55
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## colmena (Intercept) 0.7418 0.8613
## Number of obs: 60, groups: colmena, 15
##
## Overdispersion parameter for nbinom2 family (): 3.08
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0490 0.4122 7.396 1.4e-13 ***
## tratamientoproteina 1.2784 0.5794 2.207 0.02735 *
## tratamientoalmibar 1.7393 0.5794 3.002 0.00269 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kiwi_glmm_cp <- glmmTMB::glmmTMB(polen ~ tratamiento + (1|colmena), data = kiwi, family = compois)
kiwi_glmm_cp
## Formula: polen ~ tratamiento + (1 | colmena)
## Data: kiwi
## AIC BIC logLik df.resid
## 613.7348 624.2066 -301.8674 55
## Random-effects (co)variances:
##
## Conditional model:
## Groups Name Std.Dev.
## colmena (Intercept) 0.8143
##
## Number of obs: 60 / Conditional model: colmena, 15
##
## Overdispersion parameter for compois family (): 22.9
##
## Fixed Effects:
##
## Conditional model:
## (Intercept) tratamientoproteina tratamientoalmibar
## 3.116 1.163 1.679
summary(kiwi_glmm_cp)
## Family: compois ( log )
## Formula: polen ~ tratamiento + (1 | colmena)
## Data: kiwi
##
## AIC BIC logLik deviance df.resid
## 613.7 624.2 -301.9 603.7 55
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## colmena (Intercept) 0.6632 0.8143
## Number of obs: 60, groups: colmena, 15
##
## Overdispersion parameter for compois family (): 22.9
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1162 0.4054 7.686 1.52e-14 ***
## tratamientoproteina 1.1634 0.5565 2.091 0.03655 *
## tratamientoalmibar 1.6786 0.5528 3.036 0.00239 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Luego hacer los mismos plots que arriba.
Ojo, no sirve para subdispersión. Tampoco sirve para “sobredispersión causada por exceso de ceros” o “modelos alterados por el cero”.
Overdispersion is common in models of count data in ecology and evolutionary biology, and can occur due to missing covariates, non-independent (aggregated) data, or an excess frequency of zeroes (zero-inflation). Accounting for overdispersion in such models is vital, as failing to do so can lead to biased parameter estimates, and false conclusions regarding hypotheses of interest. Observation-level random effects (OLRE), where each data point receives a unique level of a random effect that models the extra-Poisson variation present in the data, are commonly employed to cope with overdispersion in count data.
Agregamos el término + (1|obs_id)
al modelo. Tiene que ser una librería que admita random effects, obviamente, glmmTMB
o glmer
:
kiwi %<>% mutate(obs_id = 1:n())
kiwi_glmm_olre <- glmmTMB::glmmTMB(polen ~ tratamiento + (1|colmena) + (1|obs_id), data = kiwi, family = poisson)
kiwi_glmm_olre <- glmer(polen ~ tratamiento + (1|colmena) + (1|obs_id), data = kiwi, family = poisson)
Ahora calcular la sobredispersión no tiene sentido. En el script de clase la calculan igual pero no lo consideran:
e <- resid(kiwi_glmm_olre, type = "pearson")
k <- length(fixef(kiwi_glmm_olre)) + length(ranef(kiwi_glmm_olre))
N <- length(e)
grados_lib <- N - k
sum(e^2)/grados_lib
## [1] 0.1854162
Cuando encontramos sobredispersión y probamos diferentes modos de solucionarla, como CMP, BN y OLRE, podemos comparar las regresiones resultantes con AIC o con Anova:
AIC(
kiwi_glmm_cp,
kiwi_glmm_nb,
kiwi_glmm_olre
)
## df AIC
## kiwi_glmm_cp 5 613.7348
## kiwi_glmm_nb 5 626.5751
## kiwi_glmm_olre 5 630.3457
bbmle::AICtab(
kiwi_glmm_cp,
kiwi_glmm_nb,
kiwi_glmm_olre
)
## dAIC df
## kiwi_glmm_cp 0.0 5
## kiwi_glmm_nb 12.8 5
## kiwi_glmm_olre 16.6 5
Con anova no pudo comparar también el OLRE
anova(
kiwi_glmm_cp,
kiwi_glmm_nb
)
## Data: kiwi
## Models:
## kiwi_glmm_cp: polen ~ tratamiento + (1 | colmena), zi=~0, disp=~1
## kiwi_glmm_nb: polen ~ tratamiento + (1 | colmena), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## kiwi_glmm_cp 5 613.73 624.21 -301.87 603.73
## kiwi_glmm_nb 5 626.58 637.05 -308.29 616.58 0 0 1
Selecciona el de CMP y ahora se pone a seleccionar variables dentro del modelo y distribución elegido.
drop1(kiwi_glmm_cp, test = "Chi")
## Single term deletions
##
## Model:
## polen ~ tratamiento + (1 | colmena)
## Df AIC LRT Pr(>Chi)
## <none> 613.73
## tratamiento 2 617.28 7.5477 0.02296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No quita la variable tratamiento. Listo, una vez que nos quedamos con un modelo con ciertas variables, empezamos a comparar los efectos y ver las estimaciones con sus IC.
library(emmeans)
comp_PL <- emmeans(kiwi_glmm_cp, pairwise ~ tratamiento)
summary(comp_PL)
## $emmeans
## tratamiento emmean SE df lower.CL upper.CL
## control 3.12 0.405 55 2.30 3.93
## proteina 4.28 0.381 55 3.52 5.04
## almibar 4.79 0.379 55 4.04 5.55
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## control - proteina -1.163 0.556 55 -2.091 0.1011
## control - almibar -1.679 0.553 55 -3.036 0.0101
## proteina - almibar -0.515 0.537 55 -0.959 0.6058
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
Acá se ve que los IC de los contrastes en la escala del PL (diferencias de log) contienen al cero en dos casos: significa que esas diferencias no son significativas a nivel alfa:
confint(comp_PL)
## $emmeans
## tratamiento emmean SE df lower.CL upper.CL
## control 3.12 0.405 55 2.30 3.93
## proteina 4.28 0.381 55 3.52 5.04
## almibar 4.79 0.379 55 4.04 5.55
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL
## control - proteina -1.163 0.556 55 -2.50 0.177
## control - almibar -1.679 0.553 55 -3.01 -0.347
## proteina - almibar -0.515 0.537 55 -1.81 0.779
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
Atenti al $cond
de abajo. La magnitud del efecto descripta es condicional a estar en una colmena específica de la que se describió su alfa_j (la varianza que aporta).
“Los modelos condicionales hablan del sujeto típico, no del promedio de los sujetos” ~.
summary(kiwi_glmm_cp)$coefficients$cond[, 1:2] %>% exp %>% round(2)
## Estimate Std. Error
## (Intercept) 22.56 1.50
## tratamientoproteina 3.20 1.74
## tratamientoalmibar 5.36 1.74
Lo de arriba significa que:
plot(comp_PL, comparisons = TRUE)
comp_VR <- emmeans(kiwi_glmm_cp, pairwise ~ tratamiento, type = "response")
confint(comp_VR)
## $emmeans
## tratamiento rate SE df lower.CL upper.CL
## control 22.6 9.15 55 10.0 50.8
## proteina 72.2 27.54 55 33.6 155.1
## almibar 120.9 45.77 55 56.6 258.2
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df lower.CL upper.CL
## control / proteina 0.312 0.174 55 0.0818 1.194
## control / almibar 0.187 0.103 55 0.0493 0.707
## proteina / almibar 0.597 0.321 55 0.1637 2.180
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
## Intervals are back-transformed from the log scale
kiwi$tratamiento <- factor(kiwi$tratamiento, levels = c("almibar", "proteina", "control"))
comp_VR <- emmeans(kiwi_glmm_cp, pairwise ~ tratamiento, type = "response")
confint(comp_VR)
## $emmeans
## tratamiento rate SE df lower.CL upper.CL
## almibar 22.6 9.15 55 10.0 50.8
## proteina 72.2 27.54 55 33.6 155.1
## control 120.9 45.77 55 56.6 258.2
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df lower.CL upper.CL
## almibar / proteina 0.312 0.174 55 0.0818 1.194
## almibar / control 0.187 0.103 55 0.0493 0.707
## proteina / control 0.597 0.321 55 0.1637 2.180
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimates
## Intervals are back-transformed from the log scale
Lo de las flechitas rojas no termina de cerrarme, se supone que si se solapan los tratamientos no son significativamente diferentes:
plot(comp_VR$emmeans, comparison = T)
plot(comp_VR$contrasts, comparison = T)
El gráfico que sigue es bueno para presentar a gente que pide el analisis y que no necesariamente sabe de estadística. Tiene sentido intuitivo en la escala de la VR y el otro es en la escala del PL, menos sentido intuitivo, pero bue, nos lo muestran:
estimados <- as.data.frame(confint(comp_PL$emmeans))
ggplot(estimados) +
aes(tratamiento, emmean) +
geom_col(
width = .4,
fill = "Grey80"
) +
geom_errorbar(
aes(ymin = emmean - SE, ymax = emmean + SE),
width = .2
) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
color = "RoyalBlue",
position = position_nudge(x = .3),
width = .1
) +
geom_point(color = "Tomato", size = 4) +
geom_hline(yintercept = 0) +
geom_text(
aes(label = round(emmean, 1)),
position = position_nudge(x = -.3),
color = "Tomato"
) +
geom_text(
aes(label = c("A", "AB", "B"), y = max(emmean + SE) + .5)
) +
labs(
title = "Escala del PL, estimados y error estándar",
subtitle = "En azul: IC nivel 95%",
y = "log granos de polen"
)
estimados <- as.data.frame(confint(comp_VR$emmeans))
ggplot(estimados) +
aes(tratamiento, rate) +
geom_col(
width = .4,
fill = "Grey80"
) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
color = "RoyalBlue",
width = .1
) +
geom_point(color = "Tomato", size = 4) +
geom_hline(yintercept = 0) +
geom_text(
aes(label = round(rate, 1)),
position = position_nudge(x = -.3),
color = "Tomato"
) +
geom_text(
aes(label = c("A", "AB", "B"), y = max(rate + SE) * 1.1)
) +
labs(
title = "Escala de la VR, estimados e IC nivel 95%",
y = "Granos de polen"
)
library(sjPlot)
plot_model(kiwi_glmm_nb, type="pred")
## $tratamiento
# tab_model(m3)
ranef(kiwi_glmm_cp)
## $colmena
## (Intercept)
## 1 0.54879533
## 2 0.84958475
## 3 -1.59551794
## 4 0.14911174
## 5 0.10435745
## 6 0.47254354
## 7 -0.11943591
## 8 -0.12651859
## 9 -0.05775981
## 10 -0.09501747
## 11 0.25619691
## 12 -1.08102838
## 13 1.61553363
## 14 0.25619691
## 15 -1.05165815
El modelo propone efectos aleatorios que son iid ~ N(0, sigma^2_colmena). Testeamos esto:
shapiro.test(unlist(ranef(kiwi_glmm_cp)))
##
## Shapiro-Wilk normality test
##
## data: unlist(ranef(kiwi_glmm_cp))
## W = 0.94309, p-value = 0.4228
Los modelos que vimos son condicionales y se interpretan como sujeto-específicos: para un sujeto perteneciente a tal o cual colmena, se predice el valor esperado, pero NO para el promedio de la población sin importar la colmena.
Los modelos condicionales son para tal o cual nivel de la variable aleatoria (tipo tal o cual colmena en el ejemplo anterior). NO hablan del promedio de los sujetos! En eso difieren de los modelos marginales.
Los modelos marginales son una alternativa a los modelos mixtos (diapo 25, clase 17):
Atenti a los IDs repetidos, son el mismo paciente! Va a haber correlaciones acá! Modelamos con efecto aleatorio a nivel id
para dar cuenta de eso:
escle <- read.csv2("Clase 18 MLGM DMR/esclerosis.csv")
escle
## id tiempo tratpre edad mejora grupo
## 1 1 0 0 47 1 AZ
## 2 1 12 0 47 1 AZ
## 3 1 15 0 47 1 AZ
## 4 1 18 0 47 0 AZ
## 5 2 0 0 48 1 AZ
## 6 2 6 0 48 0 AZ
## 7 2 9 0 48 1 AZ
## 8 2 12 0 48 0 AZ
## 9 2 18 0 48 1 AZ
## 10 3 0 0 53 0 AZ
## 11 3 3 0 53 0 AZ
## 12 3 6 0 53 1 AZ
## 13 3 12 0 53 0 AZ
## 14 3 15 0 53 1 AZ
## 15 3 18 0 53 1 AZ
## 16 4 0 0 46 1 AZ
## 17 4 6 0 46 0 AZ
## 18 4 9 0 46 1 AZ
## 19 4 12 0 46 1 AZ
## 20 4 15 0 46 1 AZ
## 21 4 18 0 46 1 AZ
## 22 5 0 0 56 0 AZ
## 23 5 3 0 56 1 AZ
## 24 5 6 0 56 0 AZ
## 25 5 9 0 56 0 AZ
## 26 5 12 0 56 0 AZ
## 27 5 18 0 56 1 AZ
## 28 6 3 0 50 1 AZ
## 29 6 6 0 50 0 AZ
## 30 6 9 0 50 0 AZ
## 31 6 12 0 50 1 AZ
## 32 6 15 0 50 1 AZ
## 33 6 18 0 50 1 AZ
## 34 7 0 0 52 0 AZ
## 35 7 3 0 52 0 AZ
## 36 7 6 0 52 0 AZ
## 37 7 9 0 52 1 AZ
## 38 7 12 0 52 1 AZ
## 39 7 15 0 52 1 AZ
## 40 7 18 0 52 1 AZ
## 41 8 0 0 44 0 AZ
## 42 8 3 0 44 1 AZ
## 43 8 6 0 44 0 AZ
## 44 8 9 0 44 1 AZ
## 45 8 12 0 44 1 AZ
## 46 8 15 0 44 1 AZ
## 47 8 18 0 44 1 AZ
## 48 9 0 0 44 0 AZ
## 49 9 3 0 44 0 AZ
## 50 9 6 0 44 0 AZ
## 51 9 9 0 44 0 AZ
## 52 9 12 0 44 0 AZ
## 53 9 18 0 44 1 AZ
## 54 10 6 0 53 0 AZ
## 55 10 9 0 53 0 AZ
## 56 10 12 0 53 1 AZ
## 57 10 15 0 53 1 AZ
## 58 10 18 0 53 1 AZ
## 59 11 0 0 35 0 AZ
## 60 11 3 0 35 0 AZ
## 61 11 6 0 35 1 AZ
## 62 11 9 0 35 1 AZ
## 63 11 12 0 35 1 AZ
## 64 11 15 0 35 1 AZ
## 65 11 18 0 35 0 AZ
## 66 12 0 0 61 0 AZ
## 67 12 6 0 61 1 AZ
## 68 12 15 0 61 1 AZ
## 69 12 18 0 61 0 AZ
## 70 13 0 0 50 0 AZ
## 71 13 6 0 50 0 AZ
## 72 13 9 0 50 0 AZ
## 73 13 12 0 50 0 AZ
## 74 13 15 0 50 1 AZ
## 75 13 18 0 50 1 AZ
## 76 14 0 0 57 0 AZ
## 77 14 3 0 57 0 AZ
## 78 14 6 0 57 0 AZ
## 79 14 9 0 57 1 AZ
## 80 14 12 0 57 0 AZ
## 81 14 15 0 57 1 AZ
## 82 14 18 0 57 1 AZ
## 83 15 3 0 43 1 AZ
## 84 15 6 0 43 1 AZ
## 85 15 9 0 43 0 AZ
## 86 15 12 0 43 1 AZ
## 87 15 15 0 43 1 AZ
## 88 16 0 0 44 0 AZ
## 89 16 3 0 44 0 AZ
## 90 16 6 0 44 1 AZ
## 91 16 9 0 44 1 AZ
## 92 16 12 0 44 0 AZ
## 93 16 15 0 44 1 AZ
## 94 17 3 0 45 1 AZ
## 95 17 6 0 45 1 AZ
## 96 17 9 0 45 1 AZ
## 97 17 12 0 45 0 AZ
## 98 17 15 0 45 1 AZ
## 99 17 18 0 45 1 AZ
## 100 18 0 0 52 0 AZ
## 101 18 3 0 52 0 AZ
## 102 18 6 0 52 0 AZ
## 103 18 9 0 52 1 AZ
## 104 18 12 0 52 1 AZ
## 105 18 15 0 52 1 AZ
## 106 18 18 0 52 1 AZ
## 107 19 0 0 48 0 AZ
## 108 19 3 0 48 1 AZ
## 109 19 9 0 48 1 AZ
## 110 19 12 0 48 1 AZ
## 111 19 18 0 48 1 AZ
## 112 20 3 0 47 0 AZ
## 113 20 9 0 47 1 AZ
## 114 20 15 0 47 0 AZ
## 115 20 18 0 47 1 AZ
## 116 21 0 0 65 0 AZ
## 117 21 3 0 65 0 AZ
## 118 21 6 0 65 0 AZ
## 119 21 9 0 65 0 AZ
## 120 21 15 0 65 0 AZ
## 121 21 18 0 65 0 AZ
## 122 22 3 0 40 0 AZ
## 123 22 6 0 40 1 AZ
## 124 22 9 0 40 1 AZ
## 125 22 15 0 40 1 AZ
## 126 22 18 0 40 1 AZ
## 127 23 0 0 54 0 AZ
## 128 23 3 0 54 1 AZ
## 129 23 9 0 54 1 AZ
## 130 23 12 0 54 1 AZ
## 131 23 15 0 54 0 AZ
## 132 23 18 0 54 1 AZ
## 133 24 0 0 46 0 AZ
## 134 24 3 0 46 0 AZ
## 135 24 6 0 46 1 AZ
## 136 24 9 0 46 1 AZ
## 137 24 12 0 46 1 AZ
## 138 24 15 0 46 1 AZ
## 139 24 18 0 46 1 AZ
## 140 25 0 0 49 0 AZ
## 141 25 3 0 49 1 AZ
## 142 25 6 0 49 0 AZ
## 143 25 9 0 49 1 AZ
## 144 25 12 0 49 1 AZ
## 145 25 15 0 49 1 AZ
## 146 25 18 0 49 1 AZ
## 147 26 3 0 50 0 AZ
## 148 26 9 0 50 0 AZ
## 149 26 12 0 50 0 AZ
## 150 26 15 0 50 1 AZ
## 151 26 18 0 50 1 AZ
## 152 27 0 1 61 0 AZ
## 153 27 3 1 61 0 AZ
## 154 27 6 1 61 0 AZ
## 155 27 9 1 61 0 AZ
## 156 27 12 1 61 0 AZ
## 157 27 18 1 61 0 AZ
## 158 28 0 1 47 0 AZ
## 159 28 6 1 47 0 AZ
## 160 28 9 1 47 0 AZ
## 161 28 12 1 47 0 AZ
## 162 28 15 1 47 1 AZ
## 163 28 18 1 47 0 AZ
## 164 29 0 1 59 0 AZ
## 165 29 6 1 59 0 AZ
## 166 29 9 1 59 0 AZ
## 167 29 12 1 59 0 AZ
## 168 29 15 1 59 0 AZ
## 169 29 18 1 59 0 AZ
## 170 30 0 1 58 0 AZ
## 171 30 3 1 58 0 AZ
## 172 30 6 1 58 0 AZ
## 173 30 9 1 58 0 AZ
## 174 30 12 1 58 0 AZ
## 175 30 15 1 58 0 AZ
## 176 30 18 1 58 0 AZ
## 177 31 3 1 52 0 AZ
## 178 31 6 1 52 0 AZ
## 179 31 9 1 52 0 AZ
## 180 31 12 1 52 1 AZ
## 181 32 0 1 37 1 AZ
## 182 32 6 1 37 0 AZ
## 183 32 9 1 37 1 AZ
## 184 32 15 1 37 1 AZ
## 185 33 0 1 59 0 AZ
## 186 33 6 1 59 0 AZ
## 187 33 9 1 59 0 AZ
## 188 33 12 1 59 0 AZ
## 189 33 15 1 59 1 AZ
## 190 34 0 1 50 0 AZ
## 191 34 3 1 50 0 AZ
## 192 34 6 1 50 0 AZ
## 193 34 15 1 50 0 AZ
## 194 34 18 1 50 0 AZ
## 195 35 0 1 45 1 AZ
## 196 35 3 1 45 0 AZ
## 197 35 6 1 45 1 AZ
## 198 35 9 1 45 1 AZ
## 199 35 15 1 45 1 AZ
## 200 35 18 1 45 1 AZ
## 201 36 0 1 45 0 AZ
## 202 36 3 1 45 1 AZ
## 203 36 6 1 45 0 AZ
## 204 36 9 1 45 0 AZ
## 205 36 12 1 45 1 AZ
## 206 36 15 1 45 0 AZ
## 207 36 18 1 45 1 AZ
## 208 37 0 1 57 0 AZ
## 209 37 3 1 57 0 AZ
## 210 37 6 1 57 0 AZ
## 211 37 12 1 57 0 AZ
## 212 37 15 1 57 1 AZ
## 213 37 18 1 57 0 AZ
## 214 38 0 1 35 0 AZ
## 215 38 3 1 35 0 AZ
## 216 38 9 1 35 1 AZ
## 217 38 12 1 35 0 AZ
## 218 38 15 1 35 0 AZ
## 219 39 3 1 62 0 AZ
## 220 39 6 1 62 0 AZ
## 221 39 9 1 62 0 AZ
## 222 39 12 1 62 0 AZ
## 223 39 18 1 62 0 AZ
## 224 40 0 1 50 0 AZ
## 225 40 3 1 50 0 AZ
## 226 40 6 1 50 0 AZ
## 227 40 9 1 50 0 AZ
## 228 40 12 1 50 1 AZ
## 229 40 15 1 50 0 AZ
## 230 40 18 1 50 1 AZ
## 231 41 3 1 54 0 AZ
## 232 41 6 1 54 0 AZ
## 233 41 9 1 54 0 AZ
## 234 41 12 1 54 0 AZ
## 235 41 15 1 54 0 AZ
## 236 42 0 1 58 0 AZ
## 237 42 3 1 58 0 AZ
## 238 42 6 1 58 0 AZ
## 239 42 12 1 58 0 AZ
## 240 42 15 1 58 1 AZ
## 241 42 18 1 58 1 AZ
## 242 43 0 1 53 1 AZ
## 243 43 3 1 53 0 AZ
## 244 43 6 1 53 0 AZ
## 245 43 9 1 53 0 AZ
## 246 43 12 1 53 1 AZ
## 247 43 15 1 53 0 AZ
## 248 43 18 1 53 0 AZ
## 249 44 0 1 50 0 AZ
## 250 44 3 1 50 0 AZ
## 251 44 12 1 50 0 AZ
## 252 44 15 1 50 0 AZ
## 253 44 18 1 50 1 AZ
## 254 45 3 1 54 1 AZ
## 255 45 9 1 54 0 AZ
## 256 45 12 1 54 0 AZ
## 257 45 15 1 54 1 AZ
## 258 45 18 1 54 1 AZ
## 259 46 0 1 53 0 AZ
## 260 46 3 1 53 0 AZ
## 261 46 6 1 53 0 AZ
## 262 46 12 1 53 1 AZ
## 263 46 15 1 53 0 AZ
## 264 46 18 1 53 0 AZ
## 265 47 3 1 73 0 AZ
## 266 47 6 1 73 0 AZ
## 267 47 9 1 73 0 AZ
## 268 47 12 1 73 0 AZ
## 269 47 15 1 73 0 AZ
## 270 47 18 1 73 0 AZ
## 271 48 0 1 41 0 AZ
## 272 48 3 1 41 1 AZ
## 273 48 6 1 41 1 AZ
## 274 48 9 1 41 1 AZ
## 275 48 12 1 41 0 AZ
## 276 49 0 1 58 0 AZ
## 277 49 3 1 58 0 AZ
## 278 49 6 1 58 0 AZ
## 279 49 9 1 58 0 AZ
## 280 49 15 1 58 0 AZ
## 281 49 18 1 58 0 AZ
## 282 50 0 1 41 0 AZ
## 283 50 3 1 41 0 AZ
## 284 50 6 1 41 1 AZ
## 285 50 9 1 41 1 AZ
## 286 50 12 1 41 1 AZ
## 287 50 15 1 41 1 AZ
## 288 50 18 1 41 1 AZ
## 289 51 0 1 63 0 AZ
## 290 51 6 1 63 0 AZ
## 291 51 9 1 63 0 AZ
## 292 51 12 1 63 0 AZ
## 293 51 15 1 63 0 AZ
## 294 51 18 1 63 0 AZ
## 295 52 3 1 47 1 AZ
## 296 52 6 1 47 1 AZ
## 297 52 9 1 47 1 AZ
## 298 52 12 1 47 1 AZ
## 299 52 15 1 47 1 AZ
## 300 52 18 1 47 1 AZ
## 301 53 0 1 59 0 AZ
## 302 53 3 1 59 0 AZ
## 303 53 6 1 59 0 AZ
## 304 53 9 1 59 0 AZ
## 305 53 12 1 59 0 AZ
## 306 53 15 1 59 1 AZ
## 307 54 0 1 54 0 AZ
## 308 54 3 1 54 0 AZ
## 309 54 12 1 54 0 AZ
## 310 54 18 1 54 0 AZ
## 311 55 0 1 53 0 AZ
## 312 55 3 1 53 1 AZ
## 313 55 9 1 53 1 AZ
## 314 55 12 1 53 1 AZ
## 315 55 15 1 53 1 AZ
## 316 55 18 1 53 1 AZ
## 317 56 0 1 45 0 AZ
## 318 56 9 1 45 0 AZ
## 319 56 15 1 45 0 AZ
## 320 57 0 1 55 0 AZ
## 321 57 3 1 55 1 AZ
## 322 57 6 1 55 0 AZ
## 323 57 9 1 55 0 AZ
## 324 57 12 1 55 0 AZ
## 325 57 15 1 55 1 AZ
## 326 58 0 1 58 0 AZ
## 327 58 3 1 58 0 AZ
## 328 58 6 1 58 0 AZ
## 329 58 9 1 58 0 AZ
## 330 58 12 1 58 1 AZ
## 331 58 15 1 58 0 AZ
## 332 58 18 1 58 1 AZ
## 333 59 0 1 55 0 AZ
## 334 59 6 1 55 1 AZ
## 335 59 9 1 55 0 AZ
## 336 59 12 1 55 0 AZ
## 337 59 15 1 55 0 AZ
## 338 59 18 1 55 1 AZ
## 339 60 0 1 52 0 AZ
## 340 60 3 1 52 0 AZ
## 341 60 6 1 52 0 AZ
## 342 60 9 1 52 0 AZ
## 343 60 12 1 52 0 AZ
## 344 60 15 1 52 0 AZ
## 345 60 18 1 52 0 AZ
## 346 61 3 1 44 0 AZ
## 347 61 9 1 44 1 AZ
## 348 61 15 1 44 1 AZ
## 349 62 0 1 45 0 AZ
## 350 62 6 1 45 0 AZ
## 351 62 15 1 45 0 AZ
## 352 63 0 1 47 0 AZ
## 353 63 3 1 47 0 AZ
## 354 63 6 1 47 0 AZ
## 355 63 9 1 47 0 AZ
## 356 63 15 1 47 1 AZ
## 357 63 18 1 47 0 AZ
## 358 64 0 1 53 0 AZ
## 359 64 3 1 53 0 AZ
## 360 64 6 1 53 0 AZ
## 361 64 9 1 53 0 AZ
## 362 64 12 1 53 1 AZ
## 363 64 15 1 53 1 AZ
## 364 64 18 1 53 1 AZ
## 365 65 0 1 47 1 AZ
## 366 65 6 1 47 0 AZ
## 367 65 9 1 47 1 AZ
## 368 65 12 1 47 1 AZ
## 369 65 15 1 47 1 AZ
## 370 65 18 1 47 1 AZ
## 371 66 0 1 44 1 AZ
## 372 66 3 1 44 1 AZ
## 373 66 6 1 44 1 AZ
## 374 66 9 1 44 1 AZ
## 375 66 12 1 44 1 AZ
## 376 66 15 1 44 1 AZ
## 377 66 18 1 44 1 AZ
## 378 67 0 1 42 1 AZ
## 379 67 3 1 42 1 AZ
## 380 67 9 1 42 1 AZ
## 381 67 12 1 42 0 AZ
## 382 67 15 1 42 1 AZ
## 383 67 18 1 42 0 AZ
## 384 68 0 1 52 0 AZ
## 385 68 3 1 52 0 AZ
## 386 68 6 1 52 0 AZ
## 387 68 9 1 52 0 AZ
## 388 68 12 1 52 1 AZ
## 389 68 15 1 52 1 AZ
## 390 69 3 1 48 0 AZ
## 391 69 6 1 48 0 AZ
## 392 69 9 1 48 1 AZ
## 393 69 12 1 48 0 AZ
## 394 69 15 1 48 0 AZ
## 395 70 0 1 46 0 AZ
## 396 70 3 1 46 0 AZ
## 397 70 6 1 46 1 AZ
## 398 70 9 1 46 0 AZ
## 399 70 12 1 46 1 AZ
## 400 70 15 1 46 1 AZ
## 401 71 0 1 49 1 AZ
## 402 71 6 1 49 1 AZ
## 403 71 9 1 49 1 AZ
## 404 71 12 1 49 1 AZ
## 405 71 15 1 49 0 AZ
## 406 71 18 1 49 1 AZ
## 407 72 6 1 51 0 AZ
## 408 72 9 1 51 0 AZ
## 409 72 15 1 51 0 AZ
## 410 72 18 1 51 1 AZ
## 411 73 0 1 53 1 AZ
## 412 73 3 1 53 0 AZ
## 413 73 6 1 53 0 AZ
## 414 73 9 1 53 0 AZ
## 415 73 12 1 53 0 AZ
## 416 73 18 1 53 1 AZ
## 417 74 0 1 45 0 AZ
## 418 74 3 1 45 0 AZ
## 419 74 6 1 45 0 AZ
## 420 74 9 1 45 0 AZ
## 421 74 12 1 45 1 AZ
## 422 74 18 1 45 1 AZ
## 423 75 0 1 57 0 AZ
## 424 75 3 1 57 0 AZ
## 425 75 6 1 57 0 AZ
## 426 75 9 1 57 0 AZ
## 427 75 12 1 57 0 AZ
## 428 75 15 1 57 0 AZ
## 429 75 18 1 57 1 AZ
## 430 76 0 0 53 0 AZ+MP
## 431 76 3 0 53 1 AZ+MP
## 432 76 9 0 53 1 AZ+MP
## 433 76 12 0 53 1 AZ+MP
## 434 76 18 0 53 1 AZ+MP
## 435 77 0 0 46 0 AZ+MP
## 436 77 3 0 46 0 AZ+MP
## 437 77 6 0 46 0 AZ+MP
## 438 77 9 0 46 0 AZ+MP
## 439 77 12 0 46 1 AZ+MP
## 440 77 15 0 46 1 AZ+MP
## 441 77 18 0 46 1 AZ+MP
## 442 78 0 0 55 1 AZ+MP
## 443 78 6 0 55 1 AZ+MP
## 444 78 9 0 55 1 AZ+MP
## 445 78 12 0 55 0 AZ+MP
## 446 78 15 0 55 1 AZ+MP
## 447 79 0 0 45 0 AZ+MP
## 448 79 3 0 45 1 AZ+MP
## 449 79 6 0 45 1 AZ+MP
## 450 79 9 0 45 0 AZ+MP
## 451 79 12 0 45 1 AZ+MP
## 452 79 18 0 45 1 AZ+MP
## 453 80 0 0 61 0 AZ+MP
## 454 80 3 0 61 0 AZ+MP
## 455 80 6 0 61 1 AZ+MP
## 456 80 12 0 61 0 AZ+MP
## 457 80 15 0 61 0 AZ+MP
## 458 81 0 0 43 0 AZ+MP
## 459 81 3 0 43 1 AZ+MP
## 460 81 6 0 43 0 AZ+MP
## 461 81 15 0 43 1 AZ+MP
## 462 82 0 0 44 1 AZ+MP
## 463 82 6 0 44 1 AZ+MP
## 464 82 9 0 44 1 AZ+MP
## 465 82 12 0 44 1 AZ+MP
## 466 82 15 0 44 1 AZ+MP
## 467 82 18 0 44 1 AZ+MP
## 468 83 3 0 45 0 AZ+MP
## 469 83 6 0 45 1 AZ+MP
## 470 83 9 0 45 0 AZ+MP
## 471 83 12 0 45 1 AZ+MP
## 472 84 0 0 55 0 AZ+MP
## 473 84 6 0 55 0 AZ+MP
## 474 84 12 0 55 0 AZ+MP
## 475 84 18 0 55 0 AZ+MP
## 476 85 0 0 67 0 AZ+MP
## 477 85 3 0 67 0 AZ+MP
## 478 85 12 0 67 0 AZ+MP
## 479 85 15 0 67 0 AZ+MP
## 480 85 18 0 67 0 AZ+MP
## 481 86 0 0 52 0 AZ+MP
## 482 86 12 0 52 0 AZ+MP
## 483 86 18 0 52 1 AZ+MP
## 484 87 0 0 50 1 AZ+MP
## 485 87 9 0 50 1 AZ+MP
## 486 87 12 0 50 1 AZ+MP
## 487 87 18 0 50 1 AZ+MP
## 488 88 0 0 54 0 AZ+MP
## 489 88 9 0 54 0 AZ+MP
## 490 88 12 0 54 1 AZ+MP
## 491 88 15 0 54 1 AZ+MP
## 492 89 3 0 51 0 AZ+MP
## 493 89 6 0 51 0 AZ+MP
## 494 89 12 0 51 1 AZ+MP
## 495 89 18 0 51 1 AZ+MP
## 496 90 3 0 61 0 AZ+MP
## 497 90 6 0 61 1 AZ+MP
## 498 90 9 0 61 0 AZ+MP
## 499 90 12 0 61 0 AZ+MP
## 500 90 18 0 61 1 AZ+MP
## 501 91 0 0 52 0 AZ+MP
## 502 91 3 0 52 0 AZ+MP
## 503 91 6 0 52 1 AZ+MP
## 504 91 9 0 52 1 AZ+MP
## 505 91 12 0 52 1 AZ+MP
## 506 91 15 0 52 1 AZ+MP
## 507 92 0 0 50 0 AZ+MP
## 508 92 3 0 50 0 AZ+MP
## 509 92 6 0 50 1 AZ+MP
## 510 92 9 0 50 1 AZ+MP
## 511 92 12 0 50 1 AZ+MP
## 512 92 15 0 50 1 AZ+MP
## 513 93 0 0 45 0 AZ+MP
## 514 93 3 0 45 0 AZ+MP
## 515 93 6 0 45 0 AZ+MP
## 516 93 9 0 45 1 AZ+MP
## 517 93 12 0 45 1 AZ+MP
## 518 93 15 0 45 1 AZ+MP
## 519 94 0 0 58 0 AZ+MP
## 520 94 3 0 58 0 AZ+MP
## 521 94 6 0 58 0 AZ+MP
## 522 94 9 0 58 0 AZ+MP
## 523 94 12 0 58 0 AZ+MP
## 524 94 18 0 58 0 AZ+MP
## 525 95 0 0 51 0 AZ+MP
## 526 95 3 0 51 1 AZ+MP
## 527 95 6 0 51 0 AZ+MP
## 528 95 9 0 51 1 AZ+MP
## 529 95 12 0 51 1 AZ+MP
## 530 95 15 0 51 1 AZ+MP
## 531 95 18 0 51 1 AZ+MP
## 532 96 0 0 59 0 AZ+MP
## 533 96 6 0 59 0 AZ+MP
## 534 96 9 0 59 0 AZ+MP
## 535 96 12 0 59 1 AZ+MP
## 536 96 15 0 59 1 AZ+MP
## 537 97 0 0 49 0 AZ+MP
## 538 97 3 0 49 0 AZ+MP
## 539 97 6 0 49 0 AZ+MP
## 540 97 9 0 49 1 AZ+MP
## 541 97 15 0 49 1 AZ+MP
## 542 97 18 0 49 1 AZ+MP
## 543 98 0 0 55 0 AZ+MP
## 544 98 6 0 55 1 AZ+MP
## 545 98 9 0 55 0 AZ+MP
## 546 98 12 0 55 0 AZ+MP
## 547 98 18 0 55 0 AZ+MP
## 548 99 0 0 32 1 AZ+MP
## 549 99 3 0 32 1 AZ+MP
## 550 99 6 0 32 1 AZ+MP
## 551 99 12 0 32 1 AZ+MP
## 552 99 15 0 32 1 AZ+MP
## 553 99 18 0 32 1 AZ+MP
## 554 100 0 0 45 0 AZ+MP
## 555 100 3 0 45 1 AZ+MP
## 556 100 6 0 45 1 AZ+MP
## 557 100 9 0 45 1 AZ+MP
## 558 100 12 0 45 1 AZ+MP
## 559 100 18 0 45 1 AZ+MP
## 560 101 0 0 47 0 AZ+MP
## 561 101 3 0 47 1 AZ+MP
## 562 101 6 0 47 1 AZ+MP
## 563 101 9 0 47 1 AZ+MP
## 564 101 15 0 47 1 AZ+MP
## 565 101 18 0 47 1 AZ+MP
## 566 102 0 0 49 0 AZ+MP
## 567 102 3 0 49 1 AZ+MP
## 568 102 6 0 49 0 AZ+MP
## 569 102 9 0 49 1 AZ+MP
## 570 102 15 0 49 1 AZ+MP
## 571 102 18 0 49 1 AZ+MP
## 572 103 0 0 47 1 AZ+MP
## 573 103 3 0 47 0 AZ+MP
## 574 103 6 0 47 1 AZ+MP
## 575 103 12 0 47 1 AZ+MP
## 576 103 15 0 47 1 AZ+MP
## 577 103 18 0 47 1 AZ+MP
## 578 104 0 0 49 1 AZ+MP
## 579 104 3 0 49 0 AZ+MP
## 580 104 9 0 49 1 AZ+MP
## 581 104 12 0 49 1 AZ+MP
## 582 104 15 0 49 1 AZ+MP
## 583 104 18 0 49 1 AZ+MP
## 584 105 0 0 48 0 AZ+MP
## 585 105 3 0 48 1 AZ+MP
## 586 105 9 0 48 0 AZ+MP
## 587 105 12 0 48 0 AZ+MP
## 588 105 15 0 48 0 AZ+MP
## 589 105 18 0 48 1 AZ+MP
## 590 106 0 0 51 0 AZ+MP
## 591 106 3 0 51 1 AZ+MP
## 592 106 9 0 51 1 AZ+MP
## 593 106 15 0 51 1 AZ+MP
## 594 106 18 0 51 1 AZ+MP
## 595 107 3 0 49 1 AZ+MP
## 596 107 6 0 49 0 AZ+MP
## 597 107 9 0 49 0 AZ+MP
## 598 107 12 0 49 1 AZ+MP
## 599 107 18 0 49 1 AZ+MP
## 600 108 3 1 46 1 AZ+MP
## 601 108 6 1 46 1 AZ+MP
## 602 108 9 1 46 0 AZ+MP
## 603 108 12 1 46 1 AZ+MP
## 604 108 15 1 46 1 AZ+MP
## 605 108 18 1 46 1 AZ+MP
## 606 109 0 1 64 0 AZ+MP
## 607 109 6 1 64 0 AZ+MP
## 608 109 9 1 64 0 AZ+MP
## 609 109 15 1 64 0 AZ+MP
## 610 109 18 1 64 0 AZ+MP
## 611 110 0 1 37 0 AZ+MP
## 612 110 6 1 37 1 AZ+MP
## 613 110 9 1 37 1 AZ+MP
## 614 110 12 1 37 1 AZ+MP
## 615 110 15 1 37 1 AZ+MP
## 616 110 18 1 37 1 AZ+MP
## 617 111 0 1 52 0 AZ+MP
## 618 111 3 1 52 0 AZ+MP
## 619 111 6 1 52 0 AZ+MP
## 620 111 9 1 52 0 AZ+MP
## 621 111 12 1 52 1 AZ+MP
## 622 111 15 1 52 0 AZ+MP
## 623 111 18 1 52 0 AZ+MP
## 624 112 0 1 46 1 AZ+MP
## 625 112 3 1 46 0 AZ+MP
## 626 112 6 1 46 0 AZ+MP
## 627 112 9 1 46 0 AZ+MP
## 628 112 12 1 46 0 AZ+MP
## 629 112 15 1 46 1 AZ+MP
## 630 112 18 1 46 1 AZ+MP
## 631 113 0 1 40 1 AZ+MP
## 632 113 3 1 40 1 AZ+MP
## 633 113 6 1 40 1 AZ+MP
## 634 113 9 1 40 1 AZ+MP
## 635 113 12 1 40 1 AZ+MP
## 636 113 15 1 40 1 AZ+MP
## 637 114 0 1 43 0 AZ+MP
## 638 114 3 1 43 1 AZ+MP
## 639 114 9 1 43 1 AZ+MP
## 640 114 12 1 43 1 AZ+MP
## 641 114 15 1 43 1 AZ+MP
## 642 114 18 1 43 1 AZ+MP
## 643 115 6 1 52 0 AZ+MP
## 644 115 9 1 52 0 AZ+MP
## 645 115 12 1 52 1 AZ+MP
## 646 115 15 1 52 1 AZ+MP
## 647 115 18 1 52 1 AZ+MP
## 648 116 0 1 54 0 AZ+MP
## 649 116 3 1 54 0 AZ+MP
## 650 116 6 1 54 0 AZ+MP
## 651 116 15 1 54 1 AZ+MP
## 652 116 18 1 54 1 AZ+MP
## 653 117 0 1 50 0 AZ+MP
## 654 117 3 1 50 0 AZ+MP
## 655 117 6 1 50 1 AZ+MP
## 656 117 9 1 50 0 AZ+MP
## 657 117 12 1 50 1 AZ+MP
## 658 117 15 1 50 0 AZ+MP
## 659 117 18 1 50 1 AZ+MP
## 660 118 0 1 53 0 AZ+MP
## 661 118 6 1 53 0 AZ+MP
## 662 118 9 1 53 0 AZ+MP
## 663 118 12 1 53 0 AZ+MP
## 664 118 15 1 53 1 AZ+MP
## 665 118 18 1 53 1 AZ+MP
## 666 119 0 1 57 0 AZ+MP
## 667 119 3 1 57 0 AZ+MP
## 668 119 6 1 57 0 AZ+MP
## 669 119 9 1 57 0 AZ+MP
## 670 119 15 1 57 1 AZ+MP
## 671 119 18 1 57 1 AZ+MP
## 672 120 0 1 41 0 AZ+MP
## 673 120 6 1 41 1 AZ+MP
## 674 120 12 1 41 1 AZ+MP
## 675 120 15 1 41 1 AZ+MP
## 676 120 18 1 41 1 AZ+MP
## 677 121 3 1 48 0 AZ+MP
## 678 121 6 1 48 0 AZ+MP
## 679 121 12 1 48 1 AZ+MP
## 680 121 15 1 48 1 AZ+MP
## 681 121 18 1 48 1 AZ+MP
## 682 122 0 1 61 0 AZ+MP
## 683 122 3 1 61 0 AZ+MP
## 684 122 6 1 61 0 AZ+MP
## 685 122 9 1 61 1 AZ+MP
## 686 122 12 1 61 0 AZ+MP
## 687 122 18 1 61 0 AZ+MP
## 688 123 3 1 39 0 AZ+MP
## 689 123 6 1 39 1 AZ+MP
## 690 123 9 1 39 0 AZ+MP
## 691 123 12 1 39 1 AZ+MP
## 692 123 15 1 39 1 AZ+MP
## 693 123 18 1 39 1 AZ+MP
## 694 124 0 1 42 0 AZ+MP
## 695 124 3 1 42 0 AZ+MP
## 696 124 6 1 42 0 AZ+MP
## 697 124 9 1 42 0 AZ+MP
## 698 124 12 1 42 1 AZ+MP
## 699 124 15 1 42 1 AZ+MP
## 700 124 18 1 42 1 AZ+MP
## 701 125 0 1 57 0 AZ+MP
## 702 125 3 1 57 0 AZ+MP
## 703 125 6 1 57 0 AZ+MP
## 704 125 12 1 57 0 AZ+MP
## 705 125 15 1 57 1 AZ+MP
## 706 126 3 1 52 1 AZ+MP
## 707 126 15 1 52 1 AZ+MP
## 708 127 0 1 45 0 AZ+MP
## 709 127 3 1 45 0 AZ+MP
## 710 127 6 1 45 0 AZ+MP
## 711 127 9 1 45 0 AZ+MP
## 712 127 12 1 45 0 AZ+MP
## 713 127 15 1 45 1 AZ+MP
## 714 127 18 1 45 0 AZ+MP
## 715 128 0 1 54 0 AZ+MP
## 716 128 3 1 54 1 AZ+MP
## 717 128 6 1 54 0 AZ+MP
## 718 128 12 1 54 1 AZ+MP
## 719 128 15 1 54 1 AZ+MP
## 720 129 0 1 49 0 AZ+MP
## 721 129 3 1 49 0 AZ+MP
## 722 129 6 1 49 0 AZ+MP
## 723 129 9 1 49 1 AZ+MP
## 724 129 12 1 49 1 AZ+MP
## 725 129 18 1 49 1 AZ+MP
## 726 130 0 1 60 0 AZ+MP
## 727 130 3 1 60 0 AZ+MP
## 728 130 6 1 60 0 AZ+MP
## 729 130 9 1 60 0 AZ+MP
## 730 130 12 1 60 0 AZ+MP
## 731 130 15 1 60 0 AZ+MP
## 732 131 0 1 49 0 AZ+MP
## 733 131 3 1 49 0 AZ+MP
## 734 131 6 1 49 0 AZ+MP
## 735 131 9 1 49 1 AZ+MP
## 736 131 12 1 49 1 AZ+MP
## 737 131 15 1 49 0 AZ+MP
## 738 131 18 1 49 1 AZ+MP
## 739 132 0 1 48 0 AZ+MP
## 740 132 3 1 48 0 AZ+MP
## 741 132 6 1 48 0 AZ+MP
## 742 132 9 1 48 0 AZ+MP
## 743 132 12 1 48 1 AZ+MP
## 744 132 15 1 48 1 AZ+MP
## 745 132 18 1 48 1 AZ+MP
## 746 133 0 1 54 0 AZ+MP
## 747 133 3 1 54 0 AZ+MP
## 748 133 6 1 54 1 AZ+MP
## 749 133 12 1 54 1 AZ+MP
## 750 133 18 1 54 1 AZ+MP
## 751 134 0 1 54 0 AZ+MP
## 752 134 3 1 54 0 AZ+MP
## 753 134 9 1 54 1 AZ+MP
## 754 134 15 1 54 1 AZ+MP
## 755 134 18 1 54 1 AZ+MP
## 756 135 0 1 52 0 AZ+MP
## 757 135 3 1 52 0 AZ+MP
## 758 135 6 1 52 1 AZ+MP
## 759 135 9 1 52 0 AZ+MP
## 760 135 12 1 52 1 AZ+MP
## 761 135 15 1 52 1 AZ+MP
## 762 135 18 1 52 1 AZ+MP
## 763 136 0 1 55 0 AZ+MP
## 764 136 6 1 55 0 AZ+MP
## 765 136 9 1 55 0 AZ+MP
## 766 136 12 1 55 1 AZ+MP
## 767 136 18 1 55 1 AZ+MP
## 768 137 0 1 48 0 AZ+MP
## 769 137 3 1 48 1 AZ+MP
## 770 137 6 1 48 1 AZ+MP
## 771 137 9 1 48 1 AZ+MP
## 772 137 15 1 48 1 AZ+MP
## 773 137 18 1 48 1 AZ+MP
## 774 138 0 1 54 0 AZ+MP
## 775 138 3 1 54 0 AZ+MP
## 776 138 6 1 54 0 AZ+MP
## 777 138 9 1 54 0 AZ+MP
## 778 138 12 1 54 1 AZ+MP
## 779 138 15 1 54 0 AZ+MP
## 780 138 18 1 54 0 AZ+MP
## 781 139 0 1 54 0 AZ+MP
## 782 139 3 1 54 0 AZ+MP
## 783 139 6 1 54 0 AZ+MP
## 784 139 12 1 54 0 AZ+MP
## 785 139 18 1 54 0 AZ+MP
## 786 140 3 1 54 0 AZ+MP
## 787 140 6 1 54 0 AZ+MP
## 788 140 12 1 54 1 AZ+MP
## 789 140 15 1 54 1 AZ+MP
## 790 140 18 1 54 0 AZ+MP
## 791 141 0 1 52 0 AZ+MP
## 792 141 3 1 52 0 AZ+MP
## 793 141 9 1 52 1 AZ+MP
## 794 141 12 1 52 0 AZ+MP
## 795 141 15 1 52 1 AZ+MP
## 796 141 18 1 52 1 AZ+MP
## 797 142 6 1 52 0 AZ+MP
## 798 142 9 1 52 0 AZ+MP
## 799 142 12 1 52 0 AZ+MP
## 800 142 18 1 52 1 AZ+MP
## 801 143 0 1 56 0 AZ+MP
## 802 143 3 1 56 0 AZ+MP
## 803 143 6 1 56 0 AZ+MP
## 804 143 9 1 56 0 AZ+MP
## 805 143 12 1 56 0 AZ+MP
## 806 143 18 1 56 0 AZ+MP
## 807 144 0 1 38 1 AZ+MP
## 808 144 3 1 38 1 AZ+MP
## 809 144 6 1 38 1 AZ+MP
## 810 144 9 1 38 1 AZ+MP
## 811 144 12 1 38 1 AZ+MP
## 812 144 18 1 38 1 AZ+MP
## 813 145 0 1 57 0 AZ+MP
## 814 145 3 1 57 0 AZ+MP
## 815 145 6 1 57 0 AZ+MP
## 816 145 9 1 57 0 AZ+MP
## 817 145 12 1 57 0 AZ+MP
## 818 145 15 1 57 1 AZ+MP
## 819 145 18 1 57 1 AZ+MP
## 820 146 0 1 43 1 AZ+MP
## 821 146 3 1 43 0 AZ+MP
## 822 146 6 1 43 1 AZ+MP
## 823 146 9 1 43 1 AZ+MP
## 824 146 15 1 43 1 AZ+MP
## 825 146 18 1 43 1 AZ+MP
## 826 147 0 1 46 0 AZ+MP
## 827 147 6 1 46 0 AZ+MP
## 828 147 9 1 46 0 AZ+MP
## 829 147 12 1 46 0 AZ+MP
## 830 147 18 1 46 1 AZ+MP
## 831 148 0 1 51 0 AZ+MP
## 832 148 3 1 51 0 AZ+MP
## 833 148 9 1 51 0 AZ+MP
## 834 148 12 1 51 0 AZ+MP
## 835 148 15 1 51 1 AZ+MP
## 836 149 0 1 55 0 AZ+MP
## 837 149 6 1 55 0 AZ+MP
## 838 149 9 1 55 1 AZ+MP
## 839 149 12 1 55 0 AZ+MP
## 840 149 15 1 55 1 AZ+MP
## 841 149 18 1 55 1 AZ+MP
## 842 150 0 1 41 1 AZ+MP
## 843 150 3 1 41 1 AZ+MP
## 844 150 6 1 41 1 AZ+MP
## 845 150 9 1 41 1 AZ+MP
## 846 150 12 1 41 1 AZ+MP
## 847 150 18 1 41 1 AZ+MP
Esto implica que el modelo en la escala del PL tiene interceptos diferentes para cada paciente!
prop.table
frec <- table(escle$mejora, escle$grupo)
frec
##
## AZ AZ+MP
## 0 251 204
## 1 178 214
round(prop.table(frec, margin = 2), 2) # porcentajes de columnas ~
##
## AZ AZ+MP
## 0 0.59 0.49
## 1 0.41 0.51
summary(escle)
## id tiempo tratpre edad
## Min. : 1.0 Min. : 0.000 Min. :0.000 Min. :32.00
## 1st Qu.: 37.0 1st Qu.: 3.000 1st Qu.:0.000 1st Qu.:46.00
## Median : 75.0 Median : 9.000 Median :1.000 Median :50.00
## Mean : 75.4 Mean : 8.968 Mean :0.621 Mean :50.41
## 3rd Qu.:113.0 3rd Qu.:15.000 3rd Qu.:1.000 3rd Qu.:54.00
## Max. :150.0 Max. :18.000 Max. :1.000 Max. :73.00
## mejora grupo
## Min. :0.0000 AZ :429
## 1st Qu.:0.0000 AZ+MP:418
## Median :0.0000
## Mean :0.4628
## 3rd Qu.:1.0000
## Max. :1.0000
doBy
Sintaxis: usa una fórmula para elegir variables en función de los niveles de un factor:
doBy::summaryBy(
edad + tratpre + mejora ~ grupo, data = escle,
FUN = function(x) { round(c(m = mean(x), s = sd(x)), 2) }
)
## grupo edad.m edad.s tratpre.m tratpre.s mejora.m mejora.s
## 1 AZ 50.77 6.93 0.65 0.48 0.41 0.49
## 2 AZ+MP 50.04 6.43 0.59 0.49 0.51 0.50
escle %>%
group_by(grupo, tiempo) %>%
summarise(
count = n(),
perc = 100 * n()/nrow(escle)
) %>%
ggplot() +
aes(tiempo, count, fill = grupo) +
geom_col(position = position_dodge())
Como hay medidas repetidas por paciente y una estructura de correlación para dar cuenta allí, ya que las medidas del mismo paciente van a estar correlacionadas, usamos un efecto aleatorio:
escle_m0 <- glmer(mejora ~ grupo * tiempo + edad + tratpre + (1|id), data = escle, family = binomial)
summary(escle_m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: mejora ~ grupo * tiempo + edad + tratpre + (1 | id)
## Data: escle
##
## AIC BIC logLik deviance df.resid
## 854.6 887.8 -420.3 840.6 840
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7789 -0.5382 -0.1591 0.5539 4.0618
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.8119 0.901
## Number of obs: 847, groups: id, 150
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.51079 1.06166 7.075 1.50e-12 ***
## grupoAZ+MP -0.05404 0.38690 -0.140 0.888909
## tiempo 0.17243 0.02358 7.313 2.62e-13 ***
## edad -0.17934 0.02145 -8.360 < 2e-16 ***
## tratpre -0.92071 0.24395 -3.774 0.000161 ***
## grupoAZ+MP:tiempo 0.07107 0.03377 2.105 0.035325 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grAZ+MP tiempo edad tratpr
## grupoAZ+MP -0.196
## tiempo 0.026 0.533
## edad -0.957 0.008 -0.233
## tratpre -0.175 0.048 -0.085 0.052
## grpAZ+MP:tm 0.218 -0.787 -0.599 -0.085 -0.048
escle$edad_c <- scale(escle$edad, center = T, scale = F)
nAGQ
El parámetro nAGQ es el número de puntos de integración para aproximar la integral de la verosimilitud.
escle_m1b <- glmer(mejora ~ grupo * tiempo + edad_c + tratpre + (1|id), data = escle, family = binomial, nAGQ = 5)
glmmTMB también puede resolver este problema:
escle_m1c <- glmmTMB(mejora ~ grupo * tiempo + edad_c + tratpre + (1|id), data = escle, family = binomial)
Les profes usan drop1 suelto acá, sin especificar test:
drop1(escle_m1c)
## Single term deletions
##
## Model:
## mejora ~ grupo * tiempo + edad_c + tratpre + (1 | id)
## Df AIC
## <none> 854.60
## edad_c 1 930.54
## tratpre 1 866.63
## grupo:tiempo 1 857.13
patient_effects <- ranef(escle_m0)$id$`(Intercept)`
names(patient_effects) <- rownames(ranef(escle_m0)$id)
patient_effects
## 1 2 3 4 5
## 0.117113718 0.021200908 0.212700345 0.380614820 0.146157192
## 6 7 8 9 10
## 0.226520233 0.330004352 -0.018768412 -1.303321749 0.178915538
## 11 12 13 14 15
## -1.209356689 0.628728338 -0.554438248 0.416330997 0.072921247
## 16 17 18 19 20
## -0.459750135 0.244394800 0.330004352 0.541942949 -0.406158259
## 21 22 23 24 25
## -0.354212138 -0.156923235 0.666703870 0.159047194 0.432879471
## 26 27 28 29 30
## -0.463279471 -0.281312126 -0.784358399 -0.432750911 -0.495802347
## 31 32 33 34 35
## 0.092949093 0.019498678 0.312234565 -0.833395504 0.858642786
## 36 37 38 39 40
## -0.228654171 0.030082029 -1.421718431 -0.241585062 -0.153156772
## 41 42 43 44 45
## -0.534300908 0.596841227 0.111589606 -0.411595334 0.729174646
## 46 47 48 49 50
## -0.249786392 -0.059304788 0.086675102 -0.438188010 0.171105730
## 51 52 53 54 55
## -0.265634208 1.260904010 0.289329344 -0.508032932 1.444458247
## 56 57 58 59 60
## -0.805057300 0.532160576 0.524285971 0.316288471 -0.871291957
## 61 62 63 64 65
## 0.183346095 -0.746396770 -0.669031595 0.520350382 0.876858666
## 66 67 68 69 70
## 1.271372867 0.080690016 0.293255605 -0.440687480 0.185130824
## 71 72 73 74 75
## 1.047614315 -0.288958597 0.317619526 -0.350784111 -0.026076183
## 76 77 78 79 80
## 0.749392717 -0.889278166 0.832117072 -0.136106514 -0.059331745
## 81 82 83 84 85
## -0.406071615 0.625984444 -0.490570366 -1.018191438 -0.524973212
## 86 87 88 89 90
## -0.492031365 0.796352444 -0.031510210 -0.225100068 0.263538576
## 91 92 93 94 95
## 0.450226225 0.283259630 -0.567841831 -0.959858032 0.385631796
## 96 97 98 99 100
## 0.207624327 -0.300258201 -0.633166262 0.196045605 0.330805279
## 101 102 103 104 105
## 0.461188552 0.152377092 0.424801839 0.501877171 -0.991499089
## 106 107 108 109 110
## 0.576803213 -0.177378779 0.494543602 -0.403973663 -0.034658826
## 111 112 113 114 115
## -0.774979014 -0.441136978 0.841516082 0.441013579 0.165874007
## 116 117 118 119 120
## 0.170705199 -0.091254301 -0.221773018 0.239481189 0.124830887
## 121 122 123 124 125
## 0.000680692 0.074028948 -0.503427852 -0.790694728 0.010355975
## 126 127 128 129 130
## 0.769334228 -1.363415843 0.821403525 0.165260079 -0.451304934
## 131 132 133 134 135
## -0.178579041 -0.266034887 0.772903411 0.611964494 0.499708941
## 136 137 138 139 140
## 0.221690426 0.933263665 -0.614051959 -0.810992003 -0.038116503
## 141 142 143 144 145
## 0.183818530 -0.464915135 -0.765588068 0.709034963 0.075315992
## 146 147 148 149 150
## 0.546276974 -0.912824676 -0.442241556 0.383280473 0.902771800
escle$prediccion <- predict(escle_m0, type = "link")
escle$patient_effect <- sapply(escle$id, function(id) { patient_effects[id] })
ggplot(escle) +
facet_grid(. ~ grupo) +
aes(tiempo, prediccion + patient_effect) +
geom_line(aes(group = id), alpha = .5)
NO terminé esto en orden…
escle$prediccion_VR <- predict(escle_m0, type = "response")
ggplot(escle) +
aes(tiempo, prediccion_VR, color = grupo) +
geom_point(shape = 1, size = 2, alpha = .7) +
labs(
y = "Probabilidad de mejorar"
)
# geom_line(aes(group = grupo))
num_obs <- 100
grilla_tiempo <- seq(0, 18, length.out = num_obs)
AZ_MP <- rep(0:1, each = num_obs)
betas <- fixef(escle_m0)
PL.cond <-
betas["(Intercept)"] +
betas["grupoAZ+MP"] * AZ_MP +
betas["tiempo"] * grilla_tiempo +
betas["grupoAZ+MP:tiempo"] * AZ_MP * grilla_tiempo
datos_simu <- as.data.frame(cbind(
grilla_tiempo,
AZ_MP,
PL.cond
))
datos_simu$AZ_MP <- factor(datos_simu$AZ_MP)
ggplot(datos_simu) +
aes(grilla_tiempo, PL.cond, color = AZ_MP) +
geom_line(aes(group = AZ_MP))
PL <- predict(escle_m0, type = "link")
odds <- exp(PL)
prediccion_VR <- predict(escle_m0, type = "response")
sjPlot::plot_model(escle_m0, type = "eff")
## $grupo
##
## $tiempo
##
## $edad
##
## $tratpre
sjPlot::plot_model(escle_m0, type = "re")
sjPlot::plot_model(escle_m0, type = "est")
sjPlot::plot_model(escle_m0, type = "resid")
## Warning: Interaction terms are not supported by this plot type. Output for
## interaction terms may be inappropriate.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 0.995
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 1.005
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.005
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 1.005
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1.01
Hay más en la clase 18… pero paja.
Modelar la presencia/ausencia de depósitos auríferos a partir de la [Sb] y la cercanía a falla geológica.
Estudio observacional.
VR = presencia/ausencia de depósitos auríferos (dicotómica, ~ Bernoulli) VE = [Sb] cuantitativa continua VE = cercanía a falla geológica, categórica, 2 niveles
Modelo completo en escala PL:
\[ logit(\pi_i) = \eta_i = \beta_0 + \beta_1 Sb_i + \beta_2 Falla_i + \beta_3 Sb_i Falla_i \]
En escala odds:
\[ odds_i = \frac{\pi_i}{1 - \pi_i} = e^{\beta_0 + \beta_1 Sb_i + \beta_2 Falla_i + \beta_3 Sb_i Falla_i} \]
En la escala de la VR:
\[ E[Y_i] = \pi_i = \frac{e^{\eta_i}}{1 + e^{\eta_i}} \]
Odds: cociente entre proba de depósito de oro y proba de no depósito.
oro <- read_table2("tps/TP 8/oro.txt")
## Parsed with column specification:
## cols(
## Sb = col_double(),
## falla_geologica = col_double(),
## deposito_oro = col_double()
## )
oro$falla_geologica <- factor(oro$falla_geologica)
GGally::ggpairs(oro, progress = F)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Problema: Sb
es muy asimétrica cuando no hay falla geológica.
ggplot(oro) +
aes(falla_geologica, Sb) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(position = position_jitter(width = .05)) +
coord_flip()
ggplot(oro) +
aes(falla_geologica, log(Sb)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(position = position_jitter(width = .05)) +
coord_flip()
oro %>%
filter(Sb < 15) %>%
ggplot() +
aes(falla_geologica, Sb) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(position = position_jitter(width = .05)) +
coord_flip()
oro %>%
mutate(
Sb_d = cut(Sb, breaks = c(0, 2.5, 20), labels = c("< 2.5", "> 2.5"))
)
## # A tibble: 64 x 4
## Sb falla_geologica deposito_oro Sb_d
## <dbl> <fct> <dbl> <fct>
## 1 3.08 1 1 > 2.5
## 2 6.15 1 1 > 2.5
## 3 2.35 1 1 < 2.5
## 4 0.3 0 0 < 2.5
## 5 0.3 0 0 < 2.5
## 6 9.62 1 1 > 2.5
## 7 0.51 0 0 < 2.5
## 8 3.71 1 1 > 2.5
## 9 4.32 0 0 > 2.5
## 10 0.8 0 0 < 2.5
## # … with 54 more rows
oro_2 <- filter(oro, Sb < 15)
oro_2 <- oro # Ellos no quitaron el outlier al final...
oro_glm <- glm(deposito_oro ~ Sb * falla_geologica, data = oro_2, family = binomial)
drop1(oro_glm, test = "Chisq")
## Single term deletions
##
## Model:
## deposito_oro ~ Sb * falla_geologica
## Df Deviance AIC LRT Pr(>Chi)
## <none> 33.427 41.427
## Sb:falla_geologica 1 35.529 41.529 2.1016 0.1471
Podemos quitar la interacción!
oro_glm2 <- glm(deposito_oro ~ Sb + falla_geologica, data = oro_2, family = binomial)
drop1(oro_glm2, test = "Chisq")
## Single term deletions
##
## Model:
## deposito_oro ~ Sb + falla_geologica
## Df Deviance AIC LRT Pr(>Chi)
## <none> 35.529 41.529
## Sb 1 56.607 60.607 21.078 4.410e-06 ***
## falla_geologica 1 69.124 73.124 33.595 6.787e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No quitamos las variables del modelo aditivo.
En una diapo dicen que como la variable es categórica “no hay posibilidad de sobre o subdispersión”, en Bernoulli. … Hmmm.
Se hace un ajuste OBS vs ESP con Chi cuadrado. La H0 es que hay buen ajuste
ResourceSelection::hoslem.test(oro_2$deposito_oro, fitted(oro_glm2))
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: oro_2$deposito_oro, fitted(oro_glm2)
## X-squared = 4.3532, df = 8, p-value = 0.8239
betas <- coefficients(oro_glm2) %>% round(2)
betas
## (Intercept) Sb falla_geologica1
## -5.47 0.71 5.20
En caso de falla geológica, el predictor lineal es:
logit = eta_i = (5.2 - 5.47) + 0.71 * Sb_i = -0.27 + 0.71 * Sb_i
En caso de no falla geológica, es:
logit = eta_i = -5.47 + 0.71 * Sb_i
En ambos casos, significa que por cada aumento de unidad de [Sb], el log odds aumenta en 0.71. En caso de que además haya una falla geológica, si se mantiene la [Sb] constante, el log odds aumenta en 5.2.
exp(betas) %>% round(2)
## (Intercept) Sb falla_geologica1
## 0.00 2.03 181.27
Tanto en presencia como en ausencia de una falla geológica, por cada unidad de [Sb], los odds de tener un depósito acuífero crecen un 103% (son 2.03 veces los odds con una unidad menos de [Sb]). En el caso de que haya una falla geológica, los odds de tener un depósito crecen un 1800% (son 181 veces los odds sin la falla), manteniendo la [Sb] constante.
plot_model(oro_glm2, show.values = T)
## Waiting for profiling to be done...
# tab_model(oro_glm2)
La presencia de una falla geológica
Probabilidad de encontrar oro en un sitio sin fall y con 6 ppm [Sb]:
Falla_i <- 0
Sb_i <- 6
eta_i <- betas["(Intercept)"] + betas["Sb"] * Sb_i + betas["falla_geologica1"] * Falla_i
log_odds_i <- eta_i
proba_i <- exp(eta_i)/(1 - exp(eta_i))
proba_i[[1]]
## [1] 0.4249019
oro_2$pred_logit <- predict(oro_glm2, type = "link")
oro_2$pred_odds <- exp(oro_2$pred_logit)
oro_2$pred_prob <- predict(oro_glm2, type = "response")
oro_2$falla_geologica.f <- as.factor(oro_2$falla_geologica)
ggplot(oro_2) +
aes(Sb, pred_logit, group = falla_geologica.f, color = falla_geologica.f) +
geom_line() +
geom_point()
ggplot(oro_2) +
aes(Sb, pred_odds, color = falla_geologica.f, group = falla_geologica.f) +
geom_smooth(
se = F,
method = "glm",
method.args = list(
family = gaussian(link = "log")
)
) +
geom_point()
oro_2$pred_class <- oro_2$pred_prob >= .5
oro_2$pred_success <- oro_2$pred_class == oro_2$deposito_oro
ggplot(oro_2) +
aes(Sb, pred_prob, color = pred_success, group = falla_geologica.f) +
geom_smooth(
method = "glm",
method.args = list(
family = gaussian(link = "logit")
)
) +
geom_point(size = 4) +
geom_hline(yintercept = 0.5) # +
# scale_shape_manual(values = c(16, 1)) +
# scale_color_manual(values = c("Black", "DarkOrange"))
rm(list = ls())
evaluar si una vez que comience la implementación masiva podrán distinguir mediante imágenes satelitales aquellos campos donde se utilice el nuevo método de aquellas donde se haga riego por goteo la reflectancia en las tres bandas dentro del espectro visible (AZUL, VERDE y ROJO) en 11 campos con riego nuevo y en otros 11 campos del mismo tamaño y con el mismo cultivo, pero bajo riego por goteo
land <- read_table2("tps/TP 8/Landsat.txt")
## Parsed with column specification:
## cols(
## Riego = col_character(),
## Rojo = col_double(),
## Verde = col_double(),
## Azul = col_double()
## )
land
## # A tibble: 22 x 4
## Riego Rojo Verde Azul
## <chr> <dbl> <dbl> <dbl>
## 1 nuevo 10.8 19 10.3
## 2 nuevo 19.2 18.2 11.5
## 3 nuevo 5.4 6.6 8.7
## 4 nuevo 1.4 10 15.1
## 5 nuevo 5 12 17
## 6 nuevo 9 10 20.5
## 7 nuevo 21.8 37.6 11.9
## 8 nuevo 24 33 12.6
## 9 nuevo 21.8 17.8 10.9
## 10 nuevo 16.6 12 7.5
## # … with 12 more rows
library(GGally)
land$Riego <- as.factor(land$Riego)
ggpairs(
land,
progress = F,
lower = list(continuous = wrap("smooth", method = "lm", se = F))
)
## `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`.
- Indique el tipo de estudio, la variable respuesta (tipo y potencial distribución de probabilidades). Describa el tipo de variables explicatorias. Analice la asociación entre variables.
La correlación no es tan alta, la máxima es entre Rojo y Verde, 0.414.
- Escriba el modelo lineal generalizado a utilizar en términos del problema (modelo aditivo completo), explicitando los tres componentes del GLM. Interprete el odds en el contexto del problema.
Modelo aditivo completo:
$$ log ( ) = _i = _0 + _1 Nuevo_i + _2 Rojo_i + _3 Verde_i + _4 Azul_i \
odds_i = e^{_i} \
E[Y_i] = _i = $$
En este problema, el odds es el cociente entre la probabilidad de que un campo tenga nuevo riego y la de que tenga riego por goteo.
- Genere un modelo de regresión para la identificación del tipo de riego a partir de las variables explicatorias seleccionadas, realizando un procedimiento de selección “hacia adelante”.
formula_completa <- as.formula("Riego ~ Rojo + Verde + Azul")
land_glm0 <- glm(Riego ~ 1, data = land, family = binomial)
add1(land_glm0, formula_completa, test = "Chisq")
## Single term additions
##
## Model:
## Riego ~ 1
## Df Deviance AIC LRT Pr(>Chi)
## <none> 30.498 32.498
## Rojo 1 27.150 31.150 3.3488 0.06725 .
## Verde 1 27.530 31.530 2.9686 0.08490 .
## Azul 1 29.319 33.319 1.1796 0.27744
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
land_glm1 <- update(land_glm0, . ~ . + Verde)
add1(land_glm1, formula_completa, test = "Chisq")
## Single term additions
##
## Model:
## Riego ~ Verde
## Df Deviance AIC LRT Pr(>Chi)
## <none> 27.530 31.530
## Rojo 1 16.966 22.966 10.564 0.001153 **
## Azul 1 26.806 32.806 0.724 0.394831
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
land_glm2 <- update(land_glm1, . ~ . + Rojo)
add1(land_glm2, formula_completa, test = "Chisq")
## Single term additions
##
## Model:
## Riego ~ Verde + Rojo
## Df Deviance AIC LRT Pr(>Chi)
## <none> 16.966 22.966
## Azul 1 16.961 24.961 0.0051112 0.943
Como vemos, al variable azul no aporta. Nos quedamos con el modelo:
\[ \eta_i = \beta_0 + \beta_1 Nuevo_i + \beta_2 Rojo_i + \beta_3 Verde_i \]
Confirmamos que no hay colinealidad preocupante:
car::vif(land_glm2)
## Verde Rojo
## 2.222361 2.222361
- Analice el modelo obtenido y los supuestos asociados.
summary(land_glm2)
##
## Call:
## glm(formula = Riego ~ Verde + Rojo, family = binomial, data = land)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.22049 -0.52791 0.01653 0.61002 1.47083
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9133 1.7526 0.521 0.6023
## Verde 0.3126 0.1618 1.932 0.0533 .
## Rojo -0.3212 0.1400 -2.294 0.0218 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 30.498 on 21 degrees of freedom
## Residual deviance: 16.966 on 19 degrees of freedom
## AIC: 22.966
##
## Number of Fisher Scoring iterations: 6
land$predichos_prob <- predict(land_glm2, type = "response")
land$predichos_logit <- predict(land_glm2, type = "link")
land$predichos_odds <- exp(land$predichos_logit)
land
## # A tibble: 22 x 7
## Riego Rojo Verde Azul predichos_prob predichos_logit predichos_odds
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 nuevo 10.8 19 10.3 0.967 3.38 29.5
## 2 nuevo 19.2 18.2 11.5 0.607 0.435 1.55
## 3 nuevo 5.4 6.6 8.7 0.776 1.24 3.46
## 4 nuevo 1.4 10 15.1 0.973 3.59 36.2
## 5 nuevo 5 12 17 0.955 3.06 21.3
## 6 nuevo 9 10 20.5 0.759 1.15 3.15
## 7 nuevo 21.8 37.6 11.9 0.997 5.66 289.
## 8 nuevo 24 33 12.6 0.971 3.52 33.8
## 9 nuevo 21.8 17.8 10.9 0.372 -0.525 0.592
## 10 nuevo 16.6 12 7.5 0.339 -0.668 0.513
## # … with 12 more rows
- ¿Cómo interpreta la magnitud del efecto para cada una de las variables explicatorias seleccionadas?c. Escriba la ecuación estimada en la escala del predictor lineal para el modelo obtenido.
Ecuación estimada en escala del PL:
eta_i = 0.91 + 0.31 Verde_i - 0.32 Rojo_i
Ecuación estimada en escala de los odds:
odds_i = e^{0.91 + 0.31 Verde_i - 0.32 Rojo_i}
Probabilidad estimada:
pi_i = e^{eta_i} / (1 - e^{eta_i})
Magnitud del efecto:
coef(land_glm2) %>% round(2)
## (Intercept) Verde Rojo
## 0.91 0.31 -0.32
confint(land_glm2) %>% round(2)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) -2.59 4.76
## Verde 0.09 0.76
## Rojo -0.69 -0.10
En la escala de los odds:
exp(coef(land_glm2)) %>% round(2)
## (Intercept) Verde Rojo
## 2.49 1.37 0.73
exp(confint(land_glm2)) %>% round(2)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 0.08 117.18
## Verde 1.09 2.15
## Rojo 0.50 0.90
Grafique la probabilidad de ocurrencia del sistema de goteo nuevo en función de las variables explicatorias.
ggplot(land) +
aes(Rojo, predichos_prob) +
geom_point()
ggplot(land) +
aes(Verde, predichos_prob) +
geom_point()
ggplot(land) +
aes(Verde, predichos_prob) +
geom_point(aes(x = Verde), color = "Green4") +
geom_point(aes(x = Rojo), color = "Tomato") +
geom_smooth(aes(x = Verde), color = "Green4", se = F, method = "lm") +
geom_smooth(aes(x = Rojo), color = "Tomato", se = F, method = "lm") +
labs(x = "Reflectancia")
- Concluya en términos del problema e Informe sus resultados.
El objetivo en este caso es predecir la probabilidad de ocurrencia de sobre-regularización con respecto a la edad del niño, y si existen diferencias de adquisición según el sexo.
leng <- read_table2("tps/TP 8/lenguaje.txt")
## Parsed with column specification:
## cols(
## edad_meses = col_double(),
## sexo = col_character(),
## N_verbos_pasado = col_double(),
## N_sobrereg = col_double()
## )
leng
## # A tibble: 92 x 4
## edad_meses sexo N_verbos_pasado N_sobrereg
## <dbl> <chr> <dbl> <dbl>
## 1 56 F 90 6
## 2 51 F 56 4
## 3 44 F 32 4
## 4 32 M 44 10
## 5 55 F 71 8
## 6 44 F 15 2
## 7 14 M 3 1
## 8 14 M 10 5
## 9 37 M 60 9
## 10 33 M 45 5
## # … with 82 more rows
- Indique cuál es la variable respuesta y las variables explicatorias. ¿Qué valores toman y cómo es su distribución?
- Construya a partir de los datos una variable que sea la proporción de errores de sobre-regularización sobre el total de verbos usados en tiempo pasado.
leng$proporcion <- leng$N_sobrereg/leng$N_verbos_pasado
leng$obs_id <- 1:nrow(leng)
leng
## # A tibble: 92 x 6
## edad_meses sexo N_verbos_pasado N_sobrereg proporcion obs_id
## <dbl> <chr> <dbl> <dbl> <dbl> <int>
## 1 56 F 90 6 0.0667 1
## 2 51 F 56 4 0.0714 2
## 3 44 F 32 4 0.125 3
## 4 32 M 44 10 0.227 4
## 5 55 F 71 8 0.113 5
## 6 44 F 15 2 0.133 6
## 7 14 M 3 1 0.333 7
## 8 14 M 10 5 0.5 8
## 9 37 M 60 9 0.15 9
## 10 33 M 45 5 0.111 10
## # … with 82 more rows
Parece haber un outlier: el pibe habla MUY MAL! Lo quitamos.
ggplot(leng) +
aes(x = edad_meses, y = proporcion, color = sexo) +
geom_point() +
geom_text(
data = filter(leng, proporcion > .5 & edad_meses > 50),
mapping = aes(label = obs_id),
nudge_y = -.05
)
leng2 <- filter(leng, obs_id != 68)
- Genere un modelo de regresión de la ocurrencia de errores de sobre- regularización en función de la edad.
\[ \log odds_i = \eta_i = \beta_0 + \beta_1 Meses_i \]
leng_glm <- glm(proporcion ~ edad_meses, family = binomial, weights = N_verbos_pasado, data = leng2)
- Analice el modelo obtenido, los supuestos asociados, y grafique el modelo obtenido en función de la variable explicatoria. ¿Cuánto vale el odds para la misma?
leng2$pred_logit <- predict(leng_glm, type = "link")
leng2$pred_odds <- exp(leng2$pred_logit)
leng2$pred_prob <- predict(leng_glm, type = "response")
ggplot(leng2) +
aes(edad_meses, proporcion) +
geom_segment(aes(xend = edad_meses, yend = pred_prob), color = "Grey70") +
geom_point(aes(color = sexo)) +
# geom_point(aes(y = pred_prob), shape = 4) +
geom_line(aes(y = pred_prob))
No hay sobre/subdispersión:
rp <- residuals(leng_glm, type = "pearson")
gl <- leng_glm$df.residual
sum(rp^2)/gl
## [1] 0.9872394
Supuestos:
leng2$residuos_p <- rp
ggplot(leng2) +
aes(pred_prob, residuos_p) +
geom_hline(yintercept = 0) +
geom_point(aes(color = sexo))
leng2$cooks_D <- cooks.distance(leng_glm)
ggplot(leng2) +
aes(obs_id, cooks_D) +
geom_point()
- Escriba la ecuación del modelo obtenido.
coefficients(leng_glm) %>% round(2)
## (Intercept) edad_meses
## -0.56 -0.04
log odds i = eta_i = - 0.56 - 0.04 Edad_i
proporcion esperada = e^eta_i / (1 + e ^ eta_i)
- Incorpore al modelo la información sobre el sexo de la persona. Analice la presencia de interacción en el modelo. Grafique y escriba el modelo obtenido para ambos sexos.
leng_glm2 <- glm(proporcion ~ edad_meses * sexo,
data = leng2, family = binomial, weights = N_verbos_pasado)
drop1(leng_glm2, test = "Chisq")
## Single term deletions
##
## Model:
## proporcion ~ edad_meses * sexo
## Df Deviance AIC LRT Pr(>Chi)
## <none> 92.482 388.68
## edad_meses:sexo 1 92.787 386.99 0.3049 0.5808
No hay interacción.
leng_glm3 <- glm(proporcion ~ edad_meses + sexo,
data = leng2, family = binomial, weights = N_verbos_pasado)
drop1(leng_glm3, test = "Chisq")
## Single term deletions
##
## Model:
## proporcion ~ edad_meses + sexo
## Df Deviance AIC LRT Pr(>Chi)
## <none> 92.787 386.99
## edad_meses 1 188.308 480.51 95.521 <2e-16 ***
## sexo 1 93.036 385.24 0.249 0.6178
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El sexo no explica la adquisición del lenguaje.
Nos quedamos con el primer modelo:
coefficients(leng_glm) %>% round(2)
## (Intercept) edad_meses
## -0.56 -0.04
confint(leng_glm) %>% round(2)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) -0.86 -0.26
## edad_meses -0.05 -0.03
eta_i = -0.56 - 0.04 Edad_i
exp(coefficients(leng_glm)) %>% round(2)
## (Intercept) edad_meses
## 0.57 0.96
exp(confint(leng_glm)) %>% round(2)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 0.42 0.77
## edad_meses 0.96 0.97
sjPlot::plot_model(leng_glm)
## Waiting for profiling to be done...
sjPlot::plot_model(leng_glm, type = "pred", terms = "edad_meses [all]")
- Concluya en términos del problema e Informe sus resultados.