Bastão de Asclépio & Distribuição Normal

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

Adicionar

  • cor_multilevel

  • misty::multilevel.cor

  • Contextual 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. <br> Pekár & Brabec, 2017, p. 88.

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:

  1. 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.

  1. 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.

  2. 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:

  1. calcular medidas de resumo para cada cluster e analisar essas medidas de resumo usando métodos-padrão;
  2. usar erros-padrão robustos para corrigir erros-padrão para o agrupamento;
  3. utilizar modelos de efeitos aleatórios (GLMM) que modelam explicitamente a similaridade entre indivíduos de um mesmo cluster;
  4. 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

  1. 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.

  2. É sempre válido derivar medidas-resumo para cada agrupamento e depois analisá-las usando métodos-padrão.

  3. 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.

  4. 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.

  5. 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

Hamaker, 2013, in Mehl & Conner, p. 43-61.

Hamaker, 2013, in Mehl & Conner, p. 43-61.

Sainani, 2010, p. 861.

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

https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#/media/File:Correlation_examples2.svg

https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#/media/File:Correlation_examples2.svg

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.”

Baigorri et al., 1997: A tonometria gástrica determina PaCO<sub>2</sub> intraluminal, que se presume estar em equilíbrio com o PaCO<sub>2</sub> na mucosa gástrica. pHi pode ser calculado pela equação de Henderson-Hasselbalch, usando o valor de PaCO<sub>2</sub> determinado pela tonometria gástrica e a concentração de íon bicarbonato no sangue arterial.

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).

Análise descritiva de correlação intraparticipantes

alfa <- 0.05

Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
print(Dados)
   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
## rmcorr
fit <- rmcorr::rmcorr(participant = Subject,
                      measure1 = PaCO2,
                      measure2 = pH,
                      dataset = Dados)

r_rm <- fit$r
p_rm <- fit$p

lab_rm <- sprintf("r_rm = %.3f, p = %.3g", r_rm, p_rm)

## car
car::scatterplot(Dados$PaCO2,
                 Dados$pH,
                 groups = Dados$Subject,
                 regLine = FALSE,
                 smooth = FALSE,
                 ellipse = list(levels = c(.68),
                                robust = TRUE,
                                fill = FALSE),
                 col = "black",
                 grid = FALSE,
                 xlab = "PaCO2",
                 ylab = "pHi")
