Os alunos querem se sentir seguros e protegidos ao frequentar uma faculdade ou universidade. Em resposta à legislação, o Departamento de Educação dos Estados Unidos busca fornecer dados e garantias para alunos e pais. Todas as instituições pós-secundárias que participam de programas federais de auxílio estudantil são exigidas pela Jeanne Clery Disclosure of Campus Security Policy e Campus Crime Statistics Act e pela Higher Education Opportunity Act para coletar e relatar dados sobre crimes ocorridos no campus para o Departamento de Educação. Por sua vez, esses dados estão disponíveis publicamente no site da Secretaria de Educação Superior. Estamos interessados em verificar se existem diferenças regionais nos crimes violentos no campus, controlando as diferenças no tipo de escola.
input = ("
Enrollment type nv nvrate enroll1000 region
5590 U 30 5.366726296958855 5.59 SE
540 C 0 0 0.54 SE
35747 U 23 0.6434106358575545 35.747 W
28176 C 1 0.03549119818285065 28.176 W
10568 U 1 0.09462528387585163 10.568 SW
3127 U 0 0 3.127 SW
20675 U 7 0.3385731559854897 20.675 W
12548 C 0 0 12.548 W
30063 U 19 0.6320061204803247 30.063 C
4429 C 4 0.9031384059607135 4.429 C
12254 U 3 0.2448180186061694 12.254 NE
11136 C 6 0.5387931034482758 11.136 NE
20949 U 26 1.2411093608286792 20.949 NE
12316 U 23 1.8674894446248782 12.316 SE
10678 C 7 0.6555534744334145 10.678 SE
25743 U 15 0.5826826710173639 25.743 SE
3641 C 0 0 3.641 SE
11321 U 9 0.7949827753732003 11.321 MW
1260 C 0 0 1.26 MW
27823 U 11 0.3953563598461704 27.823 C
14410 U 2 0.13879250520471895 14.41 C
22396 U 7 0.31255581353813183 22.396 C
6723 U 1 0.14874312063067083 6.723 C
23901 U 8 0.33471402870172795 23.901 SE
15489 U 13 0.8393053134482537 15.489 SW
13404 C 0 0 13.404 SW
10698 U 1 0.09347541596560105 10.698 NE
10966 U 5 0.4559547692868868 10.966 NE
11237 U 0 0 11.237 NE
5181 C 19 3.667245705462266 5.181 NE
13348 U 1 0.0749175906502847 13.348 NE
4339 C 6 1.3828070984097718 4.339 NE
38248 U 18 0.4706128425015687 38.248 MW
9764 C 0 0 9.764 MW
46597 U 9 0.1931454814687641 46.597 MW
9380 U 2 0.21321961620469082 9.38 MW
16878 U 4 0.2369949046095509 16.878 SE
1318 C 1 0.7587253414264037 1.318 SE
12969 U 2 0.15421389467190993 12.969 C
2878 C 3 1.0423905489923557 2.878 C
22764 U 3 0.13178703215603585 22.764 C
6426 U 0 0 6.426 C
23313 U 13 0.5576287908034144 23.313 W
9697 C 0 0 9.697 W
35650 U 24 0.6732117812061711 35.65 NE
6846 C 7 1.0224948875255624 6.846 NE
23753 U 27 1.136698522291921 23.753 SW
3556 U 3 0.843644544431946 3.556 SW
4018 U 1 0.24888003982080636 4.018 NE
6483 U 1 0.15424957581366652 6.483 NE
25494 U 8 0.31379932533145055 25.494 SE
1673 C 0 0 1.673 SE
11764 U 5 0.42502550153009183 11.764 C
2292 C 0 0 2.292 C
18739 U 3 0.1600939217674369 18.739 MW
20496 C 1 0.04879000780640125 20.496 MW
25104 U 7 0.2788400254939452 25.104 SW
1961 C 0 0 1.961 SW
3239 U 0 0 3.239 NE
1902 C 1 0.5257623554153522 1.902 NE
14264 U 4 0.28042624789680315 14.264 NE
7774 C 5 0.6431695394906097 7.774 NE
23000 U 17 0.7391304347826088 23 SE
1476 U 2 1.3550135501355014 1.476 SE
9260 U 1 0.1079913606911447 9.26 C
26033 U 8 0.3073022701955211 26.033 SE
12736 C 1 0.0785175879396985 12.736 SE
7369 U 1 0.13570362328674176 7.369 SW
13313 U 1 0.07511454968827462 13.313 SW
27668 U 6 0.21685701893884632 27.668 C
22609 C 1 0.044230173824583136 22.609 C
4638 U 10 2.1561017680034498 4.638 NE
4114 C 1 0.24307243558580455 4.114 NE
21073 U 7 0.33217861718787073 21.073 W
4227 C 1 0.2365744026496333 4.227 W
22774 U 1 0.04390972161236498 22.774 NE
3072 C 3 0.9765625 3.072 NE
40922 U 13 0.3176775328674063 40.922 MW
10549 U 7 0.6635700066357001 10.549 MW
12366 U 0 0 12.366 C
2729 C 0 0 2.729 C
")
Dados = read.table(textConnection(input),header=TRUE)
str(Dados)
## 'data.frame': 81 obs. of 6 variables:
## $ Enrollment: int 5590 540 35747 28176 10568 3127 20675 12548 30063 4429 ...
## $ type : chr "U" "C" "U" "C" ...
## $ nv : int 30 0 23 1 1 0 7 0 19 4 ...
## $ nvrate : num 5.3667 0 0.6434 0.0355 0.0946 ...
## $ enroll1000: num 5.59 0.54 35.75 28.18 10.57 ...
## $ region : chr "SE" "SE" "W" "W" ...
Os fatores foram ordenados pela ordem no arquivo de dados, Caso contrário, o R irá colocá-los em ordem alfabética.
Dados$type = factor(Dados$type, levels = unique(Dados$type))
Dados$region = factor(Dados$region, levels = unique(Dados$region))
head(Dados)
## Enrollment type nv nvrate enroll1000 region
## 1 5590 U 30 5.36672630 5.590 SE
## 2 540 C 0 0.00000000 0.540 SE
## 3 35747 U 23 0.64341064 35.747 W
## 4 28176 C 1 0.03549120 28.176 W
## 5 10568 U 1 0.09462528 10.568 SW
## 6 3127 U 0 0.00000000 3.127 SW
Estes dados foram utlizados no livro Beyond Multiple Linear Regression: Applied Generalized Linear Models and Multilevel Models in R. Autores Paul Roback and Julie Legler, livro publicado em 2021.
Cada linha do arquivo de dados contém informações sobre crimes de uma instituição pós-secundária, uma faculdade ou universidade.
As variáveis incluem:
Um gráfico do número de crimes violentos, na figura abaixo, revela o padrão freqüentemente encontrado com distribuições de contagens de eventos raros. Muitas escolas relataram nenhum crime violento ou muito poucos crimes. Algumas escolas têm um grande número de crimes, o que significa uma distribuição que parece longe do normal. Portanto, a regressão Poisson deve ser usada para modelar nossos dados; Variáveis aleatórias Poisson são frequentemente usadas para representar contagens (por exemplo, número de crimes violentos) por unidade de tempo ou espaço (por exemplo, um ano).
library(lattice)
histogram( ~ nv | type, data = Dados, layout=c(2,1))
histogram( ~ nv | region, data = Dados, layout=c(6,1))
Vamos dar uma olhada em duas covariáveis de interesse para essas escolas: tipo de instituição e região. Em nossos dados, a maioria das instituições são universidades (65% das 81 escolas) e apenas 35% são faculdades. Centra o interesse em saber se as diferentes regiões tendem a ter diferentes taxas de criminalidade. A tabela a seguir contém o nome de cada região e cada coluna representa a porcentagem de escolas naquela região que são faculdades ou universidades. A proporção de faculdades varia de um mínimo de 20% no sudoeste (SW) a um máximo de 50% no oeste (W).
tabela = xtabs( ~ type + region, data = Dados)
library(gmodels)
CrossTable(tabela, prop.r = FALSE, prop.t = FALSE, prop.chisq = FALSE)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Col Total |
## |-------------------------|
##
##
## Total Observations in Table: 81
##
##
## | region
## type | SE | W | SW | C | NE | MW | Row Total |
## -------------|-----------|-----------|-----------|-----------|-----------|-----------|-----------|
## U | 9 | 4 | 8 | 12 | 13 | 7 | 53 |
## | 0.600 | 0.500 | 0.800 | 0.706 | 0.619 | 0.700 | |
## -------------|-----------|-----------|-----------|-----------|-----------|-----------|-----------|
## C | 6 | 4 | 2 | 5 | 8 | 3 | 28 |
## | 0.400 | 0.500 | 0.200 | 0.294 | 0.381 | 0.300 | |
## -------------|-----------|-----------|-----------|-----------|-----------|-----------|-----------|
## Column Total | 15 | 8 | 10 | 17 | 21 | 10 | 81 |
## | 0.185 | 0.099 | 0.123 | 0.210 | 0.259 | 0.123 | |
## -------------|-----------|-----------|-----------|-----------|-----------|-----------|-----------|
##
##
# Proporção de faculdades e universidades dentro da região no conjunto de dados de crimes do campus.
Embora um modelo de regressão de Poisson seja uma boa primeira escolha, porque as respostas são contagens por ano, é importante notar que as contagens não são diretamente comparáveis; porque vêm de escolas de tamanhos diferentes. Esse problema às vezes é referido como a necessidade de contabilizar o esforço de amostragem; em outras palavras, esperamos que escolas com mais alunos tenham mais denúncias de crimes violentos, uma vez que há mais alunos que poderiam ser afetados.
Não podemos comparar diretamente os 30 crimes violentos da primeira escola no conjunto de dados com nenhum crime violento para a segunda escola quando suas matrículas são muito diferentes: 5.590 para a escola 1 contra 540 para a escola 2. Podemos levar em consideração as diferenças nas matrículas incluindo um deslocamento em nosso modelo, que discutiremos na próxima seção. Para o restante do estudo descritivo, examinamos as contagens de crimes violentos em termos da taxa por 1.000 inscritos, ou seja \[\begin{equation*} \frac{\mbox{número de crimes violentos}}{\mbox{número de inscritos}}\times 1000\cdot \end{equation*}\]
Observe que há um outlier perceptível para uma escola do sudeste (5.4 crimes violentos por 1000 alunos) e há uma taxa observada de 0 para as faculdades do sudoeste que pode levar a alguns problemas computacionais. Portanto, combinamos SW e SE para formar uma única categoria do Sul e também removemos a observação extrema do conjunto de dados.
levels(Dados$region)[match(c("SE","SW"),levels(Dados$region))] <- "S"
Dados = Dados[-1,]
library(doBy)
summaryBy(nvrate ~ region + type, data=Dados, FUN=c(mean, var))
## region type nvrate.mean nvrate.var
## 1 S U 0.57131616 0.2778065495
## 2 S C 0.18659955 0.1047177679
## 3 W U 0.46794780 0.0246669855
## 4 W C 0.06801640 0.0129073834
## 5 C U 0.22194406 0.0349266302
## 6 C C 0.39795183 0.2780912668
## 7 NE U 0.43592725 0.3850333363
## 8 NE C 1.12498845 1.1820999739
## 9 MW U 0.40190031 0.0620748425
## 10 MW C 0.01626334 0.0007934883
A tabelas acima mostram as taxas médias de crimes violentos que geralmente são mais baixas nas faculdades de uma região (com exceção do Nordeste). Além disso, o padrão regional das taxas nas universidades parece ser diferente do das faculdades.
library(ggplot2)
ggplot(Dados, aes(x=region, y=nvrate, fill=type)) + geom_boxplot()
Embora trabalhar com as taxas observadas (por 1000 alunos) seja útil durante a análise exploratória dos dados, não usamos essas taxas explicitamente no modelo. As contagens (por ano) são as respostas de Poisson na modelagem, portanto, devemos levar em consideração a matrícula de uma forma diferente.
Nossa abordagem é incluir um termo no lado direito do modelo chamado de offset (deslocamento), que é o logaritmo da matrícula, em milhares. Existe uma heurística intuitiva para a forma do offset. Se pensarmos em \(\lambda\) como o número médio de crimes violentos por ano, então \(\lambda/\mbox{enroll1000}\) representa o número por 1000 alunos, de forma que a contagem anual seja ajustada para ser comparável em escolas de tamanhos diferentes.
Ajustar a contagem anual por matrícula é equivalente a adicionar \(\log(enroll000)\) ao lado direito da equação de regressão Poisson - essencialmente adicionando um preditor com um coeficiente fixo de 1: \[\begin{equation} \begin{array}{rcl} \log\left(\frac{\lambda}{enroll1000}\right) & = & \beta_0+\beta_1 \mbox{type} \\ \log(\lambda)-\log(enroll1000) & = & \beta_0+\beta_1 \mbox{type} \\ \log(\lambda) & = & \beta_0+\beta_1 \mbox{type}+\log(enroll1000)\cdot \end{array} \end{equation}\]
Embora essa heurística seja útil, é importante observar que não é \(\lambda/enroll1000\) que estamos modelando. Ainda estamos modelando \(\log(\lambda)\), mas estamos adicionando um deslocamento para ajustar para matrículas diferentes, onde o deslocamento tem a característica incomum de que o coeficiente é fixado em 1,0. Como resultado, nenhum coeficiente estimado para enroll1000 ou \(\log(enroll1000)\) aparecerá na saída. Isto ilustra que a modelagem do \(\log(\lambda)\) e a adição de um offset são equivalentes à modelagem de taxas e os coeficientes podem ser interpretados dessa forma.
Na tabela acima vemos que variâncias são muito maiores do que as contagens médias em alguns grupos. Assim, temos motivos para questionar a suposição da regressão Poisson de variabilidade igual à média; teremos que retornar a esse problema após alguma modelagem inicial. O fato de que a variação da taxa de crimes violentos por 1000 alunos tende a estar na mesma escala que a média nos diz que o ajuste para matrícula pode fornecer alguma ajuda, embora isso possa não resolver completamente nossos problemas com a variação excessiva.
Quanto a outras suposições do modelo, a linearidade com relação ao \(\log(\lambda)\) é difícil de discernir sem preditores contínuos e não é possível avaliar a independência sem saber como as escolas foram selecionadas.
Estamos interessados principalmente nas diferenças em crimes violentos entre os tipos institucionais que controlam as diferenças nas regiões, portanto, ajustamos um modelo com a região, tipo institucional e nossa compensação. Observe que a região central é o nível de referência em nosso modelo.
modeltr <- glm(nv ~ type + region, family = poisson(link = "log"), offset = log(enroll1000), data = Dados)
summary(modeltr)
##
## Call:
## glm(formula = nv ~ type + region, family = poisson(link = "log"),
## data = Dados, offset = log(enroll1000))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.5420 -1.7772 -0.7457 1.0884 6.7398
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.68586 0.08399 -8.166 3.18e-16 ***
## typeC -0.27956 0.13314 -2.100 0.03576 *
## regionW -0.31963 0.16296 -1.961 0.04983 *
## regionC -0.58238 0.14896 -3.910 9.24e-05 ***
## regionNE 0.19574 0.12182 1.607 0.10808
## regionMW -0.48327 0.15143 -3.191 0.00142 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 392.76 on 79 degrees of freedom
## Residual deviance: 348.68 on 74 degrees of freedom
## AIC: 573.32
##
## Number of Fisher Scoring iterations: 6
Em nosso modelo, as regiões W (Oeste), C (Central) e MW (Centro-Oeste) diferem significativamente da região Sul (p = 0.04983, p = 0.0000924 e p = 0.00142, respectivamente). O coeficiente estimado de -0.31963 significa que a taxa de crimes violentos por 1.000 no Oeste é quase 0.7 (\(\exp(-0.31963)\)) vezes a da região Sul que controla o tipo de escola.
exp(confint(modeltr))
## 2.5 % 97.5 %
## (Intercept) 0.4253319 0.5912987
## typeC 0.5777879 0.9745118
## regionW 0.5233041 0.9926145
## regionC 0.4146207 0.7442309
## regionNE 0.9568936 1.5434409
## regionMW 0.4553725 0.8254041
Um intervalo de confiança do tipo Wald para esse fator pode ser construído calculando primeiro um IC para o coeficiente (\(-0.31963 \pm 1,96\times 0.16296\)) e depois exponenciando (0.52 a 0.99).
As comparações com regiões diferentes da região Sul podem ser realizadas alterando a região de referência. Se muitas comparações forem feitas, seria melhor ajustar para comparações múltiplas usando um método como Tukey’s Honestly Significant Differences, que considera todas as comparações de pares entre regiões. Este método ajuda a controlar o grande número de falsos positivos que veríamos se executássemos vários grupos de comparação de testes t. A diferença honestamente significativa compara uma diferença média padronizada entre dois grupos com um valor crítico de uma distribuição de intervalo estudentizada.
library(multcomp)
mult_comp <- summary(glht(modeltr, mcp(region="Tukey")))
mult_comp
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = nv ~ type + region, family = poisson(link = "log"),
## data = Dados, offset = log(enroll1000))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## W - S == 0 -0.31963 0.16296 -1.961 0.2795
## C - S == 0 -0.58238 0.14896 -3.910 <0.001 ***
## NE - S == 0 0.19574 0.12182 1.607 0.4863
## MW - S == 0 -0.48327 0.15143 -3.191 0.0121 *
## C - W == 0 -0.26275 0.18753 -1.401 0.6209
## NE - W == 0 0.51537 0.16587 3.107 0.0156 *
## MW - W == 0 -0.16364 0.18942 -0.864 0.9079
## NE - C == 0 0.77813 0.15307 5.084 <0.001 ***
## MW - C == 0 0.09912 0.17752 0.558 0.9804
## MW - NE == 0 -0.67901 0.15545 -4.368 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
Em nosso caso, as Diferenças Honestamente Significativas de Tukey avaliam simultaneamente todas as 10 diferenças médias entre pares de regiões. Descobrimos que o Nordeste tem taxas significativamente mais altas de crimes violentos do que as regiões Centro, Centro-Oeste e Oeste, enquanto o Sul tem taxas significativamente mais altas de crimes violentos do que o Centro e o Centro-Oeste, controlando pelo tipo de instituição. No modelo primário, o indicador da Universidade é significativo e, após exponenciar o coeficiente, pode ser interpretado como um aumento de aproximadamente (\(\exp(0.28)\)) 32% na taxa de crimes violentos em relação às faculdades após o controle da região.
Esses resultados certamente sugerem diferenças significativas nas regiões e no tipo de instituição. No entanto, os resultados da EDA sugerem que o efeito do tipo de instituição pode variar dependendo da região, então consideramos um modelo com uma interação entre região e tipo.
modeli <- glm(nv ~ type + region + region:type, family = poisson(link = "log"),
offset = log(enroll1000), data = Dados)
summary(modeli)
##
## Call:
## glm(formula = nv ~ type + region + region:type, family = poisson(link = "log"),
## data = Dados, offset = log(enroll1000))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7292 -1.5701 -0.8148 0.9364 4.9575
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.62229 0.08513 -7.310 2.67e-13 ***
## typeC -1.00806 0.34403 -2.930 0.00339 **
## regionW -0.07891 0.16506 -0.478 0.63263
## regionC -0.65587 0.15745 -4.166 3.11e-05 ***
## regionNE -0.17275 0.14224 -1.214 0.22458
## regionMW -0.43593 0.15375 -2.835 0.00458 **
## typeC:regionW -1.59851 0.79897 -2.001 0.04542 *
## typeC:regionC 0.81212 0.51078 1.590 0.11185
## typeC:regionNE 1.88187 0.39010 4.824 1.41e-06 ***
## typeC:regionMW -1.38434 1.06525 -1.300 0.19376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 392.76 on 79 degrees of freedom
## Residual deviance: 276.70 on 70 degrees of freedom
## AIC: 509.33
##
## Number of Fisher Scoring iterations: 6
Esses resultados fornecem evidências convincentes de uma interação entre o efeito de região e o tipo de instituição. Um teste de queda no desvio confirma a significância da contribuição da interação para este modelo. Temos evidências estatisticamente significativas \(\chi^2 = 71.98\), \(df=4\) e p-valor<0.001 que a diferença entre faculdades e universidades na taxa de crimes violentos difere por região.
Por exemplo, nosso modelo estima que as taxas de crimes violentos são 6.6 (\(\exp(1.88187)\)) vezes mais altas nas universidades do Nordeste (NE) em comparação com as faculdades, enquanto no Nordeste estimamos que as taxas de crimes violentos sejam 2.4 (\(\exp(-1.00806+1.88187)\)) vezes maior nas faculdades.
anova(modeltr, modeli, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: nv ~ type + region
## Model 2: nv ~ type + region + region:type
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 74 348.68
## 2 70 276.70 4 71.981 8.664e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
O desvio residual (276.70 com 70 df) sugere falta de ajuste significativa no modelo de interação (p <0.001). Uma possibilidade é que existam outras covariáveis importantes que poderiam ser usadas para descrever as diferenças nas taxas de crimes violentos. Sem covariáveis adicionais a serem consideradas, procuramos observações extremas, mas já eliminamos a mais extrema das observações.
Na ausência de outras covariáveis ou observações extremas, consideramos a superdispersão como uma possível explicação para a falta de ajuste significativa.
A superdispersão sugere que há mais variação na resposta do que o modelo sugere. Em um modelo de Poisson, esperaríamos que as médias e variâncias da resposta fossem aproximadamente as mesmas em vários grupos. Sem ajustar a superdispersão, usamos erros padrão incorretos e artificialmente pequenos, levando a p-valores artificialmente pequenos para os coeficientes do modelo. Também podemos acabar com modelos artificialmente complexos.
Podemos levar em consideração a superdispersão de várias maneiras diferentes. O mais simples é usar um fator de dispersão estimado para aumentar os erros padrão. Outra maneira é usar um modelo de regressão binomial negativo. Começamos usando uma estimativa do parâmetro de dispersão.
modeliq <- glm(nv ~ type + region + region:type, family = quasipoisson(link = "log"),
offset = log(enroll1000), data = Dados)
summary(modeliq)
##
## Call:
## glm(formula = nv ~ type + region + region:type, family = quasipoisson(link = "log"),
## data = Dados, offset = log(enroll1000))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7292 -1.5701 -0.8148 0.9364 4.9575
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.62229 0.17950 -3.467 0.000905 ***
## typeC -1.00806 0.72545 -1.390 0.169066
## regionW -0.07891 0.34807 -0.227 0.821322
## regionC -0.65587 0.33201 -1.975 0.052160 .
## regionNE -0.17275 0.29995 -0.576 0.566516
## regionMW -0.43593 0.32422 -1.345 0.183102
## typeC:regionW -1.59851 1.68478 -0.949 0.345990
## typeC:regionC 0.81212 1.07708 0.754 0.453380
## typeC:regionNE 1.88187 0.82260 2.288 0.025179 *
## typeC:regionMW -1.38434 2.24627 -0.616 0.539707
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 4.446556)
##
## Null deviance: 392.76 on 79 degrees of freedom
## Residual deviance: 276.70 on 70 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6