GEE (Generalized Estimating Equations) e
GLMM (Generalized Linear Mixed Model)
em R
MAN5728: Tópicos Essenciais da
Bioestatística
MÓDULO 8 - ESTATÍSTICA APLICADA II:
ANÁLISE E
INTERPRETAÇÃO DE DADOS
23 fevereiro 2026 21:06h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(broom.mixed, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(estimatr, warn.conflicts=FALSE))
suppressMessages(library(gee, warn.conflicts=FALSE))
suppressMessages(library(geepack, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(ggpubr, warn.conflicts=FALSE))
suppressMessages(library(HSAUR3, warn.conflicts=FALSE))
suppressMessages(library(IDPmisc, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lme4, warn.conflicts=FALSE))
suppressMessages(library(lmerTest, warn.conflicts=FALSE))
suppressMessages(library(lsr, warn.conflicts=FALSE))
suppressMessages(library(magrittr, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(purrr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(rmcorr, warn.conflicts=FALSE))
suppressMessages(library(rsample, warn.conflicts=FALSE))
suppressMessages(library(rstatix, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))
suppressMessages(library(weights, warn.conflicts=FALSE))Material
- HTML de Rmarkdown em
RPubs - Prof. José Siqueira: ResearchGate
Adicionar
cor_multilevelmisty::multilevel.corContextual Effects Model (CEM): https://psychbruce.github.io/supp/CEM
Márcio Braga de Melo, Dimitri Daldegan-Bueno, Maria Gabriela Menezes Oliveira, Altay Lino de Souza (2022) Beyond ANOVA and MANOVA for repeated measures: Advantages of generalized estimated equations and generalized linear mixed models and its use in neuroscience research. European Journal of Neuroscience 56(12): 6089-98.
Real longitudinal data analysis for real people - Building a good enough mixed model - Cheng et al - 2009
How to analyze linguistic change using mixed models, Growth Curve Analysis and Generalized Additive Modeling - Winter & Wieling - 2016
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.
Conforme Kirkwood & Sterne, 2003, p. 355-6,
“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.”
Conforme Kirkwood & Sterne, 2003, p. 369-70,
“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.”
Análise de dados longitudinais
Conforme Hamaker (2013),
“Para ilustrar o que pode acontecer quando tentamos generalizar para o geral os resultados obtidos para o agregado, considere o seguinte exemplo. Suponha que estamos interessados na relação entre a velocidade de digitação (isto é, número de palavras digitadas por minuto) e a porcentagem de erros cometidos. Se observarmos a relação transversal (isto é, no nível populacional), poderemos encontrar uma relação negativa, no sentido de que pessoas que digitam mais rápido cometem menos erros (isso pode refletir o fato de que pessoas com habilidades de digitação mais desenvolvidas e mais experiência digitam mais rápido e cometem menos erros). Essa relação negativa entre pessoas é mostrada no painel esquerdo da Figura 3.1. Se generalizarmos esse resultado para o nível dentro da pessoa, concluiríamos que, se um indivíduo específico digita mais rápido, ele cometerá menos erros. Claramente, isso não é o que esperamos: na verdade, estamos razoavelmente certos de que, para qualquer indivíduo em particular, o número de erros aumentará se ele tentar digitar mais rápido. Isso implica uma relação positiva — e não negativa — no nível dentro da pessoa. No painel direito da Figura 3.1, a relação dentro da pessoa entre essas duas variáveis é mostrada para vários indivíduos. Também se observa que, como as médias individuais diferem na direção oposta e a variabilidade entre médias individuais é muito maior do que a variabilidade dentro das pessoas, a abordagem padrão de grandes amostras baseada em uma seção transversal resultará em uma relação negativa (ver também Nezlek, 2001). No exemplo anterior, todos são caracterizados pela mesma relação dentro da pessoa entre a porcentagem de erros e a velocidade de digitação, e isso pode, portanto, ser interpretado como uma lei geral, que diverge da relação encontrada transversalmente. Entretanto, também é possível que as relações dentro da pessoa sejam caracterizadas por diferenças individuais, de modo que nenhuma lei geral se aplique a cada pessoa. De qualquer forma, o ponto central é que estudar grandes amostras representativas — como é prática comum na pesquisa psicológica — não nos fornece informação sobre leis gerais. Assim, qualificar a pesquisa com grandes amostras como nomotética e, portanto, indiretamente, como mais científica do que outras abordagens por estar voltada à descoberta de leis gerais, é equivocado.”
Hamaker, 2013, in Mehl & Conner, p. 43-61.
Hamaker, 2013, in Mehl & Conner, p. 43-61.
Sainani, 2010, p. 861.
Estrutura 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 (wide).
O formato longo (long) é 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 associação linear 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
Conforme Bakdash & Marusich (2017),
Resumo: A correlação para medidas repetidas (rmcorr) é uma técnica estatística para determinar a associação comum dentro do indivíduo para medidas pareadas avaliadas em duas ou mais ocasiões em múltiplos indivíduos. Regressão/correlação simples é frequentemente aplicada a observações não independentes ou a dados agregados; isso pode produzir resultados enviesados e espúrios devido à violação da independência e/ou a padrões diferentes entre participantes versus dentro dos participantes. Diferentemente da regressão/correlação simples, o rmcorr não viola a suposição de independência das observações. Além disso, o rmcorr tende a ter poder estatístico muito maior porque nem a média nem a agregação são necessárias para uma questão de pesquisa intraindividual. O rmcorr estima a inclinação comum da regressão, a associação compartilhada entre os indivíduos. Para tornar o rmcorr acessível, fornecemos informações de base sobre suas suposições e equações, visualização, poder e compromissos em comparação com a modelagem multinível. Introduzimos o pacote R (
rmcorr) e demonstramos seu uso para estatística inferencial e visualização com dois conjuntos de dados de exemplo. Os exemplos são usados para ilustrar questões de pesquisa em diferentes níveis de análise, intraindividual e interindividual. O rmcorr é adequado para questões de pesquisa sobre a associação linear comum em dados de medidas repetidas pareadas. Todos os resultados são totalmente reprodutíveis.
Palavras-chave: correlação, medidas repetidas, diferenças individuais, intraindividual, poder estatístico, modelagem multinível
Correlação intraparticipantes: 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).
Conforme Baigorri et al., 1997,
“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.”
Conforme figura abaixo, a relação é inversa entre pHi e PaCO₂. Quando a PaCO₂ aumenta, há mais CO₂ dissolvido nos tecidos, o que leva à formação de mais ácido carbônico e liberação de íons hidrogênio. Isso torna o meio mais ácido e o pH intramucoso cai.
Portanto, aumento da PaCO₂ reduz o pH intramucoso; queda da PaCO₂ eleva o pH intramucoso. Trata-se do efeito fisiológico clássico da retenção de CO₂ produzindo acidificação.
Baigorri et al., 1997: A tonometria gástrica determina PaCO2 intraluminal, que se presume estar em equilíbrio com o PaCO2 na mucosa gástrica. 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.
Os dados usados para testar esta associação estão disponíveis em
Bland & Altman (1995a) baseado em Boyd et al. (1993) e no pacote R
em rmcorr::bland1995.
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
knitr::kable(Dados, digits = 2, align = "c")| Subject | pH | PaCO2 |
|---|---|---|
| 1 | 6.68 | 3.97 |
| 1 | 6.53 | 4.12 |
| 1 | 6.43 | 4.09 |
| 1 | 6.33 | 3.97 |
| 2 | 6.85 | 5.27 |
| 2 | 7.06 | 5.37 |
| 2 | 7.13 | 5.41 |
| 2 | 7.17 | 5.44 |
| 3 | 7.40 | 5.67 |
| 3 | 7.42 | 3.64 |
| 3 | 7.41 | 4.32 |
| 3 | 7.37 | 4.73 |
| 3 | 7.34 | 4.96 |
| 3 | 7.35 | 5.04 |
| 3 | 7.28 | 5.22 |
| 3 | 7.30 | 4.82 |
| 3 | 7.34 | 5.07 |
| 4 | 7.36 | 5.67 |
| 4 | 7.33 | 5.10 |
| 4 | 7.29 | 5.53 |
| 4 | 7.30 | 4.75 |
| 4 | 7.35 | 5.51 |
| 5 | 7.35 | 4.28 |
| 5 | 7.30 | 4.44 |
| 5 | 7.30 | 4.32 |
| 5 | 7.37 | 3.23 |
| 5 | 7.27 | 4.46 |
| 5 | 7.28 | 4.72 |
| 5 | 7.32 | 4.75 |
| 5 | 7.32 | 4.99 |
| 6 | 7.38 | 4.78 |
| 6 | 7.30 | 4.73 |
| 6 | 7.29 | 5.12 |
| 6 | 7.33 | 4.93 |
| 6 | 7.31 | 5.03 |
| 6 | 7.33 | 4.93 |
| 7 | 6.86 | 6.85 |
| 7 | 6.94 | 6.44 |
| 7 | 6.92 | 6.52 |
| 8 | 7.19 | 5.28 |
| 8 | 7.29 | 4.56 |
| 8 | 7.21 | 4.34 |
| 8 | 7.25 | 4.32 |
| 8 | 7.20 | 4.41 |
| 8 | 7.19 | 3.69 |
| 8 | 6.77 | 6.09 |
| 8 | 6.82 | 5.58 |
Correlação intraparticipantes global
A estimativa do valor absoluto da correlação de Pearson para medidas repetidas é a raiz quadrada do eta parcial ao quadrado de PaCO₂ obtida pelo modelo linear geral (GLM).
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.
O coeficiente de correlação para medidas repetidas \(r_{\text{rm}} \approx −0.5, \, p <0.001\) (significante) indica associação linear inversa de magnitude moderada entre PaCO₂ e pHi dentro dos mesmos pacientes ao longo do tempo.
| Método | \(r_{\text{rm}}\) | IC95% inf | IC95% sup | Valor p |
|---|---|---|---|---|
GLM: rmcorr::rmcorr |
-0.507 | -0.707 | -0.232 | 8.471e-04 |
GLM: lm |
-0.507 | -0.680 | -0.234 | 8.471e-04 |
GLMM: lmerTest::lmer |
-0.484 | -0.660 | -0.214 | 1.087e-03 |
Em termos clínicos: aumento da PaCO₂ tende a reduzir o pHi, e redução da PaCO₂ tende a elevar o pHi, em concordância com a fisiologia.
## =========================================================
## Tabela: r_rm e IC95% (rmcorr vs GLM vs GLMM) — ordem correta + p científico
## =========================================================
alfa <- 0.05
pkgs <- c("rmcorr", "lmerTest", "lme4", "effectsize", "knitr")
for (p in pkgs) {
if (!requireNamespace(p, quietly = TRUE)) stop("Pacote ausente: ", p)
suppressPackageStartupMessages(library(p, character.only = TRUE))
}
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fmt_p <- function(p) formatC(p, format = "e", digits = 3)
get_rmcorr_ci <- function(fit) {
if (!is.null(fit$CI) && length(fit$CI) >= 2) return(c(lwr = fit$CI[1], upr = fit$CI[2]))
if (!is.null(fit$ci) && length(fit$ci) >= 2) return(c(lwr = fit$ci[1], upr = fit$ci[2]))
s <- capture.output(print(fit))
i <- grep("confidence interval", s, ignore.case = TRUE)
if (length(i)) {
nums <- regmatches(s[i[1]], gregexpr("-?\\d+\\.\\d+|-?\\d+", s[i[1]]))[[1]]
nums <- suppressWarnings(as.numeric(nums))
nums <- nums[is.finite(nums)]
if (length(nums) >= 2) return(c(lwr = nums[1], upr = nums[2]))
}
c(lwr = NA_real_, upr = NA_real_)
}
extract_eta_ci_row <- function(eta_tbl, term) {
et <- as.data.frame(eta_tbl)
row <- et[et$Parameter == term, , drop = FALSE]
eta_col <- names(row)[grepl("Eta2", names(row), ignore.case = TRUE) &
grepl("partial", names(row), ignore.case = TRUE)][1]
if (is.na(eta_col)) {
cand <- names(row)[grepl("Eta2", names(row), ignore.case = TRUE)]
cand <- cand[sapply(row[cand], function(x) is.numeric(x) && all(is.finite(x)))]
eta_col <- cand[1]
}
eta <- as.numeric(row[[eta_col]])
low_col <- names(row)[grepl("CI.*low|CI_low|CI\\.low", names(row), ignore.case = TRUE)][1]
upr_col <- names(row)[grepl("CI.*high|CI_high|CI\\.high|CI.*up", names(row), ignore.case = TRUE)][1]
if (!is.na(low_col) && !is.na(upr_col)) {
lwr <- as.numeric(row[[low_col]])
upr <- as.numeric(row[[upr_col]])
return(list(eta = eta, lwr = lwr, upr = upr))
}
if ("CI" %in% names(row)) {
nums <- regmatches(as.character(row$CI),
gregexpr("-?\\d+\\.\\d+|-?\\d+", as.character(row$CI)))[[1]]
nums <- suppressWarnings(as.numeric(nums))
nums <- nums[is.finite(nums)]
if (length(nums) >= 2) return(list(eta = eta, lwr = nums[1], upr = nums[2]))
}
list(eta = eta, lwr = NA_real_, upr = NA_real_)
}
r_from_eta <- function(eta, beta) sign(beta) * sqrt(eta)
## ---------- 1) rmcorr ----------
fit_rm <- rmcorr::rmcorr(participant = Subject,
measure1 = PaCO2,
measure2 = pH,
dataset = Dados,
CI.level = 1 - alfa)
r_rm <- fit_rm$r
p_rm <- fit_rm$p
ci_rm <- get_rmcorr_ci(fit_rm)
## ---------- 2) GLM ----------
modelo <- lm(pH ~ Subject + PaCO2, data = Dados)
p_lm <- summary(modelo)$coefficients["PaCO2", "Pr(>|t|)"]
b_lm <- coef(modelo)["PaCO2"]
eta_lm <- effectsize::eta_squared(modelo,
alternative = "two.sided",
ci = 1 - alfa,
partial = TRUE)
lm_row <- extract_eta_ci_row(eta_lm, "PaCO2")
r_lm <- r_from_eta(lm_row$eta, b_lm)
r_lm_lwr <- r_from_eta(lm_row$lwr, b_lm)
r_lm_upr <- r_from_eta(lm_row$upr, b_lm)
## corrigir ordem
ci_lm_sorted <- sort(c(r_lm_lwr, r_lm_upr))
r_lm_lwr <- ci_lm_sorted[1]
r_lm_upr <- ci_lm_sorted[2]
## ---------- 3) GLMM ----------
m <- lmerTest::lmer(pH ~ PaCO2 + (1 | Subject), data = Dados)
p_lmer <- summary(m)$coefficients["PaCO2", "Pr(>|t|)"]
b_lmer <- lme4::fixef(m)["PaCO2"]
eta_lmer <- effectsize::eta_squared(m,
alternative = "two.sided",
ci = 1 - alfa,
partial = TRUE)
lmer_row <- extract_eta_ci_row(eta_lmer, "PaCO2")
r_lmer <- r_from_eta(lmer_row$eta, b_lmer)
r_lmer_lwr <- r_from_eta(lmer_row$lwr, b_lmer)
r_lmer_upr <- r_from_eta(lmer_row$upr, b_lmer)
## corrigir ordem
ci_lmer_sorted <- sort(c(r_lmer_lwr, r_lmer_upr))
r_lmer_lwr <- ci_lmer_sorted[1]
r_lmer_upr <- ci_lmer_sorted[2]
## ---------- tabela ----------
tab <- data.frame(
"Method" = c("GLM: rmcorr::rmcorr",
"GLM: lm",
"GLMM: lmerTest::lmer"),
"r_rm" = round(c(r_rm, r_lm, r_lmer), 3),
"CI95% lwr" = round(c(ci_rm["lwr"], r_lm_lwr, r_lmer_lwr), 3),
"CI95% upr" = round(c(ci_rm["upr"], r_lm_upr, r_lmer_upr), 3),
"p-value" = c(fmt_p(p_rm), fmt_p(p_lm), fmt_p(p_lmer)),
check.names = FALSE
)
knitr::kable(tab, align = "lrrrr")| Method | r_rm | CI95% lwr | CI95% upr | p-value |
|---|---|---|---|---|
| GLM: rmcorr::rmcorr | -0.507 | -0.707 | -0.232 | 8.471e-04 |
| GLM: lm | -0.507 | -0.680 | -0.234 | 8.471e-04 |
| GLMM: lmerTest::lmer | -0.484 | -0.660 | -0.214 | 1.087e-03 |
Correlação intraparticipantes individual: IC95% por bootstrapping
A estimação por bootstrapping (reamostragem) das correlações intraparticipantes tem utilidade inferencial quando \(n\) por indivíduo é pequeno, a distribuição de \(r\) é assimétrica e as suposições paramétricas são frágeis. O reamostramento empírico aproxima a distribuição amostral de cada correlação individual sem normalidade, fornecendo IC95% robustos e menos dependentes de modelo.
O método não estima \(r_{\mathrm{rm}}\) porque não combina sujeitos nem modela explicitamente a dependência entre medidas; ele quantifica a associação dentro de cada participante. Sua vantagem é revelar heterogeneidade interindividual: alguns sujeitos exibem correlação fortemente negativa, outros próxima de zero ou positiva, o que seria ocultado por uma única medida agregada.
Os IC bootstrap permitem avaliar estabilidade e incerteza de cada \(r_i\): intervalos largos indicam baixa informação (poucas observações/alta variabilidade), enquanto intervalos concentrados indicam associação mais estável. A mediana bootstrap é menos sensível a outliers do que a estimativa pontual.
Embora conceitualmente distinto, o padrão obtido é coerente com as
correlações intraparticipantes do rmcorr::rmcorr, isto é,
predominância de associações negativas, corroborando a relação inversa
PaCO₂–pHi ao nível individual, sem recorrer a pressupostos paramétricos
fortes. Note que \(r_7\) não foi
estimada por reamostragem.
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.58## 1) Bootstrap por Subject (sem piping)
boot_corr <- Dados
boot_corr <- tidyr::nest(boot_corr, data = -c(Subject))
boot_corr <- dplyr::mutate(
boot_corr,
boots = purrr::map(
data,
function(df) {
rsample::bootstraps(df, times = 1e2, apparent = FALSE)
}
)
)Warning: There were 5 warnings in `dplyr::mutate()`.
The first warning was:
ℹ In argument: `boots = purrr::map(...)`.
Caused by warning in `rsample::bootstraps()`:
! Some assessment sets contained zero rows.
ℹ Run `dplyr::last_dplyr_warnings()` to see the 4 remaining warnings.boot_corr <- tidyr::unnest(boot_corr, cols = c(boots))
boot_corr <- dplyr::mutate(
boot_corr,
correlations = purrr::map(
splits,
function(sp) {
rstatix::cor_test(PaCO2, pH, data = rsample::analysis(sp))
}
)
)
corr <- boot_corr
corr <- tidyr::unnest(corr, cols = c(correlations))
corr <- dplyr::select(corr, -data, -splits, -id)
corr <- IDPmisc::NaRV.omit(corr)
## 2) Correlação pontual por Subject (sem piping)
point <- Dados
point <- dplyr::group_by(point, Subject)
point <- rstatix::cor_test(point, PaCO2, pH, method = "pearson")
## 3) IC bootstrap por Subject (sem piping)
CI <- corr
CI <- dplyr::group_by(CI, Subject)
CI <- dplyr::summarise(
CI,
lwr_CI = stats::quantile(cor, alfa/2, na.rm = TRUE),
.estimate = stats::median(cor, na.rm = TRUE),
upr_CI = stats::quantile(cor, 1 - alfa/2, na.rm = TRUE),
.groups = "drop"
)
knitr::kable(CI, digits = 3, align = "crrr")| Subject | lwr_CI | .estimate | upr_CI |
|---|---|---|---|
| 1 | -0.880 | -0.053 | 0.955 |
| 2 | 1.000 | 1.000 | 1.000 |
| 3 | -0.930 | -0.470 | 0.533 |
| 4 | -0.350 | 0.620 | 1.000 |
| 5 | -0.936 | -0.655 | 0.615 |
| 6 | -1.000 | -0.390 | 0.500 |
| 8 | -0.970 | -0.820 | 0.544 |
Correlação intraparticipantes global por bootstrapping
A estimação de \(r_{\mathrm{rm}}\) por bootstrapping é útil porque não depende de normalidade, lida melhor com \(n\) pequeno e com assimetria/heterocedasticidade, e fornece inferência robusta para a correlação dentro dos participantes mantendo a estrutura de medidas repetidas. O reamostramento por participante preserva a dependência intraindivíduo e produz IC baseados na distribuição empírica de \(r_{\mathrm{rm}}\), reduzindo sensibilidade a violações de modelo e a outliers.
Comparando os IC95% analítico \([−0.711, −0.223]\) e bootstrap \([−0.727, −0.070]\), há diferença numérica esperada porque os métodos usam aproximações distintas: o analítico depende de suposições paramétricas e linearização; o bootstrap usa quantis da distribuição reamostrada. Apesar disso, a conclusão é qualitativamente idêntica: ambos os intervalos são inteiramente negativos, excluem 0 e indicam associação inversa significativa entre PaCO₂ e pHi dentro dos participantes. O IC bootstrap é assimétrico e com maior amplitude.
# 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.7298604 -0.0720461 Correlação intraparticipantes por cor.test:
incorreta
A correlação de Pearson ingênua (\(r = −0.065, \,p = 0.663\)) é não significante e praticamente nula, contrariando a expectativa fisiológica de relação inversa entre PaCO₂ e pHi.
O problema é metodológico: o cálculo ignora a estrutura longitudinal/clusterizada, mistura variação entre pacientes com variação dentro do mesmo paciente, viola a independência e dilui a associação verdadeira. Diferenças basais de pHi e PaCO₂ entre indivíduos e heterogeneidade intraindividual mascaram a relação fisiológica.
Quando a dependência é tratada corretamente (correlação para medidas repetidas), aparece a associação inversa esperada (\(r_{\text{rm}} = −0.507, \, p = 0.000847\) (significante)). Portanto, o resultado não significante da Pearson simples é artefato de especificação incorreta do modelo, não evidência contra a fisiologia.
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 não é:
a média aritmética das correlações dos participantes;
a média ponderada das correlações dos participantes.
Bland & Altman (1995a) mostraram a maneira correta de estimar a correlação intraparticipantes.
Bland & Altman, 1995a
Correlação entre participantes por
weights::wtd.cor
Esse método estima uma correlação entre participantes, usando médias por sujeito. Ele não é adequado para recuperar a relação fisiológica PaCO₂–pHi porque responde a outra pergunta estatística.
Primeiro, ao agregar por médias, elimina-se a variação dentro do participante, que é justamente onde a fisiologia atua. A relação esperada é intraparticipante (mudanças de PaCO₂ acompanhadas de mudanças opostas de pHi no mesmo indivíduo). A correlação entre médias compara indivíduos entre si, não a dinâmica dentro de cada um.
Segundo, há confusão por heterogeneidade basal. Diferenças individuais de pHi e PaCO₂ (set-point, gravidade, compensação, erros de medida) introduzem variação entre sujeitos que não reflete o mecanismo ácido–base, podendo inverter ou anular o sinal.
Terceiro, ocorre viés ecológico (falácia ecológica). A associação entre médias pode ser próxima de zero ou até positiva mesmo quando a associação verdadeira dentro do indivíduo é negativa. O resultado obtido (\(r = 0.07, p = 0.90\)) ilustra essa diluição.
Quarto, ponderar pelo número de observações não corrige o problema estrutural: continua sendo uma correlação entre médias, sem modelar dependência nem decompor variação intra vs. intersujeitos.
Conclusão: a correlação ponderada entre médias é inadequada para
recuperar o resultado fisiológico esperado. Métodos que preservam a
estrutura longitudinal (rmcorr, modelos mistos) estimam a
associação intraparticipante correta e revelam a relação inversa
PaCO₂–pHi.
# Bland & Altman 1995b: analitica e boostrapping
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 8print(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.004277 0.5847 -0.007316 0.8984Desenhando 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
Conforme Judkins & Porter, 2016,
“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.”
Conforme Sonnet, ?,
“
estimatris 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.estimatris 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”
Regressão linear robusta ponderada (com dados médios):
estimatr::lm_robust (HC2)
Essa abordagem modela apenas diferenças entre participantes usando valores médios. Fisiologicamente, isso é inadequado porque o mecanismo ácido–base opera dentro do mesmo indivíduo: variações de PaCO₂ ao longo do tempo produzem variações opostas de pHi no próprio paciente. Ao usar médias por sujeito, a dinâmica intraparticipante é removida, e a regressão passa a comparar indivíduos entre si, não a resposta fisiológica.
A heterogeneidade basal entre pacientes (set-point de ventilação, estado ácido–base, compensações metabólicas, gravidade clínica, erro de medida) domina a variação entre médias e pode anular ou inverter o sinal esperado. O coeficiente estimado próximo de zero e positivo (\(\beta = 0.03,\, r = 0.07,\, p = 0.91\)) não reflete ausência de efeito fisiológico, mas diluição por mistura de variabilidade interindividual não causal com a variabilidade relevante.
Ponderação pelo número de observações e erro-padrão robusto (HC2) tratam apenas eficiência/heterocedasticidade; não recuperam a relação fisiológica porque o modelo continua especificado no nível errado (entre participantes). O resultado é, portanto, fisiologicamente inadequado: não capta a relação inversa PaCO₂–pHi que se manifesta intraparticipante. Modelos que preservam a estrutura longitudinal (correlação para medidas repetidas, modelos mistos) são os apropriados para recuperar o mecanismo esperado.
A correlação de Pearson entre participantes entre os valores médios
de pHi e PaCO₂ pode ser obtida por meio da regressão linear ponderada
robusta com a função estimatr::lm_robust da seguinte
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], sep="")Between Pearson's correlation coefficient = 0.06975858, p = 0.9073215Regressão linear ponderada com banda de confiança:lm e
estimatr::lm_robust (HC2)
Nesta aplicação, lm ajusta uma regressão linear
ponderada por mínimos quadrados entre as médias por participante.
estimatr::lm_robust ajusta a mesma regressão ponderada,
mas calcula erros-padrão robustos à heterocedasticidade. O tipo HC2
aplica correção de alavancagem, melhorando a inferência em amostras
pequenas ou com pontos influentes. O coeficiente estimado é o mesmo do
lm; mudam os erros-padrão, IC e valores-p,
tornando-os mais robustos.
As bandas de confiança em torno da reta incluem inclinação nula ao longo de todo o intervalo observado de PaCO₂. Isso indica que a estimativa da inclinação é compatível com zero, portanto não significante para α = 5%. Em termos gráficos, a incerteza da regressão é grande o suficiente para não excluir ausência de efeito.
Esse resultado contraria a expectativa fisiológica de relação inversa entre PaCO₂ e pHi. A razão é estrutural: o ajuste usa dados médios entre participantes, removendo a variação intraparticipante onde o mecanismo ácido–base se manifesta, e introduz heterogeneidade basal entre indivíduos. Assim, a banda de confiança larga e centrada em inclinação estatisticamente nula reflete especificação inadequada para captar a relação fisiológica.
print(
ggplot2::ggplot(data = Dados.media,
mapping = ggplot2::aes(x = PaCO2.mean,
y = pH.mean)) +
ggplot2::geom_point() +
ggplot2::stat_smooth(method = lm,
formula = y ~ x,
se = TRUE,
mapping = ggplot2::aes(weight = Number)) +
ggplot2::labs(x = "PaCO2",
y = "ipH",
title = "Weighted Linear Regression") +
ggplot2::theme_test()
)
print(
ggplot2::ggplot(data = Dados.media,
mapping = ggplot2::aes(x = PaCO2.mean,
y = pH.mean)) +
ggplot2::geom_point() +
ggplot2::stat_smooth(method = estimatr::lm_robust,
formula = y ~ x,
se = TRUE,
mapping = ggplot2::aes(weight = Number)) +
ggplot2::labs(x = "PaCO2",
y = "ipH",
title = "Robust Weighted Linear Regression") +
ggplot2::theme_test()
)
Regressão linear robusta clusterizada:
estimatr::lm_robust (CR2)
Nesta aplicação, estimatr::lm_robust com CR2 ajusta uma
regressão linear simples pH ~ PaCO2 usando todos os pontos
e calcula erros-padrão robustos clusterizados por Subject. O estimador
de coeficiente é o mesmo de MQO; o CR2 corrige a variância para
dependência intra-cluster e pequeno número de clusters.
Apesar de tratar a correlação dentro do participante, o modelo estima uma única inclinação marginal comum sem decompor variação intra vs. interparticipantes. Assim, a heterogeneidade basal entre sujeitos e a mistura de níveis podem diluir o efeito fisiológico. O resultado mostra inclinação próxima de zero (\(\beta = −0.024\)), IC95% incluindo 0 e \(p = 0.87\). Logo, não significante para α = 5%.
Portanto, embora o CR2 melhore a inferência sob dependência, a especificação continua inadequada para recuperar a relação fisiológica intraparticipante PaCO₂–pHi; métodos que modelam explicitamente a estrutura longitudinal (rmcorr, modelos mistos com efeito aleatório de intercepto e/ou inclinação) são mais apropriados.
Call:
estimatr::lm_robust(formula = pH ~ PaCO2, data = Dados, clusters = Subject)
Standard error type: CR2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 7.28803 0.6926 10.5228 0.0003502 5.4030 9.1731 4.214
PaCO2 -0.02355 0.1306 -0.1803 0.8656233 -0.3848 0.3377 4.043
Multiple R-squared: 0.004253 , Adjusted R-squared: -0.01787
F-statistic: 0.03249 on 1 and 7 DF, p-value: 0.8621GEE: Generalized Estimating Equations
Equações de Estimativa Generalizadas (GEE) foram desenvolvidas para estender o Modelo Linear Generalizado (GLM) para acomodar dados longitudinais e agrupados.
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).
GEE: relação PaCO₂–pHi
O GEE estima um efeito marginal médio na população, não a relação fisiológica intraparticipante. Fisiologicamente, o mecanismo ácido–base atua dentro do mesmo indivíduo: variações de PaCO₂ ao longo do tempo produzem variações opostas de pHi no próprio paciente. O GEE com link identidade e estrutura “exchangeable” fornece uma inclinação média populacional, misturando variação intra e interparticipantes; essa média pode ser atenuada por heterogeneidade basal entre sujeitos (set-point, compensações, gravidade), diluindo o sinal fisiológico.
Além disso, a correlação de trabalho “exchangeable” impõe mesma correlação entre quaisquer pares dentro do sujeito, o que raramente reflete a dinâmica real; a má especificação reduz eficiência e pode ampliar IC, dificultando detectar a relação esperada. O resultado não significativo (\(\beta = −0.085, p = 0.12\)) não refuta a fisiologia; indica que o parâmetro marginal médio não capta adequadamente a associação intraparticipante PaCO₂–pHi.
Portanto, para recuperar a relação fisiológica dentro do indivíduo, modelos condicionais que separam níveis (rmcorr, modelos mistos com intercepto aleatório e inclinação comum) são mais apropriados do que o GEE marginal.
Pacotes R:
geepack::geeglm,glmtoolbox::glmgeeyags,MuMIn: QIC for GEE model selection
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fitgeeglm <- geepack::geeglm(pH ~ PaCO2,
id=Subject,
family=gaussian("identity"),
corstr="exchangeable",
data=Dados)
print(anova(fitgeeglm))Analysis of 'Wald statistic' Table
Model: gaussian, link: identity
Response: pH
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
PaCO2 1 2.4088 0.1207
Call:
geepack::geeglm(formula = pH ~ PaCO2, family = gaussian("identity"),
data = Dados, id = Subject, corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) 7.54416 0.32037 554.520 <2e-16 ***
PaCO2 -0.08512 0.05485 2.409 0.121
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 0.07436 0.04749
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.6513 0.147
Number of clusters: 8 Maximum cluster size: 9 GEE: Regressão logística: geepack::geeglm: artrite
reumatóide
Referência: Norusis, 2012
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.
Método. GEE (Generalized Estimating Equations) estende o GLM para dados correlacionados. Aqui usa-se família binomial com ligação logit e estrutura de correlação de trabalho exchangeable, assumindo a mesma correlação entre as três observações dentro de cada paciente. O parâmetro estimado é marginal (médio na população), não condicional ao indivíduo. A inferência usa variância robusta tipo sanduíche, consistente mesmo com possível erro na especificação da correlação.
Modelo ajustado.
\[ \log\left(\dfrac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 \,\text{Tratamento}_{\text{Droga ativa}} \]
sendo que \(\pi\) é a probabilidade de alívio satisfatório.
Resultados. \(\hat\beta_1=1.512\), erro-padrão \(0.523\), estatística de Wald \(=8.35\), \(p=3.9\times10^{-3}\). A razão de chances é \(\text{OR}=\exp(1.512)=4.54\). Intervalo de confiança aproximado de 95% para a OR: \(\exp(1.512 \pm 1.96\times0.523)= [1.63,\;12.63]\). O parâmetro de correlação intraclasse estimado é \(\hat\alpha=0.249\), indicando correlação positiva moderada entre semanas dentro do mesmo paciente.
Interpretação. No nível populacional, a Droga ativa aumenta significantemente as chances de alívio satisfatório em relação ao Placebo. Em média ao longo das três semanas, pacientes sob Droga ativa têm cerca de 4.5 vezes mais chances de relatar alívio satisfatório do que sob Placebo. O resultado é estatisticamente significativo para \(\alpha=5\%\) e robusto à dependência intrapaciente considerada pelo GEE.
OR por modelo logístico marginal (GEE)
No modelo logístico marginal (GEE),
\[ \log\left(\dfrac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 \,\text{Tratamento}_{\text{Droga ativa}} \]
sendo que \(\pi\) é a probabilidade de alívio satisfatório.
O coeficiente \(\beta_1\) representa a diferença de log-odds de alívio satisfatório entre Droga ativa e Placebo, no nível populacional (marginal).
Estimativa. \(\hat\beta_1=1.512\).
Conversão para razão de chances (odds ratio, OR):
\[ \text{OR}=\exp(\hat\beta_1)=\exp(1.512)=4.54. \]
Interpretação. OR compara as chances (odds) de alívio satisfatório sob Droga ativa versus Placebo, mantendo o restante do modelo constante e incorporando a correlação intrapaciente via GEE. Assim, \(\text{OR}=4.54\) significa que, em média ao longo das três semanas, as chances de relatar alívio satisfatório com Droga ativa são cerca de 4.5 vezes as chances sob Placebo.
Consistência inferencial. O teste de Wald para \(\beta_1\) fornece \(p=3.9\times10^{-3}\) e o IC95% aproximado para a OR é \([1.63,\;12.63]\), ambos excluindo 1, corroborando aumento significante das chances com a Droga ativa.
Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Sexo <- factor(Dados$Sexo)
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 Sexo
1 : 3 Min. :28.0 Placebo :68 Min. : 5.00 F: 33
2 : 3 1st Qu.:45.0 Droga ativa:67 1st Qu.: 5.00 M:102
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
AlivioDor Peso
Baixo :36 Min. :2
Satisfatório:99 1st Qu.:2
Median :2
Mean :2
3rd Qu.:2
Max. :2
# A tibble: 6 × 7
Paciente Idade Tratamento Semana Sexo AlivioDor Peso
<fct> <dbl> <fct> <dbl> <fct> <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> <fct> <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)
Dados$y <- ifelse(Dados$AlivioDor == "Satisfatório", 1, 0)
fitgeeglm <- geepack::geeglm(
y ~ Tratamento,
id = Paciente,
family = binomial("logit"),
corstr = "exchangeable",
data = Dados
)
print(summary(fitgeeglm))
Call:
geepack::geeglm(formula = y ~ Tratamento, family = binomial("logit"),
data = Dados, id = Paciente, corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) 0.372 0.317 1.38 0.2400
TratamentoDroga ativa 1.512 0.523 8.35 0.0039 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 0.962 0.264
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.249 0.129
Number of clusters: 50 Maximum cluster size: 3 Analysis of 'Wald statistic' Table
Model: binomial, link: logit
Response: y
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
Tratamento 1 8.35 0.0039 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## OR e IC94% (Wald) para Tratamento = Droga ativa vs Placebo
b <- coef(fitgeeglm)["TratamentoDroga ativa"]
se <- summary(fitgeeglm)$coefficients["TratamentoDroga ativa", "Std.err"]
alfa <- 0.05
zcrit <- stats::qnorm(1 - alfa/2) # IC95% => alfa = 0.05
OR <- exp(b)
lwr <- exp(b - zcrit * se)
upr <- exp(b + zcrit * se)
cat(sprintf("\nOR (Droga ativa vs Placebo) = %.3f\n", OR))
OR (Droga ativa vs Placebo) = 4.538IC95% da OR (Wald) = [1.627, 12.658]GEE: Regressão de Poisson: geepack::geeglm:
epilepsia
Referência: Norusis, 2012
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.
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 RR (risk ratio) do treatment para o nível drug é significante e menor que 1.
GEE: regressão de Poisson para contagem de crises epilépticas
Método. GEE estende o GLM para dados longitudinais correlacionados. Aqui usa-se família Poisson com ligação log e termo de offset \(\log(\text{time})\) para modelar taxa de crises por semana. A estrutura de correlação de trabalho trata a dependência entre medidas repetidas dentro do mesmo paciente. O parâmetro estimado é marginal (médio populacional). A inferência usa variância robusta tipo sanduíche.
Modelo. Para a taxa \(\lambda\) de crises por semana,
\[ \log(\lambda)=\beta_0+\beta_1\,\text{treatment}_{\text{Drug}}+\beta_2\,\text{active}_{\text{Treatment}}+\beta_3\,(\text{treatment}\times\text{active})+\\ \log(\text{time}) \]
sendo que \(\exp(\beta)\) fornece razão de riscos (risk ratio, RR).
Resultados (modelo Poisson com interação).
- Efeito principal de
active(período de tratamento vs baseline): \(\hat\beta_2=-1.274\), \(p<2\times10^{-16}\), indicando forte redução da taxa durante o tratamento em relação ao baseline.
- Interação
treatmentDrug:activeTreatment: \(\hat\beta_3=-0.302\), \(p=1.3\times10^{-5}\), indicando efeito adicional da Droga ativa durante o período de tratamento.
- Efeito principal de
treatmentno baseline não significativo (\(p\approx0.97\)), como esperado.
Interpretação clínica via razão de riscos. O efeito protetor da Droga ativa durante o tratamento é dado por
\[ \text{RR}=\exp(\hat\beta_3)=\exp(-0.302)=0.739 \]
com IC95% aproximado
\[ [\,\exp(-0.439),\ \exp(-0.166)\,]=[\,0.645,\ 0.847\,] \]
Como \(\text{RR}<1\) e o IC95% exclui 1, a Droga ativa reduz significantemente a taxa de episódios epilépticos durante o período de tratamento, em comparação ao Placebo, no nível populacional.
Leitura prática. RR de \(0.739\) implica redução média de cerca de \(26\%\) na taxa de crises sob Droga ativa durante o tratamento. A dependência intrapaciente foi considerada pelo modelo, e a conclusão é marginal (populacional), não específica a um indivíduo.
Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$sequence <- factor(Dados$sequence)
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 Baseline:58 Min. :2.0
2 : 5 1st Qu.:23 1st Qu.:3.14 Week1 :58 1st Qu.:2.0
3 : 5 Median :28 Median :3.33 Week2 :58 Median :2.0
4 : 5 Mean :29 Mean :3.33 Week3 :58 Mean :3.2
5 : 5 3rd Qu.:35 3rd Qu.:3.56 Week4 :58 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> <fct> <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> <fct> <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.739 [ 0.645 0.847 ]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 mista para pHi e PaCO₂
Método. Modelos lineares mistos (GLMM com família gaussiana, equivalente a LMM) acomodam medidas repetidas incluindo efeitos aleatórios por sujeito. Aqui:
\[ \text{pH}_{ij}=\beta_0+\beta_1\,\text{PaCO2}_{ij}+b_{0i}+\varepsilon_{ij} \]
sendo que \(b_{0i}\sim N(0,\sigma_b^2)\) é o intercepto aleatório do sujeito \(i\) (heterogeneidade basal) e \(\varepsilon_{ij}\sim N(0,\sigma^2)\) é o erro residual. O coeficiente \(\beta_1\) estima a associação intraparticipante comum (inclinação média), controlando diferenças entre sujeitos. Inferência pode usar Wald, F com Satterthwaite/Kenward–Roger, ou ML/REML; resultados devem ser coerentes.
Resultados principais (consistentes entre implementações):
- Inclinação fixa: \(\hat\beta_1=-0.1033\) (EP \(0.0294\)), \(t\approx-3.52\), \(p\approx1.1\times10^{-3}\)
- Testes globais: \(F\approx12.1\)–\(12.4\), \(p\approx(4.3\times10^{-4}\) a \(1.2\times10^{-3})\)
- Intervalo de confiança (nlme): \(\beta_1\in[-0.163,\,-0.044]\)
- Variâncias: \(\sigma_b^2\approx0.097\) (DP \(0.312\)) para interceptos entre sujeitos; \(\sigma^2\approx0.0088\) (DP \(0.094\)) residual.
Interpretação. A inclinação negativa e significante indica relação inversa intraparticipante: aumento de PaCO₂ associa-se a redução do pHi dentro do mesmo paciente. IC95% de \(\beta_1\) é inteiramente negativo e exclui 0; a associação é estatisticamente significante para \(\alpha=5\%\) e de magnitude clinicamente coerente com a fisiologia ácido–base.
Comparação entre métodos.
lme4::glmercom família gaussiana (atalho paralmer) elmerTest::lmerproduzem estimativas idênticas de \(\beta_1\); diferem apenas na forma de inferência (Wald vs F com Satterthwaite).
car::Anovacom ajuste Kenward–Roger fornece \(F\) e \(p\) praticamente iguais aos delmerTest.
nlme::lmeconfirma os mesmos coeficientes e IC95% similares, com pequena variação numérica por diferenças de algoritmo/parametrização.
Conclusão. Todos os modelos mistos concordam: há associação linear inversa significativa entre PaCO₂ e pHi ao nível intraparticipante, após controlar a heterogeneidade entre sujeitos.
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fitglmer <- lme4::glmer(pH ~ PaCO2 + (1|Subject),
family=gaussian("identity"),
data=Dados)Warning in lme4::glmer(pH ~ PaCO2 + (1 | Subject), family =
gaussian("identity"), : calling glmer() with family=gaussian (identity link) as
a shortcut to lmer() is deprecated; please call lmer() directlyLinear mixed model fit by REML ['lmerMod']
Formula: pH ~ PaCO2 + (1 | Subject)
Data: Dados
REML criterion at convergence: -50.4
Random effects:
Groups Name Std.Dev.
Subject (Intercept) 0.3116
Residual 0.0937
Number of obs: 47, groups: Subject, 8
Fixed Effects:
(Intercept) PaCO2
7.631 -0.103 Linear mixed model fit by REML ['lmerMod']
Formula: pH ~ PaCO2 + (1 | Subject)
Data: Dados
REML criterion at convergence: -50.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.2791 -0.3737 0.0371 0.4872 1.7557
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.09710 0.3116
Residual 0.00879 0.0937
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 7.6314 0.1848 41.30
PaCO2 -0.1033 0.0294 -3.52
Correlation of Fixed Effects:
(Intr)
PaCO2 -0.799Analysis of Variance Table
npar Sum Sq Mean Sq F value
PaCO2 1 0.109 0.109 12.4Analysis of Deviance Table (Type II Wald chisquare tests)
Response: pH
Chisq Df Pr(>Chisq)
PaCO2 12.4 1 0.00043 ***
---
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)
PaCO2 0.109 0.109 1 40.6 12.4 0.0011 **
---
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: pH ~ PaCO2 + (1 | Subject)
Data: Dados
REML criterion at convergence: -50.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.2791 -0.3737 0.0371 0.4872 1.7557
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.09710 0.3116
Residual 0.00879 0.0937
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 7.6314 0.1848 31.9013 41.30 <2e-16 ***
PaCO2 -0.1033 0.0294 40.5650 -3.52 0.0011 **
---
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: pH
F Df Df.res Pr(>F)
PaCO2 12.1 1 40.6 0.0012 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 numDF denDF F-value p-value
(Intercept) 1 38 4095 <.0001
PaCO2 1 38 12 0.0011Linear mixed-effects model fit by REML
Data: Dados
AIC BIC logLik
-42.4 -35.2 25.2
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 0.312 0.0937
Fixed effects: pH ~ PaCO2
Value Std.Error DF t-value p-value
(Intercept) 7.63 0.1848 38 41.3 0.0000
PaCO2 -0.10 0.0294 38 -3.5 0.0011
Correlation:
(Intr)
PaCO2 -0.799
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.2791 -0.3737 0.0371 0.4872 1.7557
Number of Observations: 47
Number of Groups: 8 Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 7.257 7.631 8.0054
PaCO2 -0.163 -0.103 -0.0438
Random Effects:
Level: Subject
lower est. upper
sd((Intercept)) 0.182 0.312 0.533
Within-group standard error:
lower est. upper
0.0749 0.0937 0.1174 Analysis of Deviance Table (Type II tests)
Response: pH
Chisq Df Pr(>Chisq)
PaCO2 12.4 1 0.00043 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Relação entre correlação intraparticipantes e inclinação de regressão por GLMM
Referência: Sobreiro, MFM, Silveira, PSP, Cavenaghi, VB, da Costa, LP, de Souza, B PF, Takahashi, RES, Starek, RVM, Siqueira, JO, & Fraguas, R (2025) Long-Term Cognitive Outcomes of Esketamine Nasal Spray in Treatment-Resistant Depression: A Preliminary Report. Pharmaceuticals, 18(2): 173. https://doi.org/10.3390/ph18020173
“Resumo: Contexto/Objetivos: A cetamina/esketamina apresenta efeito antidepressivo rápido e robusto na depressão resistente ao tratamento (DRT). Contudo, seus efeitos cognitivos de longo prazo permanecem incertos. Neste estudo, investigamos os possíveis efeitos cognitivos de um spray de esketamina em uma série de pacientes com DRT. Métodos: Avaliamos o desempenho cognitivo de oito pacientes com DRT submetidos a spray nasal de esketamina como tratamento adjuvante por seis meses. As avaliações cognitivas foram realizadas antes do início do tratamento (T0) e aos três (T3) e seis (T6) meses por um neuropsicólogo experiente, utilizando uma bateria neuropsicológica abrangente. A gravidade da depressão foi avaliada pela Escala de Avaliação da Depressão de Montgomery–Åsberg. As mudanças no desempenho cognitivo foram analisadas determinando-se o viés entre os tempos de avaliação. Para investigar a associação entre a gravidade da depressão e o desempenho nos testes cognitivos, utilizamos correlação com correção para medidas repetidas e análise de regressão com modelo linear misto geral (GLMM). Utilizamos o método de Tukey para comparar três estimativas e o método de Dunnett para comparar duas estimativas. Resultados: Melhoras em pelo menos um teste entre T0 e T6 foram observadas para atenção, memória e funções executivas de memória de trabalho, mudança de conjunto (set-shifting) e controle inibitório. A maioria das melhoras ocorreu até T3, mas as melhoras em memória de trabalho e mudança de conjunto foram significantes apenas em T6. A gravidade da depressão diminuiu significantemente de T0 a T6, e a maioria das melhoras cognitivas correlacionou-se com a redução da gravidade da depressão. Nenhum teste indicou piora do desempenho cognitivo entre T0 e T6. Conclusões: Nossos resultados sugerem que o desempenho cognitivo de pacientes com DRT melhorou com tratamento adjuvante de longo prazo com spray nasal de esketamina. Estudos confirmatórios são necessários.”
“Palavras-chave: cetamina; esketamina; neuropsicológico; cognitivo; longo prazo; depressão resistente ao tratamento; depressão”
“Duas estratégias foram empregadas para testar a associação entre os escores totais da MADRS e o escore de cada um dos testes neuropsicológicos. A primeira foi uma medida de correlação com correção para medidas repetidas [52], implementada em R utilizando a função rmcorr do pacote rmcorr [53]. A segunda estratégia envolveu uma análise de regressão obtida a partir de um modelo linear misto geral (GLMM), no qual o paciente foi utilizado como seu próprio controle (efeito aleatório) devido às medidas repetidas. Essa regressão foi implementada com a função
lmerdo pacotelmerTestdo R [54].”
## Dados (extraídos da tabela)
r <- c(0.079, 0.432, 0.694, 0.771, 0.140, -0.424, -0.526, -0.301, -0.599,
-0.192, -0.435, -0.788, -0.784, -0.024, -0.314, -0.574, -0.092,
0.221, 0.120, 0.315, -0.033, -0.657)
slope <- c(0.020, 0.086, 0.245, 0.519, 0.011, -0.212, -0.589, -0.053, -0.108,
-0.034, -0.039, -0.717, -0.610, 0.016, -0.085, -0.068, -0.025,
0.237, 0.076, 0.150, -0.009, -0.379)
## Correlação de Pearson
ct <- cor.test(r, slope, method = "pearson")
print(ct)
Pearson's product-moment correlation
data: r and slope
t = 8, df = 20, p-value = 2e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.694 0.942
sample estimates:
cor
0.863 ## Diagrama de dispersão + reta OLS
df <- data.frame(r = r, slope = slope)
print(
ggplot2::ggplot(df, ggplot2::aes(x = r, y = slope)) +
ggplot2::geom_point(size = 2) +
ggplot2::geom_smooth(method = "lm", formula = y ~ x, se = TRUE, color = "black") +
ggplot2::labs(x = "R.M. correlation (r)",
y = "GLMM slope",
title = sprintf("Pearson r = %.3f, p = %.3g", ct$estimate, ct$p.value)) +
ggplot2::theme_test()
)
GLMM: Regressão logística: lme4::glmer: artrite
reumatóide
Referência: Norusis, 2012
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.
Método. GLMM com família binomial e ligação logit acomoda medidas repetidas incluindo um intercepto aleatório por paciente. O modelo é
\[ \log\left(\dfrac{\pi_{ij}}{1-\pi_{ij}}\right)=\beta_0+\beta_1\,\text{Tratamento}_{\text{Droga ativa}}+b_{0i} \]
sendo que \(\pi_{ij}\) é a probabilidade de alívio satisfatório na observação \(j\) do paciente \(i\), e \(b_{0i}\sim N(0,\sigma_b^2)\) captura a heterogeneidade basal entre pacientes. A estimação é por máxima verossimilhança (aproximação de Laplace). O coeficiente \(\beta_1\) mede o efeito condicional ao paciente (within-subject). A razão de chances é \(\text{OR}=\exp(\beta_1)\).
Resultados. Teste de Wald para Tratamento:
\(\chi^2=7.7\), \(p=0.0055\). Estimativa fixa: \(\hat\beta_1=2.137\) (EP \(0.770\)), \(z=2.78\), \(p=0.0055\). A heterogeneidade entre
pacientes é substancial: \(\hat\sigma_b^2=2.43\) (DP \(1.56\)). OR estimada é
\[ \text{OR}=\exp(2.137)=8.47 \]
com IC95%
\[ [\,1.87,\;38.3\,] \]
Interpretação. Condicionalmente ao paciente, a Droga ativa aumenta significantemente as chances de alívio satisfatório em relação ao Placebo. Em média, para um mesmo paciente, as chances de alívio sob Droga ativa são cerca de \(8.5\) vezes as chances sob Placebo. Como o IC95% exclui \(1\) e \(p<0.05\), o efeito é estatisticamente significante.
Comparação conceitual. Diferente do GEE (efeito marginal populacional), o GLMM fornece efeito condicional ao indivíduo, explicitando a variabilidade entre pacientes via intercepto aleatório. Ambos podem indicar benefício do tratamento, mas OR do GLMM tende a ser maior (efeito condicional) quando há heterogeneidade entre sujeitos.
Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Sexo <- factor(Dados$Sexo)
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 Sexo
1 : 3 Min. :28.0 Placebo :68 Min. : 5.00 F: 33
2 : 3 1st Qu.:45.0 Droga ativa:67 1st Qu.: 5.00 M:102
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
AlivioDor Peso
Baixo :36 Min. :2
Satisfatório:99 1st Qu.:2
Median :2
Mean :2
3rd Qu.:2
Max. :2
# A tibble: 6 × 7
Paciente Idade Tratamento Semana Sexo AlivioDor Peso
<fct> <dbl> <fct> <dbl> <fct> <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> <fct> <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.47 [ 1.87 38.3 ]OR de Tratamento Placebo = 1GLMM: Regressão de Poisson: lme4::glmer: epilepsia
Referência: Norusis, 2012
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.
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
Método. GLMM de Poisson com ligação log modela taxas de contagem sob dependência intrapaciente via intercepto aleatório. Com offset \(\log(\text{time})\), o modelo é
\[ \log(\lambda_{ij})=\beta_0+\beta_1\,\text{treatment}_{\text{Drug}}+\beta_2\,\text{active}_{\text{Treatment}}+\beta_3\,(\text{treatment}\times\text{active})+\\b_{0i}+\log(\text{time}_{ij}) \]
sendo que \(\lambda_{ij}\) é a taxa de crises por semana do paciente \(i\) no período \(j\), \(b_{0i}\sim N(0,\sigma_b^2)\) captura heterogeneidade entre pacientes, e \(\exp(\beta)\) fornece razão de riscos (risk ratio, RR). A estimação é por máxima verossimilhança (aproximação de Laplace). Os efeitos são condicionais ao paciente.
Resultados.
- Efeito principal de
treatmentno baseline: \(\hat\beta_1=-0.008\), \(p=0.97\) (sem diferença no baseline).
- Efeito principal de
active(período de tratamento vs baseline): \(\hat\beta_2=-1.274\), \(p<2\times10^{-16}\) (forte redução geral da taxa durante tratamento).
- Interação
treatmentDrug:activeTreatment: \(\hat\beta_3=-0.302\), \(p=1.3\times10^{-5}\) (efeito adicional da Droga ativa durante tratamento).
- Heterogeneidade entre pacientes: \(\hat\sigma_b^2=0.514\) (DP \(0.717\)).
Interpretação via RR. O efeito protetor da Droga ativa durante o tratamento é
\[ \text{RR}=\exp(\hat\beta_3)=\exp(-0.302)=0.739 \]
com IC95%
\[ [\,\exp(-0.439),\ \exp(-0.166)\,]=[\,0.645,\ 0.847\,] \]
Como \(\text{RR}<1\) e o IC95% exclui \(1\), a Droga ativa reduz significantemente a taxa de crises no período de tratamento, comparada ao Placebo, condicionalmente ao paciente. RR de \(0.739\) implica redução média de cerca de \(26\%\) na taxa de crises sob Droga ativa durante o tratamento.
Leitura prática. Não há diferença entre grupos no baseline. Ao iniciar tratamento, a taxa cai em ambos, porém cai mais com Droga ativa (efeito de interação). O modelo controla a dependência intrapaciente e quantifica o efeito protetor como razão de taxas condicional.
Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$sequence <- factor(Dados$sequence)
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 Baseline:58 Min. :2.0
2 : 5 1st Qu.:23 1st Qu.:3.14 Week1 :58 1st Qu.:2.0
3 : 5 Median :28 Median :3.33 Week2 :58 Median :2.0
4 : 5 Mean :29 Mean :3.33 Week3 :58 Mean :3.2
5 : 5 3rd Qu.:35 3rd Qu.:3.56 Week4 :58 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> <fct> <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> <fct> <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.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 ]GEE vs. GLMM
Ideia central. Ambos analisam dados com medidas repetidas. A diferença é o tipo de pergunta clínica que cada um responde.
GEE (modelo marginal, média populacional).
Pergunta: em média, na população, o tratamento muda o
desfecho ao longo do tempo?
Interpretação: efeito médio do tratamento considerando
todos os pacientes juntos.
Propriedade: erros-padrão robustos; mesmo se a correlação
temporal for mal especificada, a inferência do efeito médio permanece
válida.
Leitura clínica: “no conjunto dos pacientes, o tratamento
melhora/piora o desfecho?”
GLMM (modelo misto, específico do paciente).
Pergunta: dentro de um mesmo paciente, o tratamento muda o
desfecho? Os pacientes diferem entre si na resposta?
Interpretação: efeito condicional ao paciente e
quantificação da variabilidade entre pacientes (heterogeneidade).
Vantagens: lida bem com dados faltantes e desbalanceados;
permite modelar diferenças individuais (ex.: alguns melhoram muito,
outros pouco).
Leitura clínica: “para um paciente típico, qual é o efeito do
tratamento? Os pacientes respondem de forma diferente?”
Relação entre resultados. Com correlação “exchangeable”, GEE e GLMM costumam dar coeficientes muito parecidos; os erros-padrão podem diferir um pouco. Em geral:
- GEE: efeito médio populacional.
- GLMM: efeito dentro do paciente + variabilidade entre pacientes.
Tradução prática.
Se a decisão clínica depende do efeito médio na população
(política/efetividade), use GEE.
Se interessa a resposta individual e a heterogeneidade
(prognóstico/medicina personalizada), use GLMM.
Um exemplo de semelhança dos resultados de GEE e GLMM está em Getting Started with Generalized Estimating Equations.
GEE vs GLMM: Regressão: pHi e PaCO2
Pergunta clínica. A relação fisiológica esperada é inversa: quando PaCO₂ aumenta, o pHi diminui dentro do mesmo paciente. Logo, o método deve captar a associação intraparticipante, controlando diferenças basais entre sujeitos.
GEE (modelo marginal, média populacional).
Modelo: \(pH=\beta_0+\beta_1\,\text{PaCO2}\) com
correlação “exchangeable” por sujeito.
Resultado: \(\hat\beta_1=-0.085\), erro-padrão \(0.055\), \(p=0.12\) (não significante).
Leitura: o GEE estima o efeito médio na
população. Com forte heterogeneidade entre sujeitos e poucos
clusters (\(n=8\)), o efeito médio pode
ser atenuado (mistura variação entre e dentro de sujeitos), reduzindo
poder para detectar a relação fisiológica.
GLMM/LMM (modelo misto, específico do paciente).
Modelo: \(pH_{ij}=\beta_0+\beta_1\,\text{PaCO2}_{ij}+b_{0i}+\varepsilon_{ij}\),
com intercepto aleatório \(b_{0i}\) por
sujeito.
Resultados consistentes entre lme4::glmer,
lmerTest::lmer:
\(\hat\beta_1=-0.1033\), \(t\approx-3.52\), \(F\approx12.1\)–\(12.4\), \(p\approx1.1\times10^{-3}\)–\(1.2\times10^{-3}\) (significante).
Observação. O ajuste via
lme4::glmer(..., family=gaussian("identity")) não é um
modelo diferente neste caso; é um atalho para lmer()
(deprecated). As pequenas diferenças de \(p\) decorrem do método usado para graus de
liberdade e teste (Wald, Satterthwaite, Kenward–Roger), não do
coeficiente estimado, que é o mesmo.
Leitura: o GLMM estima a associação intraparticipante
comum, removendo a variabilidade basal entre sujeitos (\(\sigma_b^2\approx0.097\)). O efeito
negativo é claro e coerente com a fisiologia ácido–base.
Resultados consistentes entre lme4::lmer e
lmerTest::lmer (e também nlme::lme):
\(\hat\beta_1=-0.1033\) (EP \(0.0294\)), \(t\approx-3.52\), \(F\approx12.1\)–\(12.4\), \(p\approx(1.1\times10^{-3}\) a \(1.2\times10^{-3})\) (significante).
Comparação e aderência fisiológica.
- GEE: efeito médio populacional, aqui não significante,
possivelmente diluído pela heterogeneidade entre sujeitos e pequeno
número de clusters.
- GLMM: efeito dentro do paciente, significante e negativo, alinhado com a fisiologia (PaCO₂↑ → pHi↓).
Conclusão. Para esta pergunta fisiológica (relação intraparticipante), o GLMM é mais apropriado e recupera a associação esperada; o GEE, por focar no efeito médio populacional, pode subestimar a relação quando há forte variabilidade entre sujeitos e poucos clusters.
Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fitgeeglm <- geepack::geeglm(pH ~ PaCO2,
id=Subject,
family=gaussian("identity"),
corstr="exchangeable",
data=Dados)
print(anova(fitgeeglm))Analysis of 'Wald statistic' Table
Model: gaussian, link: identity
Response: pH
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
PaCO2 1 2.41 0.12
Call:
geepack::geeglm(formula = pH ~ PaCO2, family = gaussian("identity"),
data = Dados, id = Subject, corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) 7.5442 0.3204 554.52 <2e-16 ***
PaCO2 -0.0851 0.0548 2.41 0.12
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 0.0744 0.0475
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.651 0.147
Number of clusters: 8 Maximum cluster size: 9 Warning in lme4::glmer(pH ~ PaCO2 + (1 | Subject), data = Dados): calling
glmer() with family=gaussian (identity link) as a shortcut to lmer() is
deprecated; please call lmer() directlyAnalysis of Variance Table
npar Sum Sq Mean Sq F value
PaCO2 1 0.109 0.109 12.4Linear mixed model fit by REML ['lmerMod']
Formula: pH ~ PaCO2 + (1 | Subject)
Data: Dados
REML criterion at convergence: -50.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.2791 -0.3737 0.0371 0.4872 1.7557
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.09710 0.3116
Residual 0.00879 0.0937
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error t value
(Intercept) 7.6314 0.1848 41.30
PaCO2 -0.1033 0.0294 -3.52Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: pH
F Df Df.res Pr(>F)
PaCO2 12.1 1 40.6 0.0012 **
---
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)
PaCO2 0.109 0.109 1 40.6 12.4 0.0011 **
---
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: pH ~ PaCO2 + (1 | Subject)
Data: Dados
REML criterion at convergence: -50.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.2791 -0.3737 0.0371 0.4872 1.7557
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 0.09710 0.3116
Residual 0.00879 0.0937
Number of obs: 47, groups: Subject, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 7.6314 0.1848 31.9013 41.30 <2e-16 ***
PaCO2 -0.1033 0.0294 40.5650 -3.52 0.0011 **
---
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: pH
F Df Df.res Pr(>F)
PaCO2 12.1 1 40.6 0.0012 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1GEE vs GLMM: Regressão logística: artrite reumatóide
Pergunta clínica. O tratamento ativo aumenta a probabilidade de alívio satisfatório? A diferença está em como o efeito é definido: média populacional (GEE) versus efeito dentro do paciente (GLMM).
Método GEE (modelo marginal, média populacional). O modelo é \[ \log\{\pi/(1-\pi)\}=\beta_0+\beta_1\,\text{Tratamento}_{\text{Droga ativa}}, \] com correlação intrapaciente “exchangeable” e variância robusta. O efeito é médio na população.
Resultados: \(\hat\beta_1=1.512\) (EP \(0.523\)), Wald \(=8.35\), \(p=3.9\times10^{-3}\). A razão de chances é \[ \text{OR}_{\text{GEE}}=\exp(1.512)=4.54,\quad \text{IC95%}=[1.63,\;12.66] \] Interpretação: em média na população, a Droga ativa multiplica por cerca de 4.5 as chances de alívio satisfatório em relação ao Placebo.
Método GLMM (modelo misto, específico do paciente). O modelo é \[ \log\{\pi_{ij}/(1-\pi_{ij})\}=\beta_0+\beta_1\,\text{Tratamento}_{\text{Droga ativa}}+b_{0i}\\ b_{0i}\sim N(0,\sigma_b^2) \] A estimação é por máxima verossimilhança (Laplace). O efeito é condicional ao paciente e a variabilidade entre pacientes é modelada explicitamente.
Resultados: \(\hat\beta_1=2.137\) (EP \(0.770\)), \(z=2.78\), \(p=5.5\times10^{-3}\), variância aleatória \(\hat\sigma_b^2=2.43\). A razão de chances é \[ \text{OR}_{\text{GLMM}}=\exp(2.137)=8.47,\quad \text{IC95%}=[1.87,\;38.3] \] Interpretação: para um mesmo paciente, a Droga ativa multiplica por cerca de 8.5 as chances de alívio satisfatório em relação ao Placebo.
Comparação. Ambos indicam benefício do tratamento. OR do GLMM é maior que a do GEE porque o efeito condicional (dentro do paciente) tende a exceder o efeito marginal (médio populacional) quando há heterogeneidade entre pacientes. GEE responde ao efeito médio na população; GLMM responde ao efeito dentro do paciente, controlando a linha de base individual. As duas abordagens são coerentes e complementares.
Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Sexo <- factor(Dados$Sexo)
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)
Dados$y <- ifelse(Dados$AlivioDor == "Satisfatório", 1, 0)
fitgeeglm <- geepack::geeglm(
y ~ Tratamento,
id = Paciente,
family = binomial("logit"),
corstr = "exchangeable",
data = Dados
)
print(summary(fitgeeglm))
Call:
geepack::geeglm(formula = y ~ Tratamento, family = binomial("logit"),
data = Dados, id = Paciente, corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) 0.372 0.317 1.38 0.2400
TratamentoDroga ativa 1.512 0.523 8.35 0.0039 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 0.962 0.264
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.249 0.129
Number of clusters: 50 Maximum cluster size: 3 Analysis of 'Wald statistic' Table
Model: binomial, link: logit
Response: y
Terms added sequentially (first to last)
Df X2 P(>|Chi|)
Tratamento 1 8.35 0.0039 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## OR e IC94% (Wald) para Tratamento = Droga ativa vs Placebo
b <- coef(fitgeeglm)["TratamentoDroga ativa"]
se <- summary(fitgeeglm)$coefficients["TratamentoDroga ativa", "Std.err"]
alfa <- 0.05
zcrit <- stats::qnorm(1 - alfa/2) # IC95% => alfa = 0.05
OR <- exp(b)
lwr <- exp(b - zcrit * se)
upr <- exp(b + zcrit * se)
cat(sprintf("\nOR (Droga ativa vs Placebo) = %.3f\n", OR))
OR (Droga ativa vs Placebo) = 4.538IC95% da OR (Wald) = [1.627, 12.658]GEE vs GLMM: Regressão de Poisson: epilepsia
Pergunta clínica. A droga ativa reduz a taxa de crises durante o período de tratamento? A comparação é entre efeito médio populacional (GEE) e efeito dentro do paciente (GLMM), ambos modelando taxa com ligação log e offset \(\log(\text{time})\).
Método GEE (marginal, média populacional). O modelo é \[ \log(\lambda)=\beta_0+\beta_1\,\text{treatment}_{\text{Drug}}+\beta_2\,\text{active}_{\text{Treatment}}+\beta_3\,(\text{treatment}\times\text{active})+\\\log(\text{time}) \] com correlação intrapaciente “exchangeable” e variância robusta.
Resultados: interação \(\hat\beta_3=-0.302\) (EP \(0.171\)), \(p=0.077\). A razão de riscos é \[ \text{RR}_{\text{GEE}}=\exp(-0.302)=0.739,\quad \text{IC95%}=[0.529,\;1.03] \] Interpretação: na média populacional há redução estimada de cerca de \(26\%\) na taxa sob droga ativa durante o tratamento, porém o IC95% inclui \(1\) e o resultado não é significante a \(\alpha=5\%\).
Método GLMM (misto, específico do paciente). O modelo é \[ \log(\lambda_{ij})=\beta_0+\beta_1\,\text{treatment}_{\text{Drug}}+\beta_2\,\text{active}_{\text{Treatment}}+\beta_3\,(\text{treatment}\times\text{active})+\\b_{0i}+\log(\text{time}_{ij})\\ b_{0i}\sim N(0,\sigma_b^2) \] A estimação é por máxima verossimilhança (Laplace).
Resultados: interação \(\hat\beta_3=-0.302\) (EP \(0.069\)), \(p=1.3\times10^{-5}\), heterogeneidade entre pacientes \(\hat\sigma_b^2=0.514\). A razão de riscos é \[ \text{RR}_{\text{GLMM}}=\exp(-0.302)=0.739,\quad \text{IC95%}=[0.645,\;0.847] \]
Interpretação: dentro do mesmo paciente, a droga ativa reduz significativamente a taxa de crises durante o tratamento em cerca de \(26\%\); o IC95% exclui \(1\).
Comparação. Os coeficientes são praticamente idênticos, mas a inferência difere. GEE estima o efeito médio populacional e, aqui, com maior variância e possível sobredispersão, produz IC mais largo e não significante. GLMM estima o efeito condicional ao paciente, modela a heterogeneidade entre sujeitos e fornece maior precisão, recuperando redução significante da taxa. Ambas as abordagens apontam redução semelhante em magnitude; a evidência estatística é mais forte no GLMM por capturar a variabilidade intrapaciente.
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)) # funciona as vezes...
print(summary(fitgeeglm, corr=FALSE))
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.
- Boyd, O; Mackay, CJ; Lamb, G; Bland, JM; Grounds, RM; & Bennett, ED (1993) Comparison of clinical information gained from routine blood-gas analysis and from gastric tonometry for intramural pH. Lancet 341: 142-6.
- 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.
- Sobreiro, MFM, Silveira, PSP, Cavenaghi, VB, da Costa, LP, de Souza, B PF, Takahashi, RES, Starek, RVM, Siqueira, JO, & Fraguas, R (2025) Long-Term Cognitive Outcomes of Esketamine Nasal Spray in Treatment-Resistant Depression: A Preliminary Report. Pharmaceuticals, 18(2) 173. https://doi.org/10.3390/ph18020173
- 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.