Warning in scatterplot.default(Dados$PaCO2, Dados$pH, groups = Dados$Subject, : number of groups exceeds number of available colors
  colors are recycled

## posição abaixo de todos
y_pos <- min(Dados$pH, na.rm = TRUE) - 0.05 * diff(range(Dados$pH, na.rm = TRUE))

## ggplot único
print(
  ggplot2::ggplot(data = Dados,
                  mapping = ggplot2::aes(x = PaCO2,
                                         y = pH,
                                         color = Subject)) +
    ggplot2::geom_point() +
    ggplot2::geom_smooth(method = "lm",
                         formula = y ~ x,
                         se = FALSE) +
    ggpubr::stat_cor(method = "pearson",
                     cor.coef.name = "r",
                     na.rm = TRUE,
                     r.digits = 3,
                     p.digits = 2,
                     label.x = 6.5) +
    ggplot2::annotate("text",
                      x = 6.5,
                      y = y_pos,
                      label = lab_rm,
                      hjust = 0,
                      vjust = 1,
                      color = "black") +
    ggplot2::labs(x = "PaCO2", y = "pHi") +
    ggplot2::xlim(3, 9) +
    ggplot2::theme_test()
)

Correlação intraparticipantes por cor.test: incorreta

print(cor.test(Dados$PaCO2, Dados$pH))

    Pearson's product-moment correlation

data:  Dados$PaCO2 and Dados$pH
t = -0.43843, df = 45, p-value = 0.6632
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3459063  0.2261851
sample estimates:
        cor 
-0.06521774

Para um delineamento intraparticipantes, é possível estimar e testar a o coeficiente de correlação de Pearson implícito entre duas variáveis intervalares. Esta correlação implícita intraparticipantes nao é: (i) a média aritmética das correlações dos participantes; (ii) a média ponderada das correlações dos participantes. Bland & Altman (1995a) mostraram a maneira correta de estimar a correlação intra.

Bland & Altman, 1995a

Bland & Altman, 1995a

IC95% por bootstrapping das correlações intraparticipantes

alfa <- 0.05

## 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.88 -0.053 0.960
2 1.00 1.000 1.000
3 -0.95 -0.610 0.691
4 -0.20 0.490 1.000
5 -0.93 -0.690 0.641
6 -1.00 -0.430 0.271
8 -0.97 -0.840 0.179

Correlação intraparticipantes por bootstrapping: correta

A estimativa do valor absoluto da correlação de Pearson para medidas repetidas é a raiz quadrada do eta parcial ao quadrado de PaCO2.

O sinal da correlação é obtido por meio do sinal da estimativa do parâmetro beta da regressão entre PaCO2 e pH intramural controlada pelo fator Subject.

# Bland & Altman, 1995a
print(rmcorr::rmcorr(participant=Subject, 
                     measure1=PaCO2, 
                     measure2=pH, 
                     dataset=Dados,
                     CI.level=1-alfa,
                     CIs="bootstrap",
                     nreps=1e4), digits=4)

Repeated measures correlation

r
-0.5067697

degrees of freedom
38

p-value
0.0008471081

95% confidence interval
-0.7287354 -0.06696856 
# 95% confidence interval (analytic)  = -0.7112297 -0.223255 
# 95% confidence interval (bootstrap) = -0.727479  -0.06999879

Medida média: “entre participantes”

Dados.media <- aggregate(Dados, x=cbind(PaCO2, pH)~Subject, FUN=mean)
Dados.n <- aggregate(Dados, x=cbind(Subject)~Subject, FUN=length)
Dados.media <- cbind(Dados.media,Dados.n[,2])
names(Dados.media) <- c(names(Dados.media)[1:(length(names(Dados.media))-1)], 
                        "Number")
names(Dados.media) <- c("Subject", "PaCO2.mean", "pH.mean", "Number")
print(Dados.media, digits=4)
  Subject PaCO2.mean pH.mean Number
1       1      4.037   6.492      4
2       2      5.373   7.053      4
3       3      4.830   7.357      9
4       4      5.312   7.326      5
5       5      4.399   7.314      8
6       6      4.920   7.323      6
7       7      6.603   6.907      3
8       8      4.784   7.115      8

Correlação “entre participantes” por cor.test: incorreta

print(cor.test(Dados.media$PaCO2.mean, Dados.media$pH.mean))

    Pearson's product-moment correlation

data:  Dados.media$PaCO2.mean and Dados.media$pH.mean
t = 0.20098, df = 6, p-value = 0.8474
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6609881  0.7435975
sample estimates:
       cor 
0.08177314

Correlação “entre participantes” por weights::wtd.cor: correta

# Bland & Altman 1995b: analitica e boostrapping
print(weights::wtd.cor(Dados.media$PaCO2.mean, 
                       Dados.media$pH.mean, 
                       weight=Dados.media$Number,
                       bootp=TRUE, bootse=TRUE, bootn=1e4), digits=4)
  correlation  bootcor std.err t.value p.value
Y     0.06976 0.004722  0.5837 0.00809  0.8854

Desenhando 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.

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, ?,

“estimatr is a package in R dedicated to providing fast estimators that take into consideration designs often used by social scientists. Estimators are statistical methods for estimating quantities of interest like treatment effects or regression parameters. Many of the estimators included with the R programming language or popular R packages are slow and have default settings that lead to statistically inappropriate estimates. Certain estimators that reflect cutting-edge advances in statistics are not yet implemented in R packages for convenient use. estimatr is designed to solve these problems and provide estimators tuned for design-based inference. The current estimators we provide are:

  • lm_robust: for fitting linear models with heteroskedasticity/cluster-robust standard errors
  • lm_lin: a wrapper for lm_robust() to simplify interacting centered pre-treatment covariates with a treatment variable
  • iv_robust: two stage least squares estimation of instrumental variables regression
  • difference_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 designs
  • horvitz_thompson: for estimating average treatment effects taking into consideration treatment probabilities or sampling probabilities for simple and cluster randomized designs”

Regressão linear robusta INCORRETA: estimatr::lm_robust (HC2)

fit <- estimatr::lm_robust(pH ~ PaCO2, 
                           data=Dados)
print(summary(fit))

Call:
estimatr::lm_robust(formula = pH ~ PaCO2, data = Dados)

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
(Intercept)  7.28803    0.33926   21.48 2.733e-25   6.6047   7.9713 45
PaCO2       -0.02355    0.06541   -0.36 7.205e-01  -0.1553   0.1082 45

Multiple R-squared:  0.004253 , Adjusted R-squared:  -0.01787 
F-statistic: 0.1296 on 1 and 45 DF,  p-value: 0.7205

Regressão linear robusta ponderada (com dados médios): estimatr::lm_robust (HC2)

A correlação de Pearson entre participantes entre os valores médios de pHi e PaCO2 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.9073
sunflowerplot(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.9073215

Regressão linear ponderada com banda de confiança:lm e estimatr::lm_robust (HC2)

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)

fit <- estimatr::lm_robust(pH ~ PaCO2, 
                           clusters=Subject, 
                           data=Dados)
print(summary(fit))

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.8621

GEE: 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).

Pacotes do R

  • geepack::geeglm, glmtoolbox::glmgee

  • yags, MuMIn: QIC for GEE model selection

Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)

fitgeeglm <- geepack::geeglm(PaCO2 ~ pH,
                             id=Subject,
                             family=gaussian("identity"),
                             corstr="exchangeable",
                             data=Dados)
print(anova(fitgeeglm))
Analysis of 'Wald statistic' Table
Model: gaussian, link: identity
Response: PaCO2
Terms added sequentially (first to last)

   Df      X2 P(>|Chi|)
pH  1 0.33044    0.5654
print(summary(fitgeeglm, corr=FALSE))

Call:
geepack::geeglm(formula = PaCO2 ~ pH, family = gaussian("identity"), 
    data = Dados, id = Subject, corstr = "exchangeable")

 Coefficients:
            Estimate Std.err  Wald Pr(>|W|)
(Intercept)   9.6493  8.2509 1.368    0.242
pH           -0.6538  1.1374 0.330    0.565

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)   0.5574   0.174
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha   0.3205   0.158
Number of clusters:   8  Maximum cluster size: 9

GEE: Regressão logística: geepack::geeglm: artrite reumatóide

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.

Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Tratamento <- factor(Dados$Tratamento)
Dados$Tratamento <- relevel(Dados$Tratamento,
                            ref="Placebo")
Dados$AlivioDor <- factor(Dados$AlivioDor)
Dados$AlivioDor <- relevel(Dados$AlivioDor,
                           ref="Baixo")
