Exemplo No.12
Fernando Lucambio Pérez
08/04/2021
Neste exemplo, os pesquisadores montaram parcelas de jardim com diferentes suítes de plantas, com cada suíte identificada como um nível da variável Garden abaixo. Em setembro, eles contaram o número de borboletas monarca Monarchs em cada canteiro de jardim.
Dados = data.frame(Garden=c(rep("A",8),rep("B",8),rep("C",8)), Monarchs=c(0,4,2,2,0,6,0,0,5,9,7,5,7,5,9,5,10,14,12,12,10,16,10,10))Os fatores foram ordenados pela ordem no arquivo de dados, Caso contrário, o R irá colocá-los em ordem alfabética.
Dados$Garden = factor(Dados$Garden, levels = unique(Dados$Garden))Verifique o quadro de dados
library(psych)
headTail(Dados)## Garden Monarchs
## 1 A 0
## 2 A 4
## 3 A 2
## 4 A 2
## ... <NA> ...
## 21 C 10
## 22 C 16
## 23 C 10
## 24 C 10str(Dados)## 'data.frame': 24 obs. of 2 variables:
## $ Garden : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 2 2 ...
## $ Monarchs: num 0 4 2 2 0 6 0 0 5 9 ...summary(Dados)## Garden Monarchs
## A:8 Min. : 0.000
## B:8 1st Qu.: 3.500
## C:8 Median : 6.500
## Mean : 6.667
## 3rd Qu.:10.000
## Max. :16.000Histogramas
Utilizando a opção layout indicamos colunas e linhas de parcelas individuais.
library(lattice)
histogram( ~ Monarchs | Garden, data = Dados, layout=c(3,1))
mean(Dados$Monarchs[Dados$Garden=="A"]); var(Dados$Monarchs[Dados$Garden=="A"])## [1] 1.75## [1] 5.071429mean(Dados$Monarchs[Dados$Garden=="B"]); var(Dados$Monarchs[Dados$Garden=="B"])## [1] 6.5## [1] 3.142857mean(Dados$Monarchs[Dados$Garden=="C"]); var(Dados$Monarchs[Dados$Garden=="C"])## [1] 11.75## [1] 5.071429Exemplo de regressão Poisson
A regressão Poisson faz certas suposições sobre a relação entre a média e a dispersão da variável dependente. Como essa suposição pode não ser atendida para todos os conjuntos de dados, a regressão Poisson pode não ser recomendada para uso rotineiro. Particularmente, a regressão Poisson clássica deve ser evitada se houver superdispersão nos dados ou se houver várias contagens zero na variável dependente.
Uma abordagem alternativa para dados com superdispersão e a regressão Binomial Negativa e uma abordagem alternativa para dados com muitos zeros é a regressão Poisson inflada de zero.
model.p = glm(Monarchs ~ Garden, data = Dados, family="poisson")
summary(model.p)##
## Call:
## glm(formula = Monarchs ~ Garden, family = "poisson", data = Dados)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8708 -0.6135 -0.2257 0.3045 2.5071
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5596 0.2672 2.094 0.0363 *
## GardenB 1.3122 0.3011 4.358 1.31e-05 ***
## GardenC 1.9042 0.2865 6.648 2.98e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 95.120 on 23 degrees of freedom
## Residual deviance: 28.657 on 21 degrees of freedom
## AIC: 110.86
##
## Number of Fisher Scoring iterations: 5library(car)
Anova(model.p, type="II", test="LR")## Analysis of Deviance Table (Type II tests)
##
## Response: Monarchs
## LR Chisq Df Pr(>Chisq)
## Garden 66.463 2 3.697e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(rcompanion)
nagelkerke(model.p)## $Models
##
## Model: "glm, Monarchs ~ Garden, poisson, Dados"
## Null: "glm, Monarchs ~ 1, poisson, Dados"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.387929
## Cox and Snell (ML) 0.937293
## Nagelkerke (Cragg and Uhler) 0.938037
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -2 -33.231 66.463 3.6967e-15
##
## $Number.of.observations
##
## Model: 24
## Null: 24
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"Calculemos o \(R^2_V\):
library(rsq)
rsq.v(model.p)## [1] 0.8114865library(multcompView)
library(emmeans)
marginal = emmeans(model.p, ~ Garden)
pairs(marginal, adjust="sidak")## contrast estimate SE df z.ratio p.value
## A - B -1.312 0.301 Inf -4.358 <.0001
## A - C -1.904 0.286 Inf -6.648 <.0001
## B - C -0.592 0.173 Inf -3.426 0.0018
##
## Results are given on the log (not the response) scale.
## P value adjustment: sidak method for 3 testsmultcomp::cld(marginal, alpha=0.05,
Letters=letters, ### Use letras minúsculas para .group
adjust="sidak") ### Ajuste de "Tukey" para comparações múltiplas; orque "Tukey" é apropriado apenas ### para um conjunto de comparações em pares## Garden emmean SE df asymp.LCL asymp.UCL .group
## A 0.56 0.267 Inf -0.0785 1.20 a
## B 1.87 0.139 Inf 1.5407 2.20 b
## C 2.46 0.103 Inf 2.2176 2.71 c
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## P value adjustment: sidak method for 3 tests
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.### Observe que as estimativas estão em escala logarítmicaResíduos
library(car)
residualPlots(model.p)
influenceIndexPlot(model.p)
library(ggfortify)
autoplot(model.p)
residuos = residuals(model.p, type="pearson")
plot(density(residuos), ylab = "Densidade estimada", main ="", xlim = c(-5,5), ylim = c(0,0.5))
grid()
rug(residuos, lwd = 2)
curve(dnorm(x, mean = mean(residuos), sd = sd(residuos)), add = TRUE, col = "red", lwd = 2)
Exemplo de regressão Binomial Negativa
A regressão Binomial Negativa é semelhante na aplicação à regressão de Poisson, mas permite a superdispersão na variável de contagem dependente. Lembremos que se \(Y\sim BN(n,p)\) então: \(\mbox{E}(Y)=\mu = n(1-p)/p\) e \(\mbox{Var}(Y)=n(1-p)/p^2\).
Este exemplo usará a função glm.nb no pacote MASS. A função Anova no pacote car será usada para uma análise de desvio e a função nagelkerke será usada para determinar um valor p e pseudo R-quadrado para o modelo. A análise post-hoc pode ser conduzida com o pacote emmeans.
library(MASS)
model.nb = glm.nb(Monarchs ~ Garden, data=Dados, control = glm.control(maxit=10000))
model.nb##
## Call: glm.nb(formula = Monarchs ~ Garden, data = Dados, control = glm.control(maxit = 10000),
## init.theta = 510466474900, link = log)
##
## Coefficients:
## (Intercept) GardenB GardenC
## 0.5596 1.3122 1.9042
##
## Degrees of Freedom: 23 Total (i.e. Null); 21 Residual
## Null Deviance: 95.12
## Residual Deviance: 28.66 AIC: 112.8car::Anova(model.nb, type="II", test="LR")## Analysis of Deviance Table (Type II tests)
##
## Response: Monarchs
## LR Chisq Df Pr(>Chisq)
## Garden 66.464 2 3.694e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1rcompanion::nagelkerke(model.nb)## $Models
##
## Model: "glm.nb, Monarchs ~ Garden, Dados, glm.control(maxit = 10000), 510466474900, log"
## Null: "glm.nb, Monarchs ~ 1, Dados, glm.control(maxit = 10000), 1.749060129, log"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.255141
## Cox and Snell (ML) 0.776007
## Nagelkerke (Cragg and Uhler) 0.778217
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -2 -17.954 35.907 1.5952e-08
##
## $Number.of.observations
##
## Model: 24
## Null: 24
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"rsq::rsq.v(model.nb)## [1] 0.8114857marginal = emmeans(model.nb, ~ Garden)
pairs(marginal, adjust="sidak")## contrast estimate SE df z.ratio p.value
## A - B -1.312 0.301 Inf -4.358 <.0001
## A - C -1.904 0.286 Inf -6.647 <.0001
## B - C -0.592 0.173 Inf -3.426 0.0018
##
## Results are given on the log (not the response) scale.
## P value adjustment: sidak method for 3 testsmultcomp::cld(marginal,
alpha = 0.05,
Letters = letters,
type = "response", ### Relatório emmeans na escala original
adjust = "sidak") ## Garden response SE df asymp.LCL asymp.UCL .group
## A 1.75 0.468 Inf 0.924 3.31 a
## B 6.50 0.901 Inf 4.668 9.05 b
## C 11.75 1.212 Inf 9.185 15.03 c
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## Intervals are back-transformed from the log scale
## P value adjustment: sidak method for 3 tests
## Tests are performed on the log scale
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.Vejamos um poucoos resíduos.
residualPlots(model.nb)
influenceIndexPlot(model.nb)
autoplot(model.nb)
residuos = residuals(model.nb, type="pearson")
plot(density(residuos), ylab = "Densidade estimada", main ="", xlim = c(-5,5), ylim = c(0,0.5))
grid()
rug(residuos, lwd = 2)
curve(dnorm(x, mean = mean(residuos), sd = sd(residuos)), add = TRUE, col = "red", lwd = 2)
Exemplo de regressão inflacionada com zeros
A regressão inflada de zeros é semelhante na aplicação à regressão Poisson, mas permite uma abundância de zeros na variável de contagem dependente.
Este exemplo usará a função zeroinfl no pacote pscl. A função Anova no pacote car será usada para uma análise de desvio assim como outras funções.
library(pscl)
model.zi = zeroinfl(Monarchs ~ Garden, data = Dados, dist = "poisson") ### dist = "negbin" pode ser usado
summary(model.zi)##
## Call:
## zeroinfl(formula = Monarchs ~ Garden, data = Dados, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.8156 -0.5883 -0.2188 0.3112 1.9807
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2182 0.2847 4.278 1.89e-05 ***
## GardenB 0.6536 0.3167 2.064 0.039 *
## GardenC 1.2457 0.3029 4.113 3.90e-05 ***
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.046e-02 7.363e-01 -0.096 0.924
## GardenB -2.057e+01 1.071e+04 -0.002 0.998
## GardenC -2.057e+01 1.071e+04 -0.002 0.998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 9
## Log-likelihood: -48.14 on 6 DfAnova(model.zi, type="II", test="Chisq")## Analysis of Deviance Table (Type II tests)
##
## Response: Monarchs
## Df Chisq Pr(>Chisq)
## Garden 2 23.914 6.414e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1rcompanion::nagelkerke(model.zi)## $Models
##
## Model: "zeroinfl, Monarchs ~ Garden, Dados, poisson"
## Null: "zeroinfl, Monarchs ~ 1, Dados, poisson"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.284636
## Cox and Snell (ML) 0.797356
## Nagelkerke (Cragg and Uhler) 0.800291
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -4 -19.156 38.311 9.6649e-08
##
## $Number.of.observations
##
## Model: 24
## Null: 24
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"marginal = emmeans(model.zi, ~ Garden)
pairs(marginal, adjust="sidak")## contrast estimate SE df z.ratio p.value
## A - B -4.75 1.18 Inf -4.032 0.0002
## A - C -10.00 1.43 Inf -6.994 <.0001
## B - C -5.25 1.51 Inf -3.476 0.0015
##
## P value adjustment: sidak method for 3 testsmultcomp::cld(marginal, alpha=0.05, Letters=letters, adjust="sidak") ## Garden emmean SE df asymp.LCL asymp.UCL .group
## A 1.75 0.759 Inf -0.0614 3.56 a
## B 6.50 0.901 Inf 4.3477 8.65 b
## C 11.75 1.212 Inf 8.8563 14.64 c
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## P value adjustment: sidak method for 3 tests
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.residuos = residuals(model.zi, type="pearson")
plot(density(residuos), ylab = "Densidade estimada", main ="", xlim = c(-5,5), ylim = c(0,0.6))
grid()
rug(residuos, lwd = 2)
curve(dnorm(x, mean = mean(residuos), sd = sd(residuos)), add = TRUE, col = "red", lwd = 2)
A regressão Poisson robusta
Todos os métodos estatísticos dependem explícita ou implicitamente de uma série de suposições. Essas suposições geralmente visam formalizar o que o estatístico sabe ou conjectura sobre o problema de análise de dados ou modelagem estatística que ele enfrenta e, ao mesmo tempo, visam tornar o modelo resultante gerenciável do ponto de vista teórico e computacional. No entanto, é geralmente entendido que os modelos formais resultantes são simplificações da realidade e que sua validade é, na melhor das hipóteses, aproximada.
A abordagem robusta de modelagem estatística e análise de dados visa derivar métodos que produzam estimativas de parâmetros confiáveis e testes associados e intervalos de confiança, não apenas quando os dados seguem exatamente uma determinada distribuição, mas também quando isso acontece apenas aproximadamente no sentido descrito. Uma caracterização mais informal orientada a dados de métodos robustos é que eles se ajustam bem à maior parte dos dados: se os dados não contêm outliers, o método robusto fornece aproximadamente os mesmos resultados que o método clássico, enquanto se uma pequena proporção de outliers estiver presente, o método robusto fornece aproximadamente os mesmos resultados que o método clássico. O método robusto fornece aproximadamente os mesmos resultados que o método clássico aplicado aos dados típicos, ou seja, aos dados sem os outliers. Como consequência de ajustar bem a maior parte dos dados, métodos robustos fornecem um método muito confiável de detectar outliers, mesmo em situações multivariadas de alta dimensão.
Estatísticas robustas lidam com desvios de modelos paramétricos ideais e seus perigos para os procedimentos estatísticos derivados sob o modelo assumido. Seu objetivo principal é o desenvolvimento de procedimentos que ainda são confiáveis e razoavelmente eficientes sob pequenos desvios do modelo, isto é, quando a distribuição subjacente reside em uma vizinhança do modelo assumido. Estatísticas robustas é então uma extensão de estatísticas paramétricas, levando em conta que os modelos paramétricos são, na melhor das hipóteses, apenas aproximações para a realidade. Pode-se considerar Tukey (1960), Huber (1964), e Hampel (1968) os documentos fundamentais que estabeleceram as bases de estatísticas robustas modernas. Outros livros amplamente expositivos desta teoría são Hampel, Ronchetti, Rousseeuw and Stahel (1986) e Maronna, Martin and Yohai (2006).
A estatística robusta produzir estimadores consistentes e possivelmente eficientes e estatísticas de teste com nível estável quando o modelo é levemente especificado. O erro de especificação do modelo abrange um conjunto relativamente grande de possibilidades e estatísticas robustas não podem lidar com todos os tipos de erro de especificação do modelo.
Por um pequeno erro de especificação do modelo, supomos que o processo de geração de dados está em uma vizinhança do modelo verdadeiro ou modelo postulado; aquele que é considerado útil para o problema sob investigação. Por uma vizinhança do modelo postulado consideramos \[\begin{equation*} F_\epsilon = (1 - \epsilon)F_\theta + \epsilon G, \end{equation*}\] sendo que \(F_\theta\) é o modelo postulado, \(\theta\) é um conjunto de parâmetros de interesse, G é uma distribuição arbitrária e \(0\leq \epsilon\leq 1\) captura a quantidade de especificação incorreta do modelo.
Aplicado aos Modelos Lineares Generalizados, sabemos que \[\begin{equation*}
g(\mu_i)=x_i^\top \beta,
\end{equation*}\] onde \(\mbox{E}(Y_i)=\mu_i\), \(\mbox{Var}(Y_i)=V(\mu_i)\) e \(r_i=(y_i-\mu_i)/\sqrt{\phi V(\mu_i)}\), o estimador robusto é definido como solução so sistema de equações \[\begin{equation*}
\sum_{i=1}^n \left( \psi_c(r_i)w(x_i)\dfrac{\mu_i'}{\sqrt{\phi V(\mu_i)}}-a(\beta)\right) = 0,
\end{equation*}\] onde \[\begin{equation*}
\mu_i'= \dfrac{\partial \mu_i}{\partial\beta}=\dfrac{\partial \mu_i}{\partial\eta_i} x_i
\end{equation*}\]
e \[\begin{equation*}
a(\beta)= \dfrac{1}{n}\sum_{i=1}^n \mbox{E}\big(\psi_c(r_i)\big)w(x_i)\dfrac{\mu_i'}{\sqrt{\phi V(\mu_i)}}\cdot \end{equation*}\]
O termo \(a(\beta)\) é uma correção para garantir a consistência.
A regressão robusta Poisson é robusta para outliers na variável dependente. Este exemplo usa a função glmRob no pacote robust. A função anova pode ser usada para conduzir uma análise de desvio. O p-valor para o modelo pode ser encontrado comparando o modelo a um modelo nulo.
No momento em que escrevemos estas notas, a função glmRob só pode usar as famílias de modelos Poisson e binomial. Um método alternativo é a função glmrob no pacote robustbase.
library(robust)
model.rob = glmRob(Monarchs ~ Garden, data = Dados, family = "poisson")
summary(model.rob)##
## Call: glmRob(formula = Monarchs ~ Garden, family = "poisson", data = Dados)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6930 -0.6223 -0.1019 0.5301 2.8370
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3599 0.3007 1.197 2.313e-01
## GardenB 1.5155 0.3352 4.522 6.138e-06
## GardenC 2.0630 0.3239 6.370 1.895e-10
##
## (Dispersion Parameter for poisson family taken to be 1 )
##
## Null Deviance: 430.2 on 23 degrees of freedom
##
## Residual Deviance: 29.33699 on 21 degrees of freedom
##
## Number of Iterations: 4
##
## Correlation of Coefficients:
## (Intercept) GardenB
## GardenB -0.8970
## GardenC -0.9283 0.8327anova(model.rob, test="Chisq")## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL NA NA 23 430.19850 NA
## Garden 2 400.8615 21 29.33699 4.906414e-63model.rob.null = glmRob(Monarchs ~ 1, data = Dados, family = "poisson")
anova(model.rob.null, model.rob, test="Chisq")## Terms Resid. Df Resid. Dev Test Df Deviance Pr(>Chi)
## 1 1 23 95.12606 NA NA NA
## 2 Garden 21 29.33699 2 65.78907 5.94098e-11Regressão quase-Poisson
A regressão quase-Poisson é útil porque tem um parâmetro de dispersão variável, de modo que pode modelar dados superdispersos. Pode ser melhor do que a regressão binomial negativa em algumas circunstâncias (Verhoef and Boveng, 2007).
No momento que escrevemos estas notas, a regressão Quasi-Poisson não tem um conjunto completo de funções de suporte em R. Usando a opção de família quasipoisson na função glm, os resultados terão os mesmos coeficientes de parâmetros da opção poisson, mas as estatísticas de inferência são ajustados na função de resumo.
model.qp = glm(Monarchs ~ Garden, data=Dados, family="quasipoisson")
summary(model.qp)##
## Call:
## glm(formula = Monarchs ~ Garden, family = "quasipoisson", data = Dados)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8708 -0.6135 -0.2257 0.3045 2.5071
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5596 0.3013 1.857 0.077349 .
## GardenB 1.3122 0.3395 3.866 0.000896 ***
## GardenC 1.9042 0.3230 5.896 7.48e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.271146)
##
## Null deviance: 95.120 on 23 degrees of freedom
## Residual deviance: 28.657 on 21 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Anova(model.qp, type="II", test="LR")## Analysis of Deviance Table (Type II tests)
##
## Response: Monarchs
## LR Chisq Df Pr(>Chisq)
## Garden 52.286 2 4.429e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1marginal = emmeans(model.qp, ~ Garden)
pairs(marginal, adjust="tukey")## contrast estimate SE df z.ratio p.value
## A - B -1.312 0.339 Inf -3.866 0.0003
## A - C -1.904 0.323 Inf -5.896 <.0001
## B - C -0.592 0.195 Inf -3.038 0.0067
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatesmultcomp::cld(marginal, alpha=0.05, Letters=letters, adjust="tukey")## Garden emmean SE df asymp.LCL asymp.UCL .group
## A 0.56 0.301 Inf -0.16 1.28 a
## B 1.87 0.156 Inf 1.50 2.25 b
## C 2.46 0.116 Inf 2.19 2.74 c
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.rsq::rsq.v(model.qp)## [1] 0.8114865Regressão de Hermite
A distribuição Hermite generalizada é uma distribuição mais geral que pode lidar com a superdispersão ou multimodalidade (Moriña and others, 2015). Isso torna a regressão Hermite generalizada uma ferramenta poderosa e flexível para modelar dados de contagem. É implementado com o pacote
Ajustar modelos com o pacote hermite pode ser um pouco difícil. Um problema é que o ajuste do modelo pode falhar sem que alguns parâmetros sejam especificados. Freqüentemente, especificar um valor apropriado para a opção m ajudará.
Uma outra dificuldade com esta abordagem é que, no momento da escrita, o pacote não é suportado pela função anova para comparar modelos, a função Anova para testar efeitos ou outras funções úteis como emmeans para efeitos de fator.
-
Moriña, D., M. Higueras, P. Puig, and M. Oliveira. 2015. Generalized Hermite Distribution Modelling with the R Package hermite.
The R Journal 7(2):263–274. journal.r-project.org/archive/2015-2/morina-higueras-puig-etal.pdf. - Wilson, P. 2015. The Misuse of the Vuong Test for Non-Nested Models to Test for Zero-Inflation. Economic Letters 127: 51–53. cybermetrics.wlv.ac.uk/paperdata/misusevuong.pdf.
library(hermite)
model.he = glm.hermite(Monarchs ~ Garden, data = Dados, link = "log", m=3)
summary(model.he)## Call:
## glm.hermite(formula = Monarchs ~ Garden, data = Dados, link = "log",
## m = 3)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -11.7955649 -10.7955649 -5.5414900 -1.6621439 -0.6621439
##
## Coefficients:
## Estimate Std. Error z value p-value
## (Intercept) 0.5081083 0.3251349 1.5627612 1.181088e-01
## GardenB 1.3700567 0.3641379 3.7624662 1.682461e-04
## GardenC 1.9596153 0.3476326 5.6370291 1.730089e-08
## dispersion.index 1.0820807 0.2877977 0.1281707 3.601681e-01
## order 3.0000000 NA NA NA
## (Likelihood ratio test against Poisson is reported by *z value* for *dispersion.index*)
##
## AIC: 112.7762Medianas e intervalos de confiança
library(rcompanion)
Sum = groupwiseMedian(Monarchs ~ Garden,
data=Dados,
conf=0.95,
R=5000,
percentile=TRUE,
bca=FALSE,
digits=3)
Sum## Garden n Median Conf.level Percentile.lower Percentile.upper
## 1 A 8 1 0.95 0 4
## 2 B 8 6 0.95 5 9
## 3 C 8 11 0.95 10 14Gráfico de medianas e intervalos de confiança
O quadro de dados Sum criado acima será passado para o ggplot para plotagem. No final do código, annotate é usado para adicionar texto ao gráfico para indicar quais medianas são significativamente diferentes umas das outras.
library(ggplot2)
ggplot(Sum, ### O arquivo de dados a ser utilizado.
aes(x = Garden,
y = Median)) +
geom_errorbar(aes(ymin = Percentile.lower,
ymax = Percentile.upper),
width = 0.05,
size = 1) +
geom_point(shape = 15,
size = 5) +
theme_bw() +
theme(axis.title = element_text(face = "bold")) +
ylab("Contagem mediana de borboletas monarca") +
annotate("text",
x = 1:3,
y = c(5, 10, 15),
label = c("Grupo 3", "Grupo 2", "Grupo 1"))
Código opcional para teste de ajuste perfeito do qui-quadrado
Uma abordagem alternativa para lidar com dados de contagem é somar as contagens para tratamentos e usar um teste de qui-quadrado ou teste relacionado. Aqui, um teste de qualidade de ajuste qui-quadrado é usado para ver se as contagens diferem das proporções iguais “esperadas”.
Tabela = xtabs(Monarchs ~ Garden,
data = Dados)
Tabela## Garden
## A B C
## 14 52 94chisq.test(Tabela)##
## Chi-squared test for given probabilities
##
## data: Tabela
## X-squared = 60.05, df = 2, p-value = 9.127e-14teste qui-quadrado.
Garden.A = sum(Dados$Monarchs[Dados$Garden=="A"])
Garden.B = sum(Dados$Monarchs[Dados$Garden=="B"])
Garden.C = sum(Dados$Monarchs[Dados$Garden=="C"])
observed = c(Garden.A, Garden.B) # frequências observadas
expected = c(1/2, 1/2) # proporções esperadas
chisq.test(x = observed, p = expected)##
## Chi-squared test for given probabilities
##
## data: observed
## X-squared = 21.879, df = 1, p-value = 2.904e-06observed = c(Garden.B, Garden.C) # frequências observadas
expected = c(1/2, 1/2) # proporções esperadas
chisq.test(x = observed, p = expected)##
## Chi-squared test for given probabilities
##
## data: observed
## X-squared = 12.082, df = 1, p-value = 0.0005091observed = c(Garden.A, Garden.C) # frequências observadas
expected = c(1/2, 1/2) # proporções esperadas
chisq.test(x = observed, p = expected)##
## Chi-squared test for given probabilities
##
## data: observed
## X-squared = 59.259, df = 1, p-value = 1.382e-14Teste de Vuong para comparar modelos Poisson, Binomial Negativo e inflacionado de zeros
O teste Vuong, implementado pelo pacote pscl, pode testar dois modelos não aninhados. Ele funciona com negbin, zeroinfl e alguns objetos de modelo glm que são ajustados aos mesmos dados.
A hipótese nula é que não há diferença nos modelos. A função produz três testes, um teste “bruto”, um corrigido pelo AIC e um corrigido pelo BIC, qualquer um dos quais pode ser usado.
Foi sugerido que o teste Vuong não seja usado para testar a inflação zeros (Wilson, 2015).
# Teste Vuong
library(pscl)
vuong(model.p,
model.nb,
digits = 4)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw 0.03324988 model1 > model2 0.48674
## AIC-corrected 0.03324988 model1 > model2 0.48674
## BIC-corrected 0.03324988 model1 > model2 0.48674### A estatística z de Vuong positiva sugere que o modelo 1 é superior,
### mas, neste caso, a diferença não é significativa,
### e o valor da estatística é provavelmente muito pequeno para ser significativo.
vuong(model.p,
model.zi,
digits = 4)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -1.4424725 model2 > model1 0.074585
## AIC-corrected -0.4335210 model2 > model1 0.332318
## BIC-corrected 0.1607786 model1 > model2 0.436134### A estatística z de Vuong negativa sugere que o modelo 2 é superior.
### Se a estatística bruta for usada, p = 0.07 fornece alguma evidência
### que o modelo zi é superior.
vuong(model.nb,
model.zi,
digits = 4)## Vuong Non-Nested Hypothesis Test-Statistic:
## (test-statistic is asymptotically distributed N(0,1) under the
## null that the models are indistinguishible)
## -------------------------------------------------------------
## Vuong z-statistic H_A p-value
## Raw -1.4424725 model2 > model1 0.074585
## AIC-corrected -0.4335210 model2 > model1 0.332318
## BIC-corrected 0.1607786 model1 > model2 0.436134### A estatística z de Vuong negativa sugere que o modelo 2 é superior.
### Se a estatística bruta for usada, p = 0.07 fornece alguma evidência
### que o modelo zi é superior.