MAN5728: Tópicos Essenciais da Bioestatística
Módulo 6: GEE (Generalized Estimating Equations) e GLMM (Generalized Linear Mixed Model)
07 julho 2025 18:47h
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))options(warn=-1)
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(HSAUR3, warn.conflicts=FALSE))
suppressMessages(library(gee, warn.conflicts=FALSE))
suppressMessages(library(lme4, warn.conflicts=FALSE))
suppressMessages(library(lmerTest, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(rmcorr, warn.conflicts=FALSE))
suppressMessages(library(weights, warn.conflicts=FALSE))
suppressMessages(library(lsr, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(rstatix, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(magrittr, warn.conflicts=FALSE))
suppressMessages(library(purrr, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))
suppressMessages(library(rsample, warn.conflicts=FALSE))
suppressMessages(library(ggpubr, warn.conflicts=FALSE))
suppressMessages(library(IDPmisc, warn.conflicts=FALSE))
suppressMessages(library(estimatr, warn.conflicts=FALSE))
suppressMessages(library(broom.mixed, warn.conflicts=FALSE))
suppressMessages(library(geepack, warn.conflicts=FALSE))Material
- HTML de R Markdown em
RPubs - Prof. José Siqueira: ResearchGate
Conteúdo
- GEE (Generalized Estimating Equations)
- GLMM (Generalized Linear Mixed Model)
Análise de dados agrupados e longitudinais
Estudos longitudinais são caracterizados por mensurações realizadas nas mesmas unidades experimentais que na Medicina consistem, por exemplo, em pacientes, profissionais da saúde, serviços médicos e unidades familiares, ao longo do tempo, permitindo o estudo direto de mudanças decorrentes da evolução de uma doença ou do efeito sequencial de uma intervenção.
O objetivo principal de um estudo longitudinal é caracterizar os fatores que influenciam as mudanças observadas de maneira mais precisa do que se pode obter com uma única observação.
Uma outra característica dos dados é que eles podem ser agrupados hierarquicamente. Em estudos com agrupamentos hierarquizados, os agrupamentos são compostos por medidas repetidas de uma unidade experimental.
As observações dentro de um nível hierárquico de agrupamento exibem frequentemente correlações positivas que devem ser levadas em consideração na análise, pois esperamos que medidas em unidades experimentais de um agrupamento sejam mais similares que as medidas em unidades de diferentes agrupamentos.
Famílias, serviços hospitalares, práticas médicas, bairros e escolas, por exemplo, são instâncias nas quais naturalmente ocorre o efeito de agrupamento.
Um estudo pode ser longitudinal e com agrupamentos hierarquizados, o que exige conhecimento do pesquisador para, corretamente, obter inferências sobre os fatores que estuda.
No decorrer da década de 80 ocorreram importantes inovações computacionais, principalmente em termos de aprimoramentos de algoritmos iterativos, que propiciaram o desenvolvimento de novas classes de modelos, mais complexos e gerais.
Dentre esses modelos, dois se destacam:
- GEE (Generalized Estimating Equations)
- GLMM (Generalized Linear Mixed Models)
Ambos os modelos deram consideráveis saltos de qualidade para a modelagem de dados relacionados aos problemas da área da saúde, em particular, e da ciência, em geral.
Os dois modelos surgiram em departamentos de Bioestatística:
- GEE no Departamento de Bioestatística da Johns Hopkins University em 1986;
- GLMM no Departamento de Bioestatística da Harvard School of Public Health em 1982.
Os sistemas computacionais, no entanto, permitem que profissionais da saúde tirem proveito de tais técnicas estatísticas, desde que as entendam conceitualmente e possam escolher adequadamente o que utilizar e como interpretar seus resultados.
Figure 1: Comparison of the numbers of papers that used three
statistical methods published in five journals in the behavioural
sciences between 2001 and 2014.
Pekár & Brabec, 2017, p. 88.
“Os métodos estatísticos tradicionais são baseados na suposição de que as observações em uma amostra são independentes umas das outras, ou seja, o valor de uma observação não é influenciado pelo valor de outra.
Essa suposição de independência é violada se os dados são agrupados (clustered), ou seja, se as observações em um agrupamento tenderem a ser mais semelhantes entre si do que com os indivíduos do restante da amostra.
Os dados agrupados surgem de três maneiras principais:
- Medidas repetidas em estudos longitudinais: Neste caso os clusters são os sujeitos; observações repetidas sobre o mesmo sujeito serão mais semelhantes entre si do que observações sobre outros sujeitos.
Por exemplo:
em estudos de asma ou outras doenças crônicas, os episódios da doença podem ocorrer em mais de uma ocasião no mesmo sujeito;
em estudos longitudinais de doenças comuns da infância em países em desenvolvimento, as crianças podem apresentar vários episódios de diarreia, malária ou infecções respiratórias agudas durante o estudo;
em um estudo de doenças cardiovasculares e obesidade, as medições da pressão arterial, índice de massa corporal e níveis de colesterol podem ser repetidas a cada 3 meses.
Várias medidas sobre o mesmo sujeito: Por exemplo, em pesquisas odontológicas são feitas observações em mais de um dente no mesmo sujeito. Neste caso, os agrupamentos (clusters) são novamente os sujeitos.
Estudos em que os sujeitos são agrupados: Isso ocorre, por exemplo, em:
ensaios randomizados em cluster, nos quais grupos, em vez de indivíduos, são randomizados para receber as diferentes intervenções em estudo. Por exemplo, a unidade de randomização pode ser prática geral (https://en.wikipedia.org/wiki/General_practice), com todos os pacientes registrados em uma prática recebendo a mesma intervenção. Como os pacientes em uma prática geral podem ser mais semelhantes entre si do que com pacientes em outras práticas gerais, por exemplo, porque algumas áreas tendem a ser mais carentes do que outras ou devido à exposição a um risco ambiental comum, os dados são agrupados. Neste caso o cluster é o grupo de pacientes cadastrados na prática geral;
estudos de família, uma vez que indivíduos de uma mesma família tendem a ser mais semelhantes entre si do que indivíduos de famílias diferentes, pois compartilham genes e ambiente semelhantes. Neste caso o cluster é a família;
pesquisas em que a amostragem por conglomerados é empregada. Por exemplo, para estimar a porcentagem de jovens de 14 anos em Londres que trabalham nos fins de semana, podemos selecionar 1.000 crianças por amostragem aleatória de 20 escolas de todas as escolas de Londres e, em seguida, amostrar aleatoriamente 50 crianças de cada uma das escolas selecionadas. Como as crianças de uma escola podem ser mais parecidas entre si do que as crianças de escolas diferentes, os dados são agrupados. Neste caso, os clusters são as escolas.
É essencial que a presença de agrupamento seja permitida nas análises estatísticas.
A principal razão para isso, como veremos, é que os erros padrão podem ser muito pequenos se não levarem em conta o agrupamento nos dados.
Isso levará a intervalos de confiança muito estreitos e valores P muito pequenos.
Discutiremos quatro maneiras apropriadas de analisar dados agrupados:
- calcular medidas de resumo para cada cluster e analisar essas medidas de resumo usando métodos-padrão;
- usar erros-padrão robustos para corrigir erros-padrão para o agrupamento;
- utilizar modelos de efeitos aleatórios (GLMM) que modelam explicitamente a similaridade entre indivíduos de um mesmo cluster;
- usar equações de estimação generalizadas (GEE) que ajustam tanto os erros-padrão quanto as estimativas de parâmetros para permitir o agrupamento.”
Kirkwood & Sterne, 2003, p. 355-6
“RESUMO DAS ABORDAGENS PARA A ANÁLISE DE DADOS AGRUPADOS
Se os dados estiverem agrupados, é essencial que o agrupamento seja permitido nas análises. Em particular, a falha em permitir o agrupamento pode significar que os erros-padrão das estimativas dos parâmetros são muito pequenos, de modo que os intervalos de confiança são muito estreitos e os valores p são muito pequenos.
É sempre válido derivar medidas-resumo para cada agrupamento e depois analisá-las usando métodos-padrão.
O provável efeito do agrupamento nos erros-padrão pode ser avaliado pela especificação de erros-padrão robustos que permitem o agrupamento. As estimativas de parâmetros não são afetadas. Para que esses erros-padrão robustos sejam confiáveis, precisamos de um número razoável de agrupamentos (pelo menos 30). Testes de Wald, em vez de testes de razão de verossimilhança (LR), devem ser usados.
Os modelos de efeitos aleatórios (multinível), GLMM, permitem a presença de agrupamento modificando o preditor linear por uma quantidade constante no agrupamento. Supõe-se que os efeitos aleatórios variem aleatoriamente entre os agrupamentos. Modelos de efeitos aleatórios funcionam bem para desfechos normalmente distribuídos e regressão de Poisson, mas a estimativa de modelos logísticos de efeitos aleatórios é difícil e computacionalmente exigente.
Equações de estimação generalizadas (GEE) modificam tanto as estimativas de parâmetros quanto os erros-padrão para permitir o agrupamento. Novamente, deve haver um número razoável de agrupamentos A abordagem GEE é particularmente útil em análises de regressão logística e quando o foco de interesse está no efeito de exposição estimado e onde o agrupamento não tem interesse intrínseco.”
Kirkwood & Sterne, 2003, p. 369-70
Análise de medidas repetidas
Hamaker, 2013, in Mehl & Conner, p. 43-61.
Sainani, 2010, p. 861.
Formato de arquivo de dados
Largo (wide)
- Cada linha do arquivo corresponde a uma unidade experimental distinta
Longo (long)
- Cada linha do arquivo corresponde a cada observação da VD de cada unidade experimental distinta.
O formato do arquivo de dados para análise estatística por meio do Modelo Linear Geral (GLM) é o largo.
O formato longo é usado para análises com Modelo Linear Misto Generalizado (GLMM) e Equações de Estimação Generalizadas (GEE).
Análise de correlação com dados agrupados
O coeficiente de correlação de Pearson avalia o grau de linearidade entre duas variáveis intervalares num delineamento entre participantes. A correlação de Pearson varia de -1 a 1, não depende do tamanho da amostra e é adimensional: Wikipedia: Pearson correlation coefficient
Correlação intra: pHi & PaCO2
Disóxia tecidual e a hipóxia são os maiores fatores determinantes do surgimento e propagação da falência de múltiplos órgãos em pacientes críticos.
- A disóxia tecidual ocorre quando há uma quantidade insuficiente de oxigênio disponível nos tecidos do corpo, apesar de haver uma quantidade adequada de oxigênio no sangue arterial. Esse fenômeno difere da hipóxia, que se refere à deficiência de oxigênio no sangue. A disóxia pode ser causada por problemas na utilização do oxigênio pelas células, devido a disfunções metabólicas ou envenenamento por cianeto, que impede que as células utilizem o oxigênio presente.
A medida do pH intramural (pHi) ou pH intramucoso gástrico (pHim) reflete o balanço ácido-básico tecidual em uma das primeiras regiões do corpo a desenvolver disóxia no choque.
O objetivo é avaliar a associação entre pHi e PaCO2 (pressão parcial de CO2 no sangue arterial) (https://en.wikipedia.org/wiki/PCO2).
“Vários estudos clínicos demonstraram que a acidose intramucosa gástrica detectada por esse procedimento prevê um aumento da mortalidade de adultos gravemente doentes em unidades de terapia intensiva (UTI) médicas e cirúrgicas, e que é um melhor preditor de mortalidade por doença crítica do que outras medidas de entrega global de oxigênio e hemodinâmica sistêmica.”
Baigorri et al., 1997
A tonometria gástrica determina o PaCO2 intraluminal, que se presume estar em equilíbrio com o PaCO2 na mucosa gástrica. O pHi pode ser calculado pela equação de Henderson-Hasselbalch, usando o valor de PaCO2 determinado pela tonometria gástrica e a concentração de íon bicarbonato no sangue arterial.
Baigorri et al., 1997
Os dados usados para testar esta associação estão em Bland & Altman (1995a).
Análise descritiva de correlação intraparticipantes
Subject pH PaCO2
1 1 6.68 3.97
2 1 6.53 4.12
3 1 6.43 4.09
4 1 6.33 3.97
5 2 6.85 5.27
6 2 7.06 5.37
7 2 7.13 5.41
8 2 7.17 5.44
9 3 7.40 5.67
10 3 7.42 3.64
11 3 7.41 4.32
12 3 7.37 4.73
13 3 7.34 4.96
14 3 7.35 5.04
15 3 7.28 5.22
16 3 7.30 4.82
17 3 7.34 5.07
18 4 7.36 5.67
19 4 7.33 5.10
20 4 7.29 5.53
21 4 7.30 4.75
22 4 7.35 5.51
23 5 7.35 4.28
24 5 7.30 4.44
25 5 7.30 4.32
26 5 7.37 3.23
27 5 7.27 4.46
28 5 7.28 4.72
29 5 7.32 4.75
30 5 7.32 4.99
31 6 7.38 4.78
32 6 7.30 4.73
33 6 7.29 5.12
34 6 7.33 4.93
35 6 7.31 5.03
36 6 7.33 4.93
37 7 6.86 6.85
38 7 6.94 6.44
39 7 6.92 6.52
40 8 7.19 5.28
41 8 7.29 4.56
42 8 7.21 4.34
43 8 7.25 4.32
44 8 7.20 4.41
45 8 7.19 3.69
46 8 6.77 6.09
47 8 6.82 5.58car::scatterplot(Dados$PaCO2,
Dados$pH,
groups=Dados$Subject,
regLine=FALSE, smooth=FALSE,
ellipse=list(levels=c(.68),
robust=TRUE,
fill=FALSE),
col="black", xlab="PaCO2", ylab="pHi")
Dados %>%
ggplot2::ggplot(ggplot2::aes(x=PaCO2, y=pH,
color=Subject)) +
ggplot2::geom_point() +
ggplot2::geom_smooth(method="lm", formula=y~x, se=FALSE) +
ggpubr::stat_cor(method="pearson", cor.coef.name="r",
na.rm=TRUE, r.digits=3, p.digits=2,
label.x=6.5) +
ggplot2::labs(x="PaCO2", y="pHi") +
ggplot2::xlim(3,9) +
ggplot2::theme_test()
Correlação intraparticipantes por cor.test:
incorreta
Pearson's product-moment correlation
data: Dados$PaCO2 and Dados$pH
t = -0.43843, df = 45, p-value = 0.6632
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3459063 0.2261851
sample estimates:
cor
-0.06521774 Para um delineamento intraparticipantes, é possível estimar e testar a o coeficiente de correlação de Pearson implícito entre duas variáveis intervalares. Esta correlação implícita intraparticipantes nao é: (i) a média aritmética das correlações dos participantes; (ii) a média ponderada das correlações dos participantes. Bland & Altman (1995a) mostraram a maneira correta de estimar a correlação intra.
Bland & Altman, 1995a
IC95% por bootstrapping das correlações intraparticipantes
boot_corr <- Dados %>%
tidyr::nest(data=-c(Subject)) %>% # grouping the Subject
dplyr::mutate(boots=purrr::map(data,
~rsample::bootstraps(.x, times=1e2,
apparent=FALSE))) %>%
tidyr::unnest(boots) %>% dplyr::mutate(correlations=purrr::map(splits, ~rstatix::cor_test(PaCO2, pH, data=rsample::analysis((.)))))
corr <- boot_corr %>%
tidyr::unnest(correlations) %>% # unnesting tidied data frames
dplyr::select(-data, -splits, -id)
corr <- IDPmisc::NaRV.omit(corr)
# bootstrap confidence interval
point <- Dados %>%
dplyr::group_by(Subject) %>%
rstatix::cor_test(PaCO2, pH, method="pearson")
CI <- corr %>%
dplyr::group_by(Subject) %>%
dplyr::summarise(lwr_CI=quantile(cor, alfa/2),
.estimate=median(cor),
upr_CI=quantile(cor, 1-alfa/2))
print(CI)# A tibble: 7 × 4
Subject lwr_CI .estimate upr_CI
<fct> <dbl> <dbl> <dbl>
1 1 -0.88 -0.053 0.96
2 2 1 1 1
3 3 -0.935 -0.435 0.690
4 4 -0.44 0.49 1
5 5 -0.9 -0.66 0.489
6 6 -1 -0.47 0.625
7 8 -0.98 -0.82 0.395Correlação intraparticipantes por bootstrapping: correta
A estimativa do valor absoluto da correlação de Pearson para medidas repetidas é a raiz quadrada do eta parcial ao quadrado de PaCO2.
O sinal da correlação é obtido por meio do sinal da estimativa do parâmetro beta da regressão entre PaCO2 e pH intramural controlada pelo fator Subject.
# Bland & Altman, 1995a
print(rmcorr::rmcorr(participant=Subject,
measure1=PaCO2,
measure2=pH,
dataset=Dados,
CI.level=1-alfa,
CIs="bootstrap",
nreps=1e4), digits=4)
Repeated measures correlation
r
-0.5067697
degrees of freedom
38
p-value
0.0008471081
95% confidence interval
-0.7262505 -0.06883244 Medida média: “entre participantes”
Dados.media <- aggregate(Dados, x=cbind(PaCO2, pH)~Subject, FUN=mean)
Dados.n <- aggregate(Dados, x=cbind(Subject)~Subject, FUN=length)
Dados.media <- cbind(Dados.media,Dados.n[,2])
names(Dados.media) <- c(names(Dados.media)[1:(length(names(Dados.media))-1)],
"Number")
names(Dados.media) <- c("Subject", "PaCO2.mean", "pH.mean", "Number")
print(Dados.media, digits=4) Subject PaCO2.mean pH.mean Number
1 1 4.037 6.492 4
2 2 5.373 7.053 4
3 3 4.830 7.357 9
4 4 5.312 7.326 5
5 5 4.399 7.314 8
6 6 4.920 7.323 6
7 7 6.603 6.907 3
8 8 4.784 7.115 8Correlação “entre participantes” por cor.test:
incorreta
Pearson's product-moment correlation
data: Dados.media$PaCO2.mean and Dados.media$pH.mean
t = 0.20098, df = 6, p-value = 0.8474
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.6609881 0.7435975
sample estimates:
cor
0.08177314 Correlacao “entre participantes” por weights::wtd.cor:
correta
# Bland & Altman 1995b: analitica e boostrapping
print(weights::wtd.cor(Dados.media$PaCO2.mean,
Dados.media$pH.mean,
weight=Dados.media$Number,
bootp=TRUE, bootse=TRUE, bootn=1e4), digits=4) correlation bootcor std.err t.value p.value
Y 0.06976 -0.005953 0.5862 -0.01016 0.8982Desenhando as correlações intra e entre participantes
Within correlation
r = -0.51
degrees of freedom = 38
p-value = 0.00085
95% confidence interval (analytic) = -0.71 -0.22
95% confidence interval (bootstrap) = -0.73 -0.07
Between correlation
correlation bootcor std.err t.value p.value
0.067 0.004 0.584 0.0066 0.900
Hamaker, 2013, in Mehl & Conner, p. 43-61.
rcGLM: robust clustered General Linear Model
“Dada esta revisão da literatura e nossas novas simulações, achamos que a maioria dos pesquisadores pode continuar usando confortavelmente o software OLS padrão, de preferência com os erros-padrão robustos, em vez de mudar para os métodos semiparamétricos mais recentes.”
Judkins DR & Porter KE, 2016
“estimatr is a package in R dedicated to providing fast
estimators that take into consideration designs often used by social
scientists.
Estimators are statistical methods for estimating quantities of interest like treatment effects or regression parameters.
Many of the estimators included with the R programming language or popular R packages are slow and have default settings that lead to statistically inappropriate estimates.
Certain estimators that reflect cutting-edge advances in statistics are not yet implemented in R packages for convenient use.
estimatr is designed to solve these problems and provide
estimators tuned for design-based inference.”
The current estimators we provide are:
lm_robust: for fitting linear models with heteroskedasticity/cluster-robust standard errorslm_lin: a wrapper for lm_robust() to simplify interacting centered pre-treatment covariates with a treatment variableiv_robust: two stage least squares estimation of instrumental variables regressiondifference_in_means: for estimating differences in means with appropriate standard errors for unit-randomized, cluster-randomized, block-randomized, matched-pair randomized, and matched-pair clustered designshorvitz_thompson: for estimating average treatment effects taking into consideration treatment probabilities or sampling probabilities for simple and cluster randomized designs
Sonnet, ?
Regressão linear ponderada: estimatr::lm_robust
Dados.media %>%
ggplot2::ggplot(ggplot2::aes(x=PaCO2.mean, y=pH.mean)) +
ggplot2::geom_point() +
ggplot2::stat_smooth(method="lm", formula=y~x, se=TRUE,
ggplot2::aes(weight=Number)) +
ggplot2::labs(x="PaCO2", y="ipH",
title="Weighted Linear Regression") +
ggplot2::theme_test()
Regressão linear ponderada robusta com dados médios:
estimatr::lm_robust
A correlação de Pearson entre participantes entre os valores médios
de pH e PaCO2 pode ser obtida por meio da regressão linear ponderada
robusta com a função estimatr::lm_robust da seguite
forma:
wfit <- estimatr::lm_robust(pH.mean ~ PaCO2.mean,
weights=Number,
data=Dados.media)
print(out <- summary(wfit))
Call:
estimatr::lm_robust(formula = pH.mean ~ PaCO2.mean, data = Dados.media,
weights = Number)
Weighted, Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 7.02590 1.2397 5.6676 0.001298 3.9926 10.0593 6
PaCO2.mean 0.02991 0.2463 0.1214 0.907322 -0.5728 0.6326 6
Multiple R-squared: 0.004866 , Adjusted R-squared: -0.161
F-statistic: 0.01474 on 1 and 6 DF, p-value: 0.9073sunflowerplot(x=Dados.media$PaCO2.mean,
y=Dados.media$pH.mean,
number=Dados.media$Number, seg.col="black",
xlab="PaCO2.mean", ylab="pH.mean",
main="Robust Weighted Linear Regression")
abline(reg=wfit)
rw <- sign(out$coefficients[2])*sqrt(out$r.squared)
cat("Between Pearson's correlation coefficient =", rw, ", p =", wfit$p.value[2])Between Pearson's correlation coefficient = 0.06975858 , p = 0.9073215Regressão linear ponderada robusta:
estimatr::lm_robust
Call:
estimatr::lm_robust(formula = PaCO2 ~ pH, data = Dados, clusters = Subject)
Standard error type: CR2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 6.1988 10.966 0.5653 0.6232 -36.324 48.722 2.248
pH -0.1806 1.508 -0.1197 0.9141 -5.814 5.453 2.357
Multiple R-squared: 0.004253 , Adjusted R-squared: -0.01787
F-statistic: 0.01434 on 1 and 7 DF, p-value: 0.9081Regressão linear robusta INCORRETA:
estimatr::lm_robust
Call:
estimatr::lm_robust(formula = PaCO2 ~ pH, data = Dados)
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 6.1988 4.185 1.4811 0.1455 -2.231 14.6284 45
pH -0.1806 0.576 -0.3135 0.7553 -1.341 0.9796 45
Multiple R-squared: 0.004253 , Adjusted R-squared: -0.01787
F-statistic: 0.0983 on 1 and 45 DF, p-value: 0.7553GEE: Generalized Estimating Equations
As Equações de Estimativa Generalizadas (GEE) foram desenvolvidas para estender o Modelo Linear Generalizado (GLM) para acomodar dados longitudinais e agrupados.
O GEE destina-se a agrupamento simples ou medidas repetidas. Ele não pode acomodar facilmente delineamentos mais complexos, como grupos aninhados ou cruzados; por exemplo, medidas repetidas aninhadas dentro de um sujeito ou grupo. Isso é algo mais adequado para um modelo de efeito misto (GLMM).
Pacotes do R
geepack::geeglm,glmtoolbox::glmgeeyags,MuMIn: QIC for GEE model selection
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fitgeeglm <- geepack::geeglm(PaCO2 ~ pH,
id=Subject,
family=gaussian("identity"),
corstr="exchangeable",
data=Dados)
print(anova(fitgeeglm))Analysis of 'Wald statistic' Table
Model: gaussian, link: identity
Response: PaCO2
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
pH 1 0.33044 0.5654
Call:
geepack::geeglm(formula = PaCO2 ~ pH, family = gaussian("identity"),
data = Dados, id = Subject, corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) 9.6493 8.2509 1.368 0.242
pH -0.6538 1.1374 0.330 0.565
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 0.5574 0.174
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.3205 0.158
Number of clusters: 8 Maximum cluster size: 9 GEE: Regressão logística: geepack::geeglm: artrite
reumatóide
Uma nova droga é testada para o alívio da dor provocada por artrite reumatóide (AR).
O arquivo contém dois graus de alívio (Baixo, Satisfatório) após a aplicação de dois tratamentos (Placebo, Droga ativa) medidos em 3 semanas em 51 pacientes adultos do sexo feminino e masculino. Cada paciente recebeu um dos tratamentos randomizadamente.
Norusis, 2012
Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Tratamento <- factor(Dados$Tratamento)
Dados$Tratamento <- relevel(Dados$Tratamento,
ref="Placebo")
Dados$AlivioDor <- factor(Dados$AlivioDor)
Dados$AlivioDor <- relevel(Dados$AlivioDor,
ref="Baixo")
Dados$Semana <- as.numeric(Dados$Semana)
print(summary(Dados)) Paciente Idade Tratamento Semana
1 : 3 Min. :28.0 Placebo :68 Min. : 5.00
2 : 3 1st Qu.:45.0 Droga ativa:67 1st Qu.: 5.00
3 : 3 Median :55.0 Median : 9.00
4 : 3 Mean :50.9 Mean : 8.82
5 : 3 3rd Qu.:60.0 3rd Qu.:13.00
6 : 3 Max. :65.0 Max. :13.00
(Other):117
Sexo AlivioDor Peso
Length:135 Baixo :36 Min. :2
Class :character Satisfatório:99 1st Qu.:2
Mode :character Median :2
Mean :2
3rd Qu.:2
Max. :2
# A tibble: 6 × 7
Paciente Idade Tratamento Semana Sexo AlivioDor Peso
<fct> <dbl> <fct> <dbl> <chr> <fct> <dbl>
1 1 48 Droga ativa 5 M Satisfatório 2
2 1 48 Droga ativa 9 M Satisfatório 2
3 1 48 Droga ativa 13 M Satisfatório 2
4 2 29 Droga ativa 5 M Satisfatório 2
5 2 29 Droga ativa 9 M Satisfatório 2
6 2 29 Droga ativa 13 M Satisfatório 2# A tibble: 6 × 7
Paciente Idade Tratamento Semana Sexo AlivioDor Peso
<fct> <dbl> <fct> <dbl> <chr> <fct> <dbl>
1 50 58 Droga ativa 5 M Satisfatório 2
2 50 58 Droga ativa 9 M Satisfatório 2
3 50 58 Droga ativa 13 M Satisfatório 2
4 51 37 Placebo 5 M Baixo 2
5 51 37 Placebo 9 M Satisfatório 2
6 51 37 Placebo 13 M Baixo 2Dados <- na.omit(Dados)
fitgeeglm <- try(geepack::geeglm(AlivioDor ~ Tratamento,
id=Paciente,
family=binomial("logit"),
corstr="exchangeable",
data=Dados))Error in lm.fit(zsca, qlf(pr2), offset = soffset) : NA/NaN/Inf in 'y'GEE: Regressão de Poisson: geepack::geeglm:
epilepsia
Para 59 pacientes, o número total de episódios epilépticos num período de referência de 8 semanas foi registrado, sendo que um dos pacientes foi considerado inválido e foi excluído da análise, restando 58 casos.
Após o período de referência, os pacientes foram atribuídos ou ao grupo placebo ou ao grupo da droga ativa, e as contagens dos episódios para os 4 períodos de 2 semanas subsequentes foram obtidos.
O objetivo do estudo é determinar se o tratamento com a droga ativa é efetivo na redução dos episódios epilépticos.
Norusis, 2012
VD
- Número de episódios epilépticos por período: seizure
- log(episódios por semana) = log(#episódios/#semanas)
- Número de semanas: time
- Variável offset: log(time) = logtime
VI de interesse
- Fator: treatment
VI de controle
- Período de tratamento: active
Se o fator treatment diminui o número de episódios epilépticos, i.e., é protetivo, então o RR (risk ratio) do treatment para o nível drug é significante e menor que 1.
Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$treatment <- factor(Dados$treatment)
Dados$treatment <- relevel(Dados$treatment, ref="Placebo")
Dados$active <- factor(Dados$active)
Dados$active <- relevel(Dados$active, ref="Baseline")
print(summary(Dados)) ID age agelog sequence time
1 : 5 Min. :18 Min. :2.89 Length:290 Min. :2.0
2 : 5 1st Qu.:23 1st Qu.:3.14 Class :character 1st Qu.:2.0
3 : 5 Median :28 Median :3.33 Mode :character Median :2.0
4 : 5 Mean :29 Mean :3.33 Mean :3.2
5 : 5 3rd Qu.:35 3rd Qu.:3.56 3rd Qu.:2.0
6 : 5 Max. :57 Max. :4.04 Max. :8.0
(Other):260
logtime active treatment seizure visit4
Min. :0.693 Baseline : 58 Placebo:140 Min. : 0.0 Min. :0.0
1st Qu.:0.693 Treatment:232 Drug :150 1st Qu.: 3.0 1st Qu.:0.0
Median :0.693 Median : 6.0 Median :0.0
Mean :0.970 Mean : 11.5 Mean :0.2
3rd Qu.:0.693 3rd Qu.: 14.0 3rd Qu.:0.0
Max. :2.079 Max. :111.0 Max. :1.0
baseline cseizure logcseizure
Min. : 6.0 Min. : 0.00 Min. :0.000
1st Qu.: 12.0 1st Qu.: 1.50 1st Qu.:0.916
Median : 22.0 Median : 2.38 Median :1.216
Mean : 29.2 Mean : 3.57 Mean :1.264
3rd Qu.: 41.0 3rd Qu.: 4.50 3rd Qu.:1.705
Max. :111.0 Max. :38.00 Max. :3.664
# A tibble: 6 × 13
ID age agelog sequence time logtime active treatment seizure visit4
<fct> <dbl> <dbl> <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
1 1 31 3.43 Baseline 8 2.08 Baseline Placebo 11 0
2 1 31 3.43 Week1 2 0.693 Treatment Placebo 5 0
3 1 31 3.43 Week2 2 0.693 Treatment Placebo 3 0
4 1 31 3.43 Week3 2 0.693 Treatment Placebo 3 0
5 1 31 3.43 Week4 2 0.693 Treatment Placebo 3 1
6 2 30 3.40 Baseline 8 2.08 Baseline Placebo 11 0
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl># A tibble: 6 × 13
ID age agelog sequence time logtime active treatment seizure visit4
<fct> <dbl> <dbl> <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
1 58 36 3.58 Week4 2 0.693 Treatment Drug 0 1
2 59 37 3.61 Baseline 8 2.08 Baseline Drug 12 0
3 59 37 3.61 Week1 2 0.693 Treatment Drug 1 0
4 59 37 3.61 Week2 2 0.693 Treatment Drug 4 0
5 59 37 3.61 Week3 2 0.693 Treatment Drug 3 0
6 59 37 3.61 Week4 2 0.693 Treatment Drug 2 1
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl>fitglmer <- lme4::glmer(seizure ~ treatment*active + (1|ID),
data=Dados,
family=poisson("log"))
print(car::Anova(fitglmer))Analysis of Deviance Table (Type II Wald chisquare tests)
Response: seizure
Chisq Df Pr(>Chisq)
treatment 0.64 1 0.42
active 1667.49 1 < 2e-16 ***
treatment:active 18.97 1 1.3e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: poisson ( log )
Formula: seizure ~ treatment * active + (1 | ID)
Data: Dados
AIC BIC logLik -2*log(L) df.resid
1908 1926 -949 1898 285
Scaled residuals:
Min 1Q Median 3Q Max
-3.872 -0.848 -0.172 0.570 9.894
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.514 0.717
Number of obs: 290, groups: ID, 58
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.11551 0.14111 22.08 < 2e-16 ***
treatmentDrug -0.00815 0.19631 -0.04 0.97
activeTreatment -1.27446 0.04671 -27.28 < 2e-16 ***
treatmentDrug:activeTreatment -0.30239 0.06944 -4.35 1.3e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.5974 0.8782
(Intercept) 2.8331 3.3950
treatmentDrug -0.3990 0.3828
activeTreatment -1.3663 -1.1825
treatmentDrug:activeTreatment -0.4392 -0.1659RR de treatment(drug) = 0.7391 [ 0.6446 0.8472 ]GLMM: Generalized Linear Mixed Model
GLMM é o modelo linear generalizado de efeitos fixos e aleatórios para uma VD.
Fator fixo
- é a população dos níveis do fator;
- contém todos os níveis de interesse do pesquisador;
- efeito de interesse do pesquisador a ser testado estatisticamente.
Fator aleatório
- amostra aleatória dos níveis do fator;
- contém os níveis de interesse do pesquisador com dependência entre os valores da VD dentro de cada nível e independência entre os níveis;
- efeito de controle usado para testar estatisticamente efeito fixo.
Na presença de fator aleatório (agrupamento), a conclusão sobre o efeito do fator fixo se aplica à população do fator aleatório.
O efeito de interação entre efeitos fixo e aleatório é aleatório.
Todo modelo com pelo menos um fator aleatório é misto, pois o intercepto é efeito fixo.
GLMM não exige independência entre as observações e homocedasticidade da VD.
Causas de dependência
Agrupamento (cluster): região, escola, hospital, loja, participante;
Medidas repetidas: mensurações periódicas das mesmas variáveis nas mesmas unidades experimentais.
Distribuição da VD (FD)
- A VD não precisa ser necessariamente normal. GLMM permite VD com distribuição da família exponencial: normal ou gaussiana, binomial (sucesso/fracasso), Poisson (eventos raros: acidente, incidência de doença), gama (valores positivos com assimetria positiva: duração de vida), gaussiana inversa (valores positivos com assimetria positiva: duração de casamento).
Função de ligação (FL)
A FL é função não-linear que estabelece uma relação linear entre a média condicionada da VD e a combinação linear das VI.
A combinação entre a distribuição da VD e a função de ligação resulta em diferentes modelos:
- GLM (ANOVA e Regressão): FD: normal, FL: identidade
- Modelo de regressão logística binária: FD: binomial, FL: logit
- Modelo loglinear: FD: Poisson, FL: log
- Modelo probit: FD: binomial, FL: probit
- Modelo de regressão ordinal: FD: multinomial, FL: logit
Modelos hierárquicos
GLMM permite a integração os dados de níveis hierárquicos na mesma análise.
Níveis de dependência por agrupamento
- Nível 1: dados sobre a unidade observacional
- Nível 2: dados sobre o agrupamento 1
- Nível 3: …
GLMM: Regressão linear: lme4::glmer: pHi e PaCO2
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fitglmer <- lme4::glmer(PaCO2 ~ pH + (1|Subject),
family=gaussian("identity"),
data=Dados)
print(fitglmer)Linear mixed model fit by REML ['lmerMod']
Formula: PaCO2 ~ pH + (1 | Subject)
Data: Dados
REML criterion at convergence: 80.5
Random effects:
Groups Name Std.Dev.
Subject (Intercept) 0.959
Residual 0.442
Number of obs: 47, groups: Subject, 8
Fixed Effects:
(Intercept) pH
17.87 -1.81 Linear mixed model fit by REML ['lmerMod']
Formula: PaCO2 ~ pH + (1 | Subject)
Data: Dados
REML criterion at convergence: 80.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.430 -0.389 0.071 0.518 2.090
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.919 0.959
Residual 0.195 0.442
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 17.867 4.161 4.29
pH -1.805 0.583 -3.10
Correlation of Fixed Effects:
(Intr)
pH -0.997Analysis of Variance Table
npar Sum Sq Mean Sq F value
pH 1 1.87 1.87 9.59Analysis of Deviance Table (Type II Wald chisquare tests)
Response: PaCO2
Chisq Df Pr(>Chisq)
pH 9.59 1 0.002 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
pH 1.87 1.87 1 43.1 9.59 0.0034 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: PaCO2 ~ pH + (1 | Subject)
Data: Dados
REML criterion at convergence: 80.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.430 -0.389 0.071 0.518 2.090
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.919 0.959
Residual 0.195 0.442
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 17.867 4.161 42.857 4.29 9.9e-05 ***
pH -1.805 0.583 43.116 -3.10 0.0034 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: PaCO2
F Df Df.res Pr(>F)
pH 8.56 1 43.4 0.0055 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 numDF denDF F-value p-value
(Intercept) 1 38 211.12 <.0001
pH 1 38 9.59 0.0037Linear mixed-effects model fit by REML
Data: Dados
AIC BIC logLik
88.5 95.73 -40.25
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 0.9587 0.442
Fixed effects: PaCO2 ~ pH
Value Std.Error DF t-value p-value
(Intercept) 17.867 4.161 38 4.294 0.0001
pH -1.805 0.583 38 -3.097 0.0037
Correlation:
(Intr)
pH -0.997
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.43015 -0.38938 0.07101 0.51771 2.09026
Number of Observations: 47
Number of Groups: 8 Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 9.444 17.867 26.2902
pH -2.985 -1.805 -0.6253
Random Effects:
Level: Subject
lower est. upper
sd((Intercept)) 0.5228 0.9587 1.758
Within-group standard error:
lower est. upper
0.3518 0.4420 0.5553 Analysis of Deviance Table (Type II tests)
Response: PaCO2
Chisq Df Pr(>Chisq)
pH 9.59 1 0.002 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1GLMM: Regressão logística: ‘lme4::glmer’: artrite reumatóide
Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Tratamento <- factor(Dados$Tratamento)
Dados$Tratamento <- relevel(Dados$Tratamento,
ref="Placebo")
Dados$AlivioDor <- factor(Dados$AlivioDor)
Dados$AlivioDor <- relevel(Dados$AlivioDor,
ref="Baixo")
Dados$Semana <- as.numeric(Dados$Semana)
print(summary(Dados)) Paciente Idade Tratamento Semana
1 : 3 Min. :28.0 Placebo :68 Min. : 5.00
2 : 3 1st Qu.:45.0 Droga ativa:67 1st Qu.: 5.00
3 : 3 Median :55.0 Median : 9.00
4 : 3 Mean :50.9 Mean : 8.82
5 : 3 3rd Qu.:60.0 3rd Qu.:13.00
6 : 3 Max. :65.0 Max. :13.00
(Other):117
Sexo AlivioDor Peso
Length:135 Baixo :36 Min. :2
Class :character Satisfatório:99 1st Qu.:2
Mode :character Median :2
Mean :2
3rd Qu.:2
Max. :2
# A tibble: 6 × 7
Paciente Idade Tratamento Semana Sexo AlivioDor Peso
<fct> <dbl> <fct> <dbl> <chr> <fct> <dbl>
1 1 48 Droga ativa 5 M Satisfatório 2
2 1 48 Droga ativa 9 M Satisfatório 2
3 1 48 Droga ativa 13 M Satisfatório 2
4 2 29 Droga ativa 5 M Satisfatório 2
5 2 29 Droga ativa 9 M Satisfatório 2
6 2 29 Droga ativa 13 M Satisfatório 2# A tibble: 6 × 7
Paciente Idade Tratamento Semana Sexo AlivioDor Peso
<fct> <dbl> <fct> <dbl> <chr> <fct> <dbl>
1 50 58 Droga ativa 5 M Satisfatório 2
2 50 58 Droga ativa 9 M Satisfatório 2
3 50 58 Droga ativa 13 M Satisfatório 2
4 51 37 Placebo 5 M Baixo 2
5 51 37 Placebo 9 M Satisfatório 2
6 51 37 Placebo 13 M Baixo 2fit <- lme4::glmer(AlivioDor ~ Tratamento + (1|Paciente),
data=Dados,
family=binomial("logit"))
print(car::Anova(fit))Analysis of Deviance Table (Type II Wald chisquare tests)
Response: AlivioDor
Chisq Df Pr(>Chisq)
Tratamento 7.7 1 0.0055 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: AlivioDor ~ Tratamento + (1 | Paciente)
Data: Dados
AIC BIC logLik -2*log(L) df.resid
140.2 148.9 -67.1 134.2 132
Scaled residuals:
Min 1Q Median 3Q Max
-2.160 -0.556 0.221 0.430 1.126
Random effects:
Groups Name Variance Std.Dev.
Paciente (Intercept) 2.43 1.56
Number of obs: 135, groups: Paciente, 50
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.546 0.456 1.20 0.2310
TratamentoDroga ativa 2.137 0.770 2.78 0.0055 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1or <- broom.mixed::tidy(fit,
conf.int=TRUE,
exponentiate=TRUE,
effects="fixed")
cat("OR de Tratamento Droga ativa =", or$estimate[2],
"[", or$conf.low[2], or$conf.high[2], "]")OR de Tratamento Droga ativa = 8.474 [ 1.873 38.33 ]OR de Tratamento Placebo = 1GLMM: Regressão de Poisson: lme4::glmer: epilepsia
Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$treatment <- factor(Dados$treatment)
Dados$treatment <- relevel(Dados$treatment, ref="Placebo")
Dados$active <- factor(Dados$active)
Dados$active <- relevel(Dados$active, ref="Baseline")
print(summary(Dados)) ID age agelog sequence time
1 : 5 Min. :18 Min. :2.89 Length:290 Min. :2.0
2 : 5 1st Qu.:23 1st Qu.:3.14 Class :character 1st Qu.:2.0
3 : 5 Median :28 Median :3.33 Mode :character Median :2.0
4 : 5 Mean :29 Mean :3.33 Mean :3.2
5 : 5 3rd Qu.:35 3rd Qu.:3.56 3rd Qu.:2.0
6 : 5 Max. :57 Max. :4.04 Max. :8.0
(Other):260
logtime active treatment seizure visit4
Min. :0.693 Baseline : 58 Placebo:140 Min. : 0.0 Min. :0.0
1st Qu.:0.693 Treatment:232 Drug :150 1st Qu.: 3.0 1st Qu.:0.0
Median :0.693 Median : 6.0 Median :0.0
Mean :0.970 Mean : 11.5 Mean :0.2
3rd Qu.:0.693 3rd Qu.: 14.0 3rd Qu.:0.0
Max. :2.079 Max. :111.0 Max. :1.0
baseline cseizure logcseizure
Min. : 6.0 Min. : 0.00 Min. :0.000
1st Qu.: 12.0 1st Qu.: 1.50 1st Qu.:0.916
Median : 22.0 Median : 2.38 Median :1.216
Mean : 29.2 Mean : 3.57 Mean :1.264
3rd Qu.: 41.0 3rd Qu.: 4.50 3rd Qu.:1.705
Max. :111.0 Max. :38.00 Max. :3.664
# A tibble: 6 × 13
ID age agelog sequence time logtime active treatment seizure visit4
<fct> <dbl> <dbl> <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
1 1 31 3.43 Baseline 8 2.08 Baseline Placebo 11 0
2 1 31 3.43 Week1 2 0.693 Treatment Placebo 5 0
3 1 31 3.43 Week2 2 0.693 Treatment Placebo 3 0
4 1 31 3.43 Week3 2 0.693 Treatment Placebo 3 0
5 1 31 3.43 Week4 2 0.693 Treatment Placebo 3 1
6 2 30 3.40 Baseline 8 2.08 Baseline Placebo 11 0
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl># A tibble: 6 × 13
ID age agelog sequence time logtime active treatment seizure visit4
<fct> <dbl> <dbl> <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
1 58 36 3.58 Week4 2 0.693 Treatment Drug 0 1
2 59 37 3.61 Baseline 8 2.08 Baseline Drug 12 0
3 59 37 3.61 Week1 2 0.693 Treatment Drug 1 0
4 59 37 3.61 Week2 2 0.693 Treatment Drug 4 0
5 59 37 3.61 Week3 2 0.693 Treatment Drug 3 0
6 59 37 3.61 Week4 2 0.693 Treatment Drug 2 1
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl>fitglmer <- lme4::glmer(seizure ~ treatment*active + (1|ID),
data=Dados,
family=poisson("log"))
print(car::Anova(fitglmer))Analysis of Deviance Table (Type II Wald chisquare tests)
Response: seizure
Chisq Df Pr(>Chisq)
treatment 0.64 1 0.42
active 1667.49 1 < 2e-16 ***
treatment:active 18.97 1 1.3e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: poisson ( log )
Formula: seizure ~ treatment * active + (1 | ID)
Data: Dados
AIC BIC logLik -2*log(L) df.resid
1908 1926 -949 1898 285
Scaled residuals:
Min 1Q Median 3Q Max
-3.872 -0.848 -0.172 0.570 9.894
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.514 0.717
Number of obs: 290, groups: ID, 58
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.11551 0.14111 22.08 < 2e-16 ***
treatmentDrug -0.00815 0.19631 -0.04 0.97
activeTreatment -1.27446 0.04671 -27.28 < 2e-16 ***
treatmentDrug:activeTreatment -0.30239 0.06944 -4.35 1.3e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.5974 0.8782
(Intercept) 2.8331 3.3950
treatmentDrug -0.3990 0.3828
activeTreatment -1.3663 -1.1825
treatmentDrug:activeTreatment -0.4392 -0.1659RR de treatment(drug) = 0.7391 [ 0.6446 0.8472 ]GEE vs. GLMM
Um exemplo de semelhança dos resultados de GEE e GLMM está em:
- Getting Started with Generalized Estimating Equations: https://data.library.virginia.edu/getting-started-with-generalized-estimating-equations/
“Mixed model [GLMM] is close to GEE with exchangeable correlation: Nearly identical coefficients and slightly different SEs.”
Shi, 2012, slide 62
“GEE (population average model):
- On average, is there a trend in score change over time?
- Robust: even if correlation model is wrong, SE still valid
[GLMM] Mixed effects (subject specific model):
- Are there any differences between subjects on the trend in score change over time? (Do all subjects have the same trend over time?)
- Advantages
- Characterization of heterogeneity
- Incomplete unbalanced data”
Shi, 2012, slide 64
GEE x GLMM: Regressão: pHi e PaCO2
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fitgeeglm <- geepack::geeglm(PaCO2 ~ pH,
id=Subject,
family=gaussian("identity"),
corstr="exchangeable",
data=Dados)
print(anova(fitgeeglm))Analysis of 'Wald statistic' Table
Model: gaussian, link: identity
Response: PaCO2
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
pH 1 0.33 0.57
Call:
geepack::geeglm(formula = PaCO2 ~ pH, family = gaussian("identity"),
data = Dados, id = Subject, corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) 9.649 8.251 1.37 0.24
pH -0.654 1.137 0.33 0.57
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 0.557 0.174
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.321 0.158
Number of clusters: 8 Maximum cluster size: 9 Analysis of Variance Table
npar Sum Sq Mean Sq F value
pH 1 1.87 1.87 9.59Linear mixed model fit by REML ['lmerMod']
Formula: PaCO2 ~ pH + (1 | Subject)
Data: Dados
REML criterion at convergence: 80.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.430 -0.389 0.071 0.518 2.090
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.919 0.959
Residual 0.195 0.442
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 17.867 4.161 4.29
pH -1.805 0.583 -3.10Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: PaCO2
F Df Df.res Pr(>F)
pH 8.56 1 43.4 0.0055 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
pH 1.87 1.87 1 43.1 9.59 0.0034 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: PaCO2 ~ pH + (1 | Subject)
Data: Dados
REML criterion at convergence: 80.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.430 -0.389 0.071 0.518 2.090
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.919 0.959
Residual 0.195 0.442
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 17.867 4.161 42.857 4.29 9.9e-05 ***
pH -1.805 0.583 43.116 -3.10 0.0034 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: PaCO2
F Df Df.res Pr(>F)
pH 8.56 1 43.4 0.0055 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1GEE x GLMM: Regressão logística: artrite reumatóide
Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Tratamento <- factor(Dados$Tratamento)
Dados$Tratamento <- relevel(Dados$Tratamento,
ref="Placebo")
Dados$AlivioDor <- factor(Dados$AlivioDor)
Dados$AlivioDor <- relevel(Dados$AlivioDor,
ref="Baixo")
Dados$Semana <- as.numeric(Dados$Semana)
fit <- lme4::glmer(AlivioDor ~ Tratamento + (1|Paciente),
data=Dados,
family=binomial("logit"))
print(car::Anova(fit))Analysis of Deviance Table (Type II Wald chisquare tests)
Response: AlivioDor
Chisq Df Pr(>Chisq)
Tratamento 7.7 1 0.0055 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: AlivioDor ~ Tratamento + (1 | Paciente)
Data: Dados
AIC BIC logLik -2*log(L) df.resid
140.2 148.9 -67.1 134.2 132
Scaled residuals:
Min 1Q Median 3Q Max
-2.160 -0.556 0.221 0.430 1.126
Random effects:
Groups Name Variance Std.Dev.
Paciente (Intercept) 2.43 1.56
Number of obs: 135, groups: Paciente, 50
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.546 0.456 1.20 0.2310
TratamentoDroga ativa 2.137 0.770 2.78 0.0055 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1or <- broom.mixed::tidy(fit,
conf.int=TRUE,
exponentiate=TRUE,
effects="fixed")
cat("OR de Tratamento Droga ativa =", or$estimate[2],
"[", or$conf.low[2], or$conf.high[2], "]")OR de Tratamento Droga ativa = 8.47 [ 1.87 38.3 ]OR de Tratamento Placebo = 1Dados <- na.omit(Dados)
fitgeeglm <- try(geepack::geeglm(AlivioDor ~ Tratamento,
id=Paciente,
family=binomial("logit"),
corstr="exchangeable",
data=Dados))Error in lm.fit(zsca, qlf(pr2), offset = soffset) : NA/NaN/Inf in 'y'GEE x GLMM: Regressão de Poisson: epilepsia
Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$treatment <- factor(Dados$treatment)
Dados$treatment <- relevel(Dados$treatment, ref="Placebo")
Dados$active <- factor(Dados$active)
Dados$active <- relevel(Dados$active, ref="Baseline")
fitgeeglm <- geepack::geeglm(seizure ~ treatment*active,
offset=logtime,
id=ID,
family=poisson("log"),
corstr="exchangeable",
data=Dados)
print(anova(fitgeeglm))Analysis of 'Wald statistic' Table
Model: poisson, link: log
Response: seizure
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
treatment 1 0.104 0.747
active 1 0.080 0.778
treatment:active 1 3.125 0.077 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
geepack::geeglm(formula = seizure ~ treatment * active, family = poisson("log"),
data = Dados, offset = logtime, id = ID, corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) 1.348 0.157 73.34 <2e-16 ***
treatmentDrug -0.107 0.194 0.30 0.581
activeTreatment 0.112 0.116 0.93 0.335
treatmentDrug:activeTreatment -0.302 0.171 3.12 0.077 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 10.4 2.28
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.598 0.0811
Number of clusters: 58 Maximum cluster size: 5 rr <- broom.mixed::tidy(fitgeeglm,
conf.int=TRUE,
exponentiate=TRUE,
effects="fixed")
cat("RR de treatment(drug) =", rr$estimate[4],
"[", rr$conf.low[4], rr$conf.high[4], "]")RR de treatment(drug) = 0.739 [ 0.529 1.03 ]fitglmer <- lme4::glmer(seizure ~ treatment*active + (1|ID),
data=Dados,
family=poisson("log"))
print(car::Anova(fitglmer))Analysis of Deviance Table (Type II Wald chisquare tests)
Response: seizure
Chisq Df Pr(>Chisq)
treatment 0.64 1 0.42
active 1667.49 1 < 2e-16 ***
treatment:active 18.97 1 1.3e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: poisson ( log )
Formula: seizure ~ treatment * active + (1 | ID)
Data: Dados
AIC BIC logLik -2*log(L) df.resid
1908 1926 -949 1898 285
Scaled residuals:
Min 1Q Median 3Q Max
-3.872 -0.848 -0.172 0.570 9.894
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.514 0.717
Number of obs: 290, groups: ID, 58
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.11551 0.14111 22.08 < 2e-16 ***
treatmentDrug -0.00815 0.19631 -0.04 0.97
activeTreatment -1.27446 0.04671 -27.28 < 2e-16 ***
treatmentDrug:activeTreatment -0.30239 0.06944 -4.35 1.3e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Computing profile confidence intervals ... 2.5 % 97.5 %
.sig01 0.597 0.878
(Intercept) 2.833 3.395
treatmentDrug -0.399 0.383
activeTreatment -1.366 -1.182
treatmentDrug:activeTreatment -0.439 -0.166RR de treatment(drug) = 0.739 [ 0.645 0.847 ]Referências
- Baigorri, F & Calvet, X & Joseph, D (1997) Equipment review: Gastric intramucosal pH measurement. Crit Care 1: 61. https://doi.org/10.1186/cc104.
- Bakdash, JZ & Marusich, LR (2017) Repeated measures correlation. Front. Psychol. 8:456. doi: 10.3389/fpsyg.2017.00456. https://osf.io/djphm/
- Bland, JM & Altman, DG (1995a) Calculating correlation coefficients with repeated observations: Part 1 - Correlation within subjects. BMJ 310(6977): 446. https://doi.org/10.1136/bmj.310.6977.446
- Bland JM & Altman DG (1995b) Calculating correlation coefficients with repeated observations: Part 2 - Correlation between subjects. BMJ 310(6980): 633. doi: 10.1136/bmj.310.6980.633. Erratum in: BMJ 1996 312(7030):572.
- Chartier, S & Faulkner, A (2008) General linear models: An integrated approach to statistics. Tutorial in Quantitative Methods for Psychology 4(2), 65-78.
- Clark, M (?) Mixed model with R. https://m-clark.github.io/mixed-models-with-R/.
- Duursma, R & Powell, J (2016) Linear mixed-effects models. http://www.hiercourse.com/mixedeffects.
- Duursma, R & Powell, J (2016) Linear mixed-effects models: Solutions. http://www.hiercourse.com/mixedeffects.
- Faraway, JJ (2016) Extending the linear model with R: Generalized linear, mixed effects and nonparametric regression models. 2 ed. USA: CRC.
- Francq, GF et al. (2019) Confidence, prediction, and tolerance in linear mixed models. Statistics in Medicine 38(30): 5603-22.
- Francq, GF et al. (2020) Erratum: Confidence, prediction, and tolerance in linear mixed models. Statistics in Medicine 39(14): 2015-2015.
- Hamaker, EL (2013) Why researchers should think “within-person”: A paradigmatic rationale. In Mehl MR &, Conner TS (Eds.) (2013) Handbook of research methods for studying daily life. NY: Guilford. p. 43-61.
- Hothorn, T & Everitt, BS (2014) A handbook of statistical analyses using R. 3rd ed. Chapman & Hall/CRC, 2014. https://rdrr.io/cran/HSAUR3/. R package HSAUR3.
- Judkins DR & Porter KE (2016) Robustness of ordinary least squares in randomized clinical trials. Stat Med. 35(11):1763-73. doi: 10.1002/sim.6839.
- Kirkwood, BR & Sterne, JAC (2003) Essential medical statistics. 2nd ed. USA: Blackwell. Chapter 13: Analysis of clustered data.
- Konrad, M (2021) Clustered standard errors with R. https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/.
- Laird, NM and Ware, JH (1982) Random-effects models for longitudinal data. Biometrics 38(4): 963-74.
- Magnusson, K (2015) Using R and lme/lmer to fit different two- and three-level longitudinal models. https://rpsychologist.com/r-guide-longitudinal-lme-lmer.
- Mixed effects logistic regression: R data analysis examples. https://stats.oarc.ucla.edu/r/dae/mixed-effects-logistic-regression/.
- Norusis, MJ (2012) IBM SPSS statistics 19 advanced statistical procedures companion. NJ: Prentice-Hall.
- Oikawa, KF (2019) Análise estatística de dados longitudinais e hierárquicos em psicologia: uma análise comparativa entre GEE e GLMM. Tese de Doutorado, Psicologia Experimental, Instituto de Psicologia, Universidade de São Paulo, São Paulo. doi:10.11606/T.47.2019.tde-18112019-183159. Recuperado em 2022-05-08, de www.teses.usp.br.
- Oliveira, AA (Coord.) (?) Modelos generalizados mistos (GLMM). Laboratório de Ecologia de Florestas Tropicais do Instituto de Biocências da USP. http://labtrop.ib.usp.br/doku.php?id=cursos:planeco:roteiro:12-glmm.
- Pekár, S & Brabec, M (2017) Generalized estimating equations: A pragmatic and flexible approach to the marginal GLM modelling of correlated data in the behavioural sciences. Ethology 124:86–93. DOI: 10.1111/eth.12713.
- Sainani, K. (2010) The importance of accounting for correlated observations. Physical Medicine and Rehabilitation 2(9):858-61. doi: 10.1016/j.pmrj.2010.07.482.
- Shi, L (2012) Application of GEE and Mixed Effects model in longitudinal data analysis. Umass Boston: Dana-Farber/Harvard Cancer Center.
- Sonnet, L (?) Getting started using estimatr. https://declaredesign.org/r/estimatr/articles/getting-started.html.
- Xu, C et al. (2019) An R package for model fitting, model selection and the simulation for longitudinal data with dropout missingness. Commun Stat Simul Comput 48(9): 2812–29. doi:10.1080/03610918.2018.1468457.
- Winter, B & Wieling, M (2016) How to analyze linguistic change using mixed models, growth curve analysis and generalized additive modeling. Journal of Language Evolution 1(1): 7–18. https://doi.org/10.1093/jole/lzv003.
- Zeger, SL, Liang, KY & Albert, PS (1988) Models for longitudinal data: A generalized estimating equation approach. Biometrics 44(4): 1049-60.
- Zuur, A et al. (2009) Mixed effects models and extensions in ecology with R. NY: Springer.