Dados$Semana <- as.numeric(Dados$Semana)
print(summary(Dados))
    Paciente       Idade            Tratamento     Semana     
 1      :  3   Min.   :28.0   Placebo    :68   Min.   : 5.00  
 2      :  3   1st Qu.:45.0   Droga ativa:67   1st Qu.: 5.00  
 3      :  3   Median :55.0                    Median : 9.00  
 4      :  3   Mean   :50.9                    Mean   : 8.82  
 5      :  3   3rd Qu.:60.0                    3rd Qu.:13.00  
 6      :  3   Max.   :65.0                    Max.   :13.00  
 (Other):117                                                  
     Sexo                  AlivioDor       Peso  
 Length:135         Baixo       :36   Min.   :2  
 Class :character   Satisfatório:99   1st Qu.:2  
 Mode  :character                     Median :2  
                                      Mean   :2  
                                      3rd Qu.:2  
                                      Max.   :2  
                                                 
print(head(Dados))
# A tibble: 6 × 7
  Paciente Idade Tratamento  Semana Sexo  AlivioDor     Peso
  <fct>    <dbl> <fct>        <dbl> <chr> <fct>        <dbl>
1 1           48 Droga ativa      5 M     Satisfatório     2
2 1           48 Droga ativa      9 M     Satisfatório     2
3 1           48 Droga ativa     13 M     Satisfatório     2
4 2           29 Droga ativa      5 M     Satisfatório     2
5 2           29 Droga ativa      9 M     Satisfatório     2
6 2           29 Droga ativa     13 M     Satisfatório     2
print(tail(Dados))
# A tibble: 6 × 7
  Paciente Idade Tratamento  Semana Sexo  AlivioDor     Peso
  <fct>    <dbl> <fct>        <dbl> <chr> <fct>        <dbl>
1 50          58 Droga ativa      5 M     Satisfatório     2
2 50          58 Droga ativa      9 M     Satisfatório     2
3 50          58 Droga ativa     13 M     Satisfatório     2
4 51          37 Placebo          5 M     Baixo            2
5 51          37 Placebo          9 M     Satisfatório     2
6 51          37 Placebo         13 M     Baixo            2
Dados <- na.omit(Dados)
fitgeeglm <- try(geepack::geeglm(AlivioDor ~ Tratamento, 
                                 id=Paciente,
                                 family=binomial("logit"),
                                 corstr="exchangeable",
                                 data=Dados))
Warning in model.response(mf, "numeric"): uso de type="numeric" com uma
resposta fator será ignorado
Warning in Ops.factor(y, mu): '-' não faz sentido para fatores
Error in lm.fit(zsca, qlf(pr2), offset = soffset) : NA/NaN/Inf in 'y'

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 o RR (risk ratio) do treatment para o nível drug é significante e menor que 1.

Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$treatment <- factor(Dados$treatment)
Dados$treatment <- relevel(Dados$treatment, ref="Placebo")
Dados$active <- factor(Dados$active)
Dados$active <- relevel(Dados$active, ref="Baseline")
print(summary(Dados))
       ID           age         agelog       sequence              time    
 1      :  5   Min.   :18   Min.   :2.89   Length:290         Min.   :2.0  
 2      :  5   1st Qu.:23   1st Qu.:3.14   Class :character   1st Qu.:2.0  
 3      :  5   Median :28   Median :3.33   Mode  :character   Median :2.0  
 4      :  5   Mean   :29   Mean   :3.33                      Mean   :3.2  
 5      :  5   3rd Qu.:35   3rd Qu.:3.56                      3rd Qu.:2.0  
 6      :  5   Max.   :57   Max.   :4.04                      Max.   :8.0  
 (Other):260                                                               
    logtime            active      treatment      seizure          visit4   
 Min.   :0.693   Baseline : 58   Placebo:140   Min.   :  0.0   Min.   :0.0  
 1st Qu.:0.693   Treatment:232   Drug   :150   1st Qu.:  3.0   1st Qu.:0.0  
 Median :0.693                                 Median :  6.0   Median :0.0  
 Mean   :0.970                                 Mean   : 11.5   Mean   :0.2  
 3rd Qu.:0.693                                 3rd Qu.: 14.0   3rd Qu.:0.0  
 Max.   :2.079                                 Max.   :111.0   Max.   :1.0  
                                                                            
    baseline        cseizure      logcseizure   
 Min.   :  6.0   Min.   : 0.00   Min.   :0.000  
 1st Qu.: 12.0   1st Qu.: 1.50   1st Qu.:0.916  
 Median : 22.0   Median : 2.38   Median :1.216  
 Mean   : 29.2   Mean   : 3.57   Mean   :1.264  
 3rd Qu.: 41.0   3rd Qu.: 4.50   3rd Qu.:1.705  
 Max.   :111.0   Max.   :38.00   Max.   :3.664  
                                                
print(head(Dados))
# A tibble: 6 × 13
  ID      age agelog sequence  time logtime active    treatment seizure visit4
  <fct> <dbl>  <dbl> <chr>    <dbl>   <dbl> <fct>     <fct>       <dbl>  <dbl>
1 1        31   3.43 Baseline     8   2.08  Baseline  Placebo        11      0
2 1        31   3.43 Week1        2   0.693 Treatment Placebo         5      0
3 1        31   3.43 Week2        2   0.693 Treatment Placebo         3      0
4 1        31   3.43 Week3        2   0.693 Treatment Placebo         3      0
5 1        31   3.43 Week4        2   0.693 Treatment Placebo         3      1
6 2        30   3.40 Baseline     8   2.08  Baseline  Placebo        11      0
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl>
print(tail(Dados))
# A tibble: 6 × 13
  ID      age agelog sequence  time logtime active    treatment seizure visit4
  <fct> <dbl>  <dbl> <chr>    <dbl>   <dbl> <fct>     <fct>       <dbl>  <dbl>
1 58       36   3.58 Week4        2   0.693 Treatment Drug            0      1
2 59       37   3.61 Baseline     8   2.08  Baseline  Drug           12      0
3 59       37   3.61 Week1        2   0.693 Treatment Drug            1      0
4 59       37   3.61 Week2        2   0.693 Treatment Drug            4      0
5 59       37   3.61 Week3        2   0.693 Treatment Drug            3      0
6 59       37   3.61 Week4        2   0.693 Treatment Drug            2      1
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl>
fitglmer <- lme4::glmer(seizure ~ treatment*active + (1|ID), 
                        data=Dados, 
                        family=poisson("log"))
print(car::Anova(fitglmer))
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: seizure
                   Chisq Df Pr(>Chisq)    
treatment           0.64  1       0.42    
active           1667.49  1    < 2e-16 ***
treatment:active   18.97  1    1.3e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
out <- summary(fitglmer, corr=FALSE)
print(out)
Generalized 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 ' ' 1
ci <- confint(fitglmer)
Computing profile confidence intervals ...
print(ci, digits=4)
                                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.1659
cat("RR de treatment(drug) =", exp(out$coefficients[4]), 
    "[", exp(as.numeric(ci[5,1:2])), "]")
RR de treatment(drug) = 0.7391 [ 0.6446 0.8472 ]

GLMM: Generalized Linear Mixed Model

GLMM é o modelo linear generalizado de efeitos fixos e aleatórios para uma VD.

Fator fixo:

  • é a população dos níveis do fator;
  • contém todos os níveis de interesse do pesquisador;
  • efeito de interesse do pesquisador a ser testado estatisticamente.

Fator aleatório:

  • amostra aleatória dos níveis do fator;
  • contém os níveis de interesse do pesquisador com dependência entre os valores da VD dentro de cada nível e independência entre os níveis;
  • efeito de controle usado para testar estatisticamente efeito fixo.

Na presença de fator aleatório (agrupamento), a conclusão sobre o efeito do fator fixo se aplica à população do fator aleatório.

O efeito de interação entre efeitos fixo e aleatório é aleatório.

Todo modelo com pelo menos um fator aleatório é misto, pois o intercepto é efeito fixo.

GLMM não exige independência entre as observações e homocedasticidade da VD.

Causas de dependência:

  • Agrupamento (cluster): região, escola, hospital, loja, participante;

  • Medidas repetidas: mensurações periódicas das mesmas variáveis nas mesmas unidades experimentais.

Distribuição da VD (FD):

  • A VD não precisa ser necessariamente normal. GLMM permite VD com distribuição da família exponencial: normal ou gaussiana, binomial (sucesso/fracasso), Poisson (eventos raros: acidente, incidência de doença), gama (valores positivos com assimetria positiva: duração de vida), gaussiana inversa (valores positivos com assimetria positiva: duração de casamento).

Função de ligação (FL):

A FL é função não-linear que estabelece uma relação linear entre a média condicionada da VD e a combinação linear das VI.

A combinação entre a distribuição da VD e a função de ligação resulta em diferentes modelos:

  • GLM (ANOVA e Regressão): FD: normal, FL: identidade
  • Modelo de regressão logística binária: FD: binomial, FL: logit
  • Modelo loglinear: FD: Poisson, FL: log
  • Modelo probit: FD: binomial, FL: probit
  • Modelo de regressão ordinal: FD: multinomial, FL: logit

Modelos hierárquicos

GLMM permite a integração os dados de níveis hierárquicos na mesma análise.

Níveis de dependência por agrupamento

  • Nível 1: dados sobre a unidade observacional
  • Nível 2: dados sobre o agrupamento 1
  • Nível 3: …

GLMM: Regressão linear: lme4::glmer: pHi e PaCO2

Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)
fitglmer <- lme4::glmer(PaCO2 ~ pH + (1|Subject),
                        family=gaussian("identity"),
                        data=Dados)
Warning in lme4::glmer(PaCO2 ~ pH + (1 | Subject), family =
gaussian("identity"), : calling glmer() with family=gaussian (identity link) as
a shortcut to lmer() is deprecated; please call lmer() directly
print(fitglmer)
Linear mixed model fit by REML ['lmerMod']
Formula: PaCO2 ~ pH + (1 | Subject)
   Data: Dados
REML criterion at convergence: 80.5
Random effects:
 Groups   Name        Std.Dev.
 Subject  (Intercept) 0.959   
 Residual             0.442   
Number of obs: 47, groups:  Subject, 8
Fixed Effects:
(Intercept)           pH  
      17.87        -1.81  
print(summary(fitglmer))
Linear mixed model fit by REML ['lmerMod']
Formula: PaCO2 ~ pH + (1 | Subject)
   Data: Dados

REML criterion at convergence: 80.5

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.430 -0.389  0.071  0.518  2.090 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.919    0.959   
 Residual             0.195    0.442   
Number of obs: 47, groups:  Subject, 8

Fixed effects:
            Estimate Std. Error t value
(Intercept)   17.867      4.161    4.29
pH            -1.805      0.583   -3.10

Correlation of Fixed Effects:
   (Intr)
pH -0.997
print(anova(fitglmer))
Analysis of Variance Table
   npar Sum Sq Mean Sq F value
pH    1   1.87    1.87    9.59
print(car::Anova(fitglmer))
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: PaCO2
   Chisq Df Pr(>Chisq)   
pH  9.59  1      0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fitlmer <- lmerTest::lmer(PaCO2 ~ pH + (1|Subject),
                          data=Dados)
print(anova(fitlmer))
Type III Analysis of Variance Table with Satterthwaite's method
   Sum Sq Mean Sq NumDF DenDF F value Pr(>F)   
pH   1.87    1.87     1  43.1    9.59 0.0034 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(fitlmer, corr=FALSE))
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: PaCO2 ~ pH + (1 | Subject)
   Data: Dados

REML criterion at convergence: 80.5

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.430 -0.389  0.071  0.518  2.090 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.919    0.959   
 Residual             0.195    0.442   
Number of obs: 47, groups:  Subject, 8

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   17.867      4.161 42.857    4.29  9.9e-05 ***
pH            -1.805      0.583 43.116   -3.10   0.0034 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(car::Anova(fitglmer, test.statistic="F"))
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: PaCO2
      F Df Df.res Pr(>F)   
pH 8.56  1   43.4 0.0055 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fitlme <- nlme::lme(PaCO2 ~ pH,
                    random=~1|Subject,
                    data=Dados)
print(anova(fitlme))
            numDF denDF F-value p-value
(Intercept)     1    38  211.12  <.0001
pH              1    38    9.59  0.0037
print(summary(fitlme, corr=FALSE))
Linear mixed-effects model fit by REML
  Data: Dados 
   AIC   BIC logLik
  88.5 95.73 -40.25

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:      0.9587    0.442

Fixed effects:  PaCO2 ~ pH 
             Value Std.Error DF t-value p-value
(Intercept) 17.867     4.161 38   4.294  0.0001
pH          -1.805     0.583 38  -3.097  0.0037
 Correlation: 
   (Intr)
pH -0.997

Standardized Within-Group Residuals:
     Min       Q1      Med       Q3      Max 
-2.43015 -0.38938  0.07101  0.51771  2.09026 

Number of Observations: 47
Number of Groups: 8 
print(nlme::intervals(fitlme))
Approximate 95% confidence intervals

 Fixed effects:
             lower   est.   upper
(Intercept)  9.444 17.867 26.2902
pH          -2.985 -1.805 -0.6253

 Random Effects:
  Level: Subject 
                 lower   est. upper
sd((Intercept)) 0.5228 0.9587 1.758

 Within-group standard error:
 lower   est.  upper 
0.3518 0.4420 0.5553 
print(car::Anova(fitlme, test.statistic="F"))
Analysis of Deviance Table (Type II tests)

Response: PaCO2
   Chisq Df Pr(>Chisq)   
pH  9.59  1      0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Relação entre correlação intraparticipantes e inclinação de 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

## 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 = 7.6, df = 20, p-value = 2e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6937 0.9419
sample estimates:
   cor 
0.8629 
## 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

Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Tratamento <- factor(Dados$Tratamento)
Dados$Tratamento <- relevel(Dados$Tratamento,
                            ref="Placebo")
Dados$AlivioDor <- factor(Dados$AlivioDor)
Dados$AlivioDor <- relevel(Dados$AlivioDor,
                           ref="Baixo")
Dados$Semana <- as.numeric(Dados$Semana)
print(summary(Dados))
    Paciente       Idade            Tratamento     Semana     
 1      :  3   Min.   :28.0   Placebo    :68   Min.   : 5.00  
 2      :  3   1st Qu.:45.0   Droga ativa:67   1st Qu.: 5.00  
 3      :  3   Median :55.0                    Median : 9.00  
 4      :  3   Mean   :50.9                    Mean   : 8.82  
 5      :  3   3rd Qu.:60.0                    3rd Qu.:13.00  
 6      :  3   Max.   :65.0                    Max.   :13.00  
 (Other):117                                                  
     Sexo                  AlivioDor       Peso  
 Length:135         Baixo       :36   Min.   :2  
 Class :character   Satisfatório:99   1st Qu.:2  
 Mode  :character                     Median :2  
                                      Mean   :2  
                                      3rd Qu.:2  
                                      Max.   :2  
                                                 
print(head(Dados))
# A tibble: 6 × 7
  Paciente Idade Tratamento  Semana Sexo  AlivioDor     Peso
  <fct>    <dbl> <fct>        <dbl> <chr> <fct>        <dbl>
1 1           48 Droga ativa      5 M     Satisfatório     2
2 1           48 Droga ativa      9 M     Satisfatório     2
3 1           48 Droga ativa     13 M     Satisfatório     2
4 2           29 Droga ativa      5 M     Satisfatório     2
5 2           29 Droga ativa      9 M     Satisfatório     2
6 2           29 Droga ativa     13 M     Satisfatório     2
print(tail(Dados))
# A tibble: 6 × 7
  Paciente Idade Tratamento  Semana Sexo  AlivioDor     Peso
  <fct>    <dbl> <fct>        <dbl> <chr> <fct>        <dbl>
1 50          58 Droga ativa      5 M     Satisfatório     2
2 50          58 Droga ativa      9 M     Satisfatório     2
3 50          58 Droga ativa     13 M     Satisfatório     2
4 51          37 Placebo          5 M     Baixo            2
5 51          37 Placebo          9 M     Satisfatório     2
6 51          37 Placebo         13 M     Baixo            2
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 ' ' 1
print(summary(fit, corr=FALSE))
Generalized 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 ' ' 1
or <- broom.mixed::tidy(fit,
                        conf.int=TRUE,
                        exponentiate=TRUE,
                        effects="fixed")
cat("OR de Tratamento Droga ativa =", or$estimate[2], 
    "[", or$conf.low[2], or$conf.high[2], "]")
OR de Tratamento Droga ativa = 8.474 [ 1.873 38.33 ]
cat("OR de Tratamento Placebo =", 1)
OR de Tratamento Placebo = 1

GLMM: Regressão de Poisson: lme4::glmer: epilepsia

Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$treatment <- factor(Dados$treatment)
Dados$treatment <- relevel(Dados$treatment, ref="Placebo")
Dados$active <- factor(Dados$active)
Dados$active <- relevel(Dados$active, ref="Baseline")
print(summary(Dados))
       ID           age         agelog       sequence              time    
 1      :  5   Min.   :18   Min.   :2.89   Length:290         Min.   :2.0  
 2      :  5   1st Qu.:23   1st Qu.:3.14   Class :character   1st Qu.:2.0  
 3      :  5   Median :28   Median :3.33   Mode  :character   Median :2.0  
 4      :  5   Mean   :29   Mean   :3.33                      Mean   :3.2  
 5      :  5   3rd Qu.:35   3rd Qu.:3.56                      3rd Qu.:2.0  
 6      :  5   Max.   :57   Max.   :4.04                      Max.   :8.0  
 (Other):260                                                               
    logtime            active      treatment      seizure          visit4   
 Min.   :0.693   Baseline : 58   Placebo:140   Min.   :  0.0   Min.   :0.0  
 1st Qu.:0.693   Treatment:232   Drug   :150   1st Qu.:  3.0   1st Qu.:0.0  
 Median :0.693                                 Median :  6.0   Median :0.0  
 Mean   :0.970                                 Mean   : 11.5   Mean   :0.2  
 3rd Qu.:0.693                                 3rd Qu.: 14.0   3rd Qu.:0.0  
 Max.   :2.079                                 Max.   :111.0   Max.   :1.0  
                                                                            
    baseline        cseizure      logcseizure   
 Min.   :  6.0   Min.   : 0.00   Min.   :0.000  
 1st Qu.: 12.0   1st Qu.: 1.50   1st Qu.:0.916  
 Median : 22.0   Median : 2.38   Median :1.216  
 Mean   : 29.2   Mean   : 3.57   Mean   :1.264  
 3rd Qu.: 41.0   3rd Qu.: 4.50   3rd Qu.:1.705  
 Max.   :111.0   Max.   :38.00   Max.   :3.664  
                                                
print(head(Dados))
# A tibble: 6 × 13
  ID      age agelog sequence  time logtime active    treatment seizure visit4
  <fct> <dbl>  <dbl> <chr>    <dbl>   <dbl> <fct>     <fct>       <dbl>  <dbl>
1 1        31   3.43 Baseline     8   2.08  Baseline  Placebo        11      0
2 1        31   3.43 Week1        2   0.693 Treatment Placebo         5      0
3 1        31   3.43 Week2        2   0.693 Treatment Placebo         3      0
4 1        31   3.43 Week3        2   0.693 Treatment Placebo         3      0
5 1        31   3.43 Week4        2   0.693 Treatment Placebo         3      1
6 2        30   3.40 Baseline     8   2.08  Baseline  Placebo        11      0
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl>
print(tail(Dados))
# A tibble: 6 × 13
  ID      age agelog sequence  time logtime active    treatment seizure visit4
  <fct> <dbl>  <dbl> <chr>    <dbl>   <dbl> <fct>     <fct>       <dbl>  <dbl>
1 58       36   3.58 Week4        2   0.693 Treatment Drug            0      1
2 59       37   3.61 Baseline     8   2.08  Baseline  Drug           12      0
3 59       37   3.61 Week1        2   0.693 Treatment Drug            1      0
4 59       37   3.61 Week2        2   0.693 Treatment Drug            4      0
5 59       37   3.61 Week3        2   0.693 Treatment Drug            3      0
6 59       37   3.61 Week4        2   0.693 Treatment Drug            2      1
# ℹ 3 more variables: baseline <dbl>, cseizure <dbl>, logcseizure <dbl>
fitglmer <- lme4::glmer(seizure ~ treatment*active + (1|ID), 
                        data=Dados, 
                        family=poisson("log"))
print(car::Anova(fitglmer))
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: seizure
                   Chisq Df Pr(>Chisq)    
treatment           0.64  1       0.42    
active           1667.49  1    < 2e-16 ***
treatment:active   18.97  1    1.3e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
out <- summary(fitglmer, corr=FALSE)
print(out)
Generalized 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 ' ' 1
print(ci <- confint(fitglmer))
Computing 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.1659
cat("RR de treatment(drug) =", exp(out$coefficients[4]), 
    "[", exp(as.numeric(ci[5,1:2])), "]")
RR de treatment(drug) = 0.7391 [ 0.6446 0.8472 ]

GEE vs. GLMM

Um exemplo de semelhança dos resultados de GEE e GLMM está em:

Conforme Shi, 2012, slide 62,

“Mixed model [GLMM] is close to GEE with exchangeable correlation: Nearly identical coefficients and slightly different SEs.”

Conforme Shi, 2012, slide 64,

“GEE (population average model):

  • On average, is there a trend in score change over time?
  • Robust: even if correlation model is wrong, SE still valid

[GLMM] Mixed effects (subject specific model):

  • Are there any differences between subjects on the trend in score change over time? (Do all subjects have the same trend over time?)
  • Advantages
    • Characterization of heterogeneity
    • Incomplete unbalanced data”

GEE x GLMM: Regressão: pHi e PaCO2

Dados <- rmcorr::bland1995
Dados$Subject <- factor(Dados$Subject)

fitgeeglm <- geepack::geeglm(PaCO2 ~ pH,
                             id=Subject,
                             family=gaussian("identity"),
                             corstr="exchangeable",
                             data=Dados)
print(anova(fitgeeglm))
Analysis of 'Wald statistic' Table
Model: gaussian, link: identity
Response: PaCO2
Terms added sequentially (first to last)

   Df   X2 P(>|Chi|)
pH  1 0.33      0.57
print(summary(fitgeeglm, corr=FALSE))

Call:
geepack::geeglm(formula = PaCO2 ~ pH, family = gaussian("identity"), 
    data = Dados, id = Subject, corstr = "exchangeable")

 Coefficients:
            Estimate Std.err Wald Pr(>|W|)
(Intercept)    9.649   8.251 1.37     0.24
pH            -0.654   1.137 0.33     0.57

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)    0.557   0.174
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha    0.321   0.158
Number of clusters:   8  Maximum cluster size: 9 
fitglmer <- lme4::glmer(PaCO2 ~ pH + (1|Subject),
                        data=Dados)
Warning in lme4::glmer(PaCO2 ~ pH + (1 | Subject), data = Dados): calling
glmer() with family=gaussian (identity link) as a shortcut to lmer() is
deprecated; please call lmer() directly
print(anova(fitglmer))
Analysis of Variance Table
   npar Sum Sq Mean Sq F value
pH    1   1.87    1.87    9.59
print(summary(fitglmer, corr=FALSE))
Linear mixed model fit by REML ['lmerMod']
Formula: PaCO2 ~ pH + (1 | Subject)
   Data: Dados

REML criterion at convergence: 80.5

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.430 -0.389  0.071  0.518  2.090 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.919    0.959   
 Residual             0.195    0.442   
Number of obs: 47, groups:  Subject, 8

Fixed effects:
            Estimate Std. Error t value
(Intercept)   17.867      4.161    4.29
pH            -1.805      0.583   -3.10
print(car::Anova(fitglmer, test.statistic="F"))
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: PaCO2
      F Df Df.res Pr(>F)   
pH 8.56  1   43.4 0.0055 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fitlmer <- lmerTest::lmer(PaCO2 ~ pH + (1|Subject),
                          data=Dados)
print(anova(fitlmer))
Type III Analysis of Variance Table with Satterthwaite's method
   Sum Sq Mean Sq NumDF DenDF F value Pr(>F)   
pH   1.87    1.87     1  43.1    9.59 0.0034 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(fitlmer, corr=FALSE))
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: PaCO2 ~ pH + (1 | Subject)
   Data: Dados

REML criterion at convergence: 80.5

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.430 -0.389  0.071  0.518  2.090 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.919    0.959   
 Residual             0.195    0.442   
Number of obs: 47, groups:  Subject, 8

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   17.867      4.161 42.857    4.29  9.9e-05 ***
pH            -1.805      0.583 43.116   -3.10   0.0034 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(car::Anova(fitlmer, test.statistic="F"))
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: PaCO2
      F Df Df.res Pr(>F)   
pH 8.56  1   43.4 0.0055 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GEE x GLMM: Regressão logística: artrite reumatóide

Dados <- readxl::read_excel("arthritis_long.xlsx")
Dados$Paciente <- factor(Dados$Paciente)
Dados$Tratamento <- factor(Dados$Tratamento)
Dados$Tratamento <- relevel(Dados$Tratamento,
                            ref="Placebo")
Dados$AlivioDor <- factor(Dados$AlivioDor)
Dados$AlivioDor <- relevel(Dados$AlivioDor,
                           ref="Baixo")
Dados$Semana <- as.numeric(Dados$Semana)

fit <- lme4::glmer(AlivioDor ~ Tratamento + (1|Paciente), 
                   data=Dados, 
                   family=binomial("logit"))
print(car::Anova(fit))
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: AlivioDor
           Chisq Df Pr(>Chisq)   
Tratamento   7.7  1     0.0055 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(fit, corr=FALSE))
Generalized 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 ' ' 1
or <- 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 ]
cat("OR de Tratamento Placebo =", 1)
OR de Tratamento Placebo = 1
Dados <- na.omit(Dados)
fitgeeglm <- try(geepack::geeglm(AlivioDor ~ Tratamento, 
                                 id=Paciente,
                                 family=binomial("logit"),
                                 corstr="exchangeable",
                                 data=Dados))
Warning in model.response(mf, "numeric"): uso de type="numeric" com uma
resposta fator será ignorado
Warning in Ops.factor(y, mu): '-' não faz sentido para fatores
Error in lm.fit(zsca, qlf(pr2), offset = soffset) : NA/NaN/Inf in 'y'

GEE x GLMM: Regressão de Poisson: epilepsia

Dados <- readxl::read_excel("epilepsy.xlsx")
Dados$ID <- factor(Dados$ID)
Dados$treatment <- factor(Dados$treatment)
Dados$treatment <- relevel(Dados$treatment, ref="Placebo")
Dados$active <- factor(Dados$active)
Dados$active <- relevel(Dados$active, ref="Baseline")

fitgeeglm <- geepack::geeglm(seizure ~ treatment*active,
                             offset=logtime,
                             id=ID,
                             family=poisson("log"),
                             corstr="exchangeable",
                             data=Dados)
# print(anova(fitgeeglm)) # 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 ' ' 1
out <- summary(fitglmer, corr=FALSE)
print(out)
Generalized 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 ' ' 1
print(ci <- confint(fitglmer))
Computing 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.166
cat("RR de treatment(drug) =", exp(out$coefficients[4]), 
    "[", exp(as.numeric(ci[5,1:2])), "]")
RR de treatment(drug) = 0.739 [ 0.645 0.847 ]

Referências

  • Baigorri, F & Calvet, X & Joseph, D (1997) Equipment review: Gastric intramucosal pH measurement. Crit Care 1: 61. https://doi.org/10.1186/cc104.
  • Bakdash, JZ & Marusich, LR (2017) Repeated measures correlation. Front. Psychol. 8:456. doi: 10.3389/fpsyg.2017.00456. https://osf.io/djphm/
  • Bland, JM & Altman, DG (1995a) Calculating correlation coefficients with repeated observations: Part 1 - Correlation within subjects. BMJ 310(6977): 446. https://doi.org/10.1136/bmj.310.6977.446
  • Bland JM & Altman DG (1995b) Calculating correlation coefficients with repeated observations: Part 2 - Correlation between subjects. BMJ 310(6980): 633. doi: 10.1136/bmj.310.6980.633. Erratum in: BMJ 1996 312(7030):572.
  • Chartier, S & Faulkner, A (2008) General linear models: An integrated approach to statistics. Tutorial in Quantitative Methods for Psychology 4(2), 65-78.
  • Clark, M (?) Mixed model with R. https://m-clark.github.io/mixed-models-with-R/.
  • Duursma, R & Powell, J (2016) Linear mixed-effects models. http://www.hiercourse.com/mixedeffects.
  • Duursma, R & Powell, J (2016) Linear mixed-effects models: Solutions. http://www.hiercourse.com/mixedeffects.
  • Faraway, JJ (2016) Extending the linear model with R: Generalized linear, mixed effects and nonparametric regression models. 2 ed. USA: CRC.
  • Francq, GF et al. (2019) Confidence, prediction, and tolerance in linear mixed models. Statistics in Medicine 38(30): 5603-22.
  • Francq, GF et al. (2020) Erratum: Confidence, prediction, and tolerance in linear mixed models. Statistics in Medicine 39(14): 2015-2015.
  • Hamaker, EL (2013) Why researchers should think “within-person”: A paradigmatic rationale. In Mehl MR &, Conner TS (Eds.) (2013) Handbook of research methods for studying daily life. NY: Guilford. p. 43-61.
  • Hothorn, T & Everitt, BS (2014) A handbook of statistical analyses using R. 3rd ed. Chapman & Hall/CRC, 2014. https://rdrr.io/cran/HSAUR3/. R package HSAUR3.
  • Judkins DR & Porter KE (2016) Robustness of ordinary least squares in randomized clinical trials. Stat Med. 35(11):1763-73. doi: 10.1002/sim.6839.
  • Kirkwood, BR & Sterne, JAC (2003) Essential medical statistics. 2nd ed. USA: Blackwell. Chapter 13: Analysis of clustered data.
  • Konrad, M (2021) Clustered standard errors with R. https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/.
  • Laird, NM and Ware, JH (1982) Random-effects models for longitudinal data. Biometrics 38(4): 963-74.
  • Magnusson, K (2015) Using R and lme/lmer to fit different two- and three-level longitudinal models. https://rpsychologist.com/r-guide-longitudinal-lme-lmer.
  • Mixed effects logistic regression: R data analysis examples. https://stats.oarc.ucla.edu/r/dae/mixed-effects-logistic-regression/.
  • Norusis, MJ (2012) IBM SPSS statistics 19 advanced statistical procedures companion. NJ: Prentice-Hall.
  • Oikawa, KF (2019) Análise estatística de dados longitudinais e hierárquicos em psicologia: uma análise comparativa entre GEE e GLMM. Tese de Doutorado, Psicologia Experimental, Instituto de Psicologia, Universidade de São Paulo, São Paulo. doi:10.11606/T.47.2019.tde-18112019-183159. Recuperado em 2022-05-08, de www.teses.usp.br.
  • Oliveira, AA (Coord.) (?) Modelos generalizados mistos (GLMM). Laboratório de Ecologia de Florestas Tropicais do Instituto de Biocências da USP. http://labtrop.ib.usp.br/doku.php?id=cursos:planeco:roteiro:12-glmm.
  • Pekár, S & Brabec, M (2017) Generalized estimating equations: A pragmatic and flexible approach to the marginal GLM modelling of correlated data in the behavioural sciences. Ethology 124:86–93. DOI: 10.1111/eth.12713.
  • Sainani, K. (2010) The importance of accounting for correlated observations. Physical Medicine and Rehabilitation 2(9):858-61. doi: 10.1016/j.pmrj.2010.07.482.
  • Shi, L (2012) Application of GEE and Mixed Effects model in longitudinal data analysis. Umass Boston: Dana-Farber/Harvard Cancer Center.
  • 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.