Modelos Lineares Generalizados

Desenho Experimental e Análise Avançada de Dados Ecológicos

João O. Santos

joao.santos@ispa.pt

ISPA

A Mesma Pergunta

Como muda Y com X?

O Que Muda?

a resposta já não é bem representada por um modelo Gaussiano

Três Peças

preditor linear · família · ligação

O Preditor Linear

\[ \eta=X\beta \]

Interceptos, declives, factores e interacções continuam cá.

A Família

Que valores pode Y tomar — e como varia a sua incerteza?

A Ligação

Em que escala o preditor linear se liga à média da resposta?

Fixação de Larvas

Exemplo simulado

n.º de larvas fixadas / n.º de larvas expostas

Uma linha = uma placa.

Há cinco placas por substrato e 50 larvas expostas por placa. Pedir que identifiquem uma placa: a combinação substrato + número distingue-a das outras. A simulação supõe placas distintas, sem um tanque partilhado. Num estudo real, que placas partilham tanque e como varia a fixação entre placas? A família binomial continua a exigir que avaliemos a independência das larvas.

Conhecer as Placas

ds <- read.csv("../../data/larval_settlement_demo.csv")
ds$proportion <- ds$settled / ds$total
head(ds)

  substrate plate total settled proportion
1         A     1    50      18       0.36
2         B     1    50      25       0.50
3         C     1    50      23       0.46
4         A     2    50      17       0.34
5         B     2    50      25       0.50
6         C     2    50      24       0.48

CSV simulado · por placa: settled = número fixado; total = número exposto.

O Denominador

7/10 e 70/100 não contêm a mesma informação

Ver as Réplicas

Os Factores Continuam Cá

options(contrasts = c("contr.sum", "contr.poly"))

Soma-zero no preditor linear, não na escala da resposta.

Definir antes dos modelos: contr.sum para factores não ordenados e contr.poly para ordenados, nesta sessão. No binomial seguinte, os coeficientes de substrato são desvios na escala logit, não diferenças de probabilidades. No Poisson com dois habitats, o coeficiente é metade da diferença de log-médias (log-taxas com o offset): exponenciá-lo não dá a razão completa entre habitats. Usar emmeans para comparações na direcção pretendida. Contrastes e GLMs.

Binomial

m_prop <- glm(proportion ~ substrate, weights = total,
              data = ds, family = binomial())

Logit

\[ \operatorname{logit}(p)= \log\left(\frac{p}{1-p}\right)=X\beta \]

Voltar à Biologia

settlement_means <- emmeans(m_prop, ~ substrate)
summary(settlement_means, type = "response")

 substrate  prob     SE  df asymp.LCL asymp.UCL
 A         0.364 0.0304 Inf     0.307     0.425
 B         0.620 0.0307 Inf     0.558     0.678
 C         0.520 0.0316 Inf     0.458     0.581

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

summary(pairs(settlement_means), type = "response", infer = TRUE)

 contrast odds.ratio     SE  df asymp.LCL asymp.UCL null z.ratio p.value
 A / B         0.351 0.0649 Inf     0.227     0.541    1  -5.660 <0.0001
 A / C         0.528 0.0964 Inf     0.344     0.810    1  -3.496  0.0014
 B / C         1.506 0.2740 Inf     0.984     2.305    1   2.254  0.0625

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
Intervals are back-transformed from the log odds ratio scale 
P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale 

probabilidades e IC 95% · contrastes: razões de odds

weights = total conserva o número de larvas por placa; não ignoramos o denominador.

Testes no log odds, com ajuste de Tukey para as 3 comparações.

Mesma Arquitectura, Outra Família

Exemplo simulado · transectos

Quantos organismos vimos?

Poisson

ds <- read.csv("../../data/transect_counts_demo.csv")
m_count <- glm(n ~ habitat, data = ds, family = poisson())

Mas o Esforço Varia

20 organismos em 10 min ≠ 20 em 30 min

Offset

ds <- ds |> mutate(log_effort = log(effort))
m_rate <- glm(n ~ habitat + offset(log_effort),
              data = ds, family = poisson())

effort tem de ser positivo; log_effort é o seu logaritmo. O offset fixa o seu coeficiente em 1: à mesma taxa, triplicar o tempo triplica a contagem esperada. Não estimamos um declive para o esforço.

Neste cenário, effort é a duração em minutos, sempre positiva. A resposta continua a ser a contagem n; o offset fixa em 1 o coeficiente de log(effort), assumindo contagem esperada proporcional ao esforço. Para a mesma taxa, triplicar o tempo triplica a contagem esperada. Distinguir esta suposição de estimar livremente um efeito da duração.

Escala Logarítmica

\[ \log E(Y|X)=X\beta \]

Exponenciar uma diferença de log-médias produz uma razão de médias.

Família Também É Variância

A escolha não serve apenas para impedir previsões impossíveis.

Extensão: Sobredispersão

mais variação do que o Poisson espera

Perguntem porquê antes de mudar de família.

Pode Ser

heterogeneidade · agrupamento · zeros · média mal especificada

GLM ≠ Modelo Misto

Mudar a família não resolve dependência entre linhas.

Demonstrar a Comparação Global

Voltamos aos dados de fixação de larvas. Demonstração: escrever a probabilidade comum à parte mostra a hipótese que a linha substrate testa.

ds <- read.csv("../../data/larval_settlement_demo.csv")
ds$proportion <- ds$settled / ds$total
m0 <- glm(proportion ~ 1, weights = total, data = ds, family = binomial())
anova(m0, m_prop, test = "Chisq")

Analysis of Deviance Table

Model 1: proportion ~ 1
Model 2: proportion ~ substrate
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        14     48.373                          
2        12     14.705  2   33.667 4.889e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Comparação: probabilidade comum vs. uma por substrato; dois parâmetros adicionais
• Teste: razão de verosimilhanças (deviâncias), não o F Gaussiano
• No script: Anova() no modelo ajustado; modelo nulo só na demonstração

A interpretação como modelos compacto/aumentado mantém-se, mas nem todas as tabelas usam a mesma estatística. Em GLMs, anova(…, test = “Chisq”) pede um teste de razão de verosimilhanças; summary() mostra testes de Wald dos coeficientes. Não prometer que os respectivos valores-p sejam idênticos. A comparação explícita serve para explicar o teste global. O script final pede o mesmo teste de razão de verosimilhanças com Anova(…, test.statistic = “LR”), sem ajustar separadamente o modelo nulo.

Uma Tabela de Frequências

Os mesmos resultados, agora como substrato × fixada/não fixada.

ds$failed <- ds$total - ds$settled
frequencies <- xtabs(cbind(settled, failed) ~ substrate, data = ds)
frequencies

         
substrate settled failed
        A      91    159
        B     155     95
        C     130    120

xtabs() soma as contagens por substrato, sem criar uma linha por larva. Agregar não cria independência: precisamos de a justificar para larvas e placas.

Qui-Quadrado de Pearson

chisq.test(frequencies)

    Pearson's Chi-squared test

data:  frequencies
X-squared = 33.291, df = 2, p-value = 5.902e-08

Mesma hipótese: a probabilidade de fixação é igual nos três substratos. Pearson compara frequências observadas e esperadas sob essa hipótese.

A Mesma Hipótese no Binomial

Ainda na demonstração: dois testes da mesma hipótese, não dois passos a repetir na análise.

anova(m0, m_prop, test = "Chisq")

Analysis of Deviance Table

Model 1: proportion ~ 1
Model 2: proportion ~ substrate
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        14     48.373                          
2        12     14.705  2   33.667 4.889e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pearson \(X^2\) ≠ razão de verosimilhanças \(G^2\). Testam a mesma hipótese e têm a mesma referência \(\chi^2_2\) assintótica; as estatísticas e os valores-p não são exactamente iguais em amostras finitas.

Aqui, o GLM de substrato é saturado para as três proporções de grupo. Agregar por substrato conserva os totais e a comparação de verosimilhanças entre estes dois modelos, sob o mesmo modelo binomial independente. A equivalência é de hipótese e, assintoticamente, de inferência; não é uma identidade numérica entre Pearson e a diferença de deviâncias. Não confundir estes testes aproximados com um teste exacto de Fisher ou binomial exacto. R documenta estas estatísticas em help(“chisq.test”) e help(“anova.glm”).

Script Binomial

Descarreguem o CSV simulado e adaptem o caminho em read.csv() à vossa pasta de trabalho. Um modelo: Anova() dá o teste global; model_parameters() descreve os coeficientes e emmeans() dá estimativas na escala biológica.

library(car)
library(dplyr)
library(emmeans)
library(ggplot2)
library(parameters)

options(contrasts = c("contr.sum", "contr.poly"))

ds <- read.csv("../../data/larval_settlement_demo.csv") |>
  mutate(proportion = settled / total)
m <- glm(proportion ~ substrate, weights = total, data = ds, family = binomial())
settlement_means <- emmeans(m, ~ substrate)
g <- ggplot(ds, aes(substrate, proportion)) + geom_point()

# Print results and plot
model_parameters(m) |> print()
Anova(m, test.statistic = "LR") |> print()
summary(settlement_means, type = "response") |> print()
pairs(settlement_means) |> summary(type = "response", infer = TRUE) |> print()
print(g)

Resultados Binomiais

Parameter     | Log-Odds |   SE |         95% CI |     z |      p
-----------------------------------------------------------------
(Intercept)   | 3.85e-03 | 0.07 | [-0.14,  0.15] |  0.05 | 0.959 
substrate [1] |    -0.56 | 0.11 | [-0.77, -0.35] | -5.28 | < .001
substrate [2] |     0.49 | 0.11 | [ 0.28,  0.70] |  4.58 | < .001

Analysis of Deviance Table (Type II tests)

Response: proportion
          LR Chisq Df Pr(>Chisq)    
substrate   33.667  2  4.889e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 substrate  prob     SE  df asymp.LCL asymp.UCL
 A         0.364 0.0304 Inf     0.307     0.425
 B         0.620 0.0307 Inf     0.558     0.678
 C         0.520 0.0316 Inf     0.458     0.581

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

 contrast odds.ratio     SE  df asymp.LCL asymp.UCL null z.ratio p.value
 A / B         0.351 0.0649 Inf     0.227     0.541    1  -5.660 <0.0001
 A / C         0.528 0.0964 Inf     0.344     0.810    1  -3.496  0.0014
 B / C         1.506 0.2740 Inf     0.984     2.305    1   2.254  0.0625

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
Intervals are back-transformed from the log odds ratio scale 
P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale 

Ler o Teste Global

Anova(m_prop, test.statistic = "LR")

Analysis of Deviance Table (Type II tests)

Response: proportion
          LR Chisq Df Pr(>Chisq)    
substrate   33.667  2  4.889e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

substrate: as três probabilidades de fixação são iguais? LR pede o teste de razão de verosimilhanças: é o da demonstração anterior. Não é necessário voltar a ajustar um modelo de probabilidade comum.

Ler os Coeficientes

model_parameters(m_prop)

Parameter     | Log-Odds |   SE |         95% CI |     z |      p
-----------------------------------------------------------------
(Intercept)   | 3.85e-03 | 0.07 | [-0.14,  0.15] |  0.05 | 0.959 
substrate [1] |    -0.56 | 0.11 | [-0.77, -0.35] | -5.28 | < .001
substrate [2] |     0.49 | 0.11 | [ 0.28,  0.70] |  4.58 | < .001

Coeficientes na escala logit, não probabilidades nem diferenças de probabilidade. Os testes z dizem respeito a esses coeficientes.

Ler na Escala Biológica

settlement_means <- emmeans(m_prop, ~ substrate)
summary(settlement_means, type = "response")

 substrate  prob     SE  df asymp.LCL asymp.UCL
 A         0.364 0.0304 Inf     0.307     0.425
 B         0.620 0.0307 Inf     0.558     0.678
 C         0.520 0.0316 Inf     0.458     0.581

Confidence level used: 0.95 
Intervals are back-transformed from the logit scale 

summary(pairs(settlement_means), type = "response", infer = TRUE)

 contrast odds.ratio     SE  df asymp.LCL asymp.UCL null z.ratio p.value
 A / B         0.351 0.0649 Inf     0.227     0.541    1  -5.660 <0.0001
 A / C         0.528 0.0964 Inf     0.344     0.810    1  -3.496  0.0014
 B / C         1.506 0.2740 Inf     0.984     2.305    1   2.254  0.0625

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
Intervals are back-transformed from the log odds ratio scale 
P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale 

• Primeiro: probabilidades previstas por substrato
• Depois: razões de odds, não diferenças de probabilidade
• Tukey: ajuste para as três comparações entre substratos

Relatar a Fixação

A probabilidade estimada de fixação foi .364 em A, 95% CI [.307, .425], .620 em B, 95% CI [.558, .678], e .520 em C, 95% CI [.458, .581]. As odds em A foram 0.351 vezes as de B, 95% CI [0.227, 0.541], com ajuste de Tukey. O teste global de razão de verosimilhanças comparou uma probabilidade comum com três probabilidades: \(\chi^2_2\) = 33.67, P < .001. Os testes z do resumo são de Wald, não este teste global. O exemplo é simulado; a inferência depende da independência assumida entre larvas e placas.

Script de Taxas

Descarreguem o CSV simulado e adaptem o caminho em read.csv() à vossa pasta de trabalho.

library(dplyr)
library(emmeans)
library(ggplot2)
library(parameters)

options(contrasts = c("contr.sum", "contr.poly"))

ds <- read.csv("../../data/transect_counts_demo.csv") |>
  mutate(log_effort = log(effort))
m <- glm(n ~ habitat + offset(log_effort), data = ds, family = poisson())
rate_means <- emmeans(m, ~ habitat, offset = 0)
g <- ggplot(ds, aes(habitat, n / effort)) + geom_point()

# Print results and plot
model_parameters(m) |> print()
summary(rate_means, type = "response") |> print()
pairs(rate_means) |> summary(type = "response", infer = TRUE) |> print()
print(g)

Resultados das Taxas

Parameter   | Log-Mean |   SE |         95% CI |      z |      p
----------------------------------------------------------------
(Intercept) |    -1.21 | 0.07 | [-1.36, -1.07] | -16.35 | < .001
habitat [1] |    -0.34 | 0.07 | [-0.48, -0.19] |  -4.54 | < .001

 habitat  rate     SE  df asymp.LCL asymp.UCL
 areia   0.212 0.0254 Inf     0.168     0.268
 rochoso 0.416 0.0366 Inf     0.350     0.495

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

 contrast        ratio     SE  df asymp.LCL asymp.UCL null z.ratio p.value
 areia / rochoso  0.51 0.0757 Inf     0.381     0.682    1  -4.539 <0.0001

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 
Tests are performed on the log scale 

offset = 0 em emmeans() pede previsões para uma unidade de esforço.

Ler Contagens Como Taxas

model_parameters(m_rate)

Parameter   | Log-Mean |   SE |         95% CI |      z |      p
----------------------------------------------------------------
(Intercept) |    -1.21 | 0.07 | [-1.36, -1.07] | -16.35 | < .001
habitat [1] |    -0.34 | 0.07 | [-0.48, -0.19] |  -4.54 | < .001

rate_means <- emmeans(m_rate, ~ habitat, offset = 0)
summary(rate_means, type = "response")

 habitat  rate     SE  df asymp.LCL asymp.UCL
 areia   0.212 0.0254 Inf     0.168     0.268
 rochoso 0.416 0.0366 Inf     0.350     0.495

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

summary(pairs(rate_means), type = "response", infer = TRUE)

 contrast        ratio     SE  df asymp.LCL asymp.UCL null z.ratio p.value
 areia / rochoso  0.51 0.0757 Inf     0.381     0.682    1  -4.539 <0.0001

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 
Tests are performed on the log scale 

A resposta ajustada continua a ser a contagem. As previsões para um minuto são taxas em organismos/minuto; o contraste é uma razão de taxas.

Relatar as Taxas

Neste exemplo simulado, as taxas estimadas foram 0.212 organismos/minuto em areia, 95% CI [0.168, 0.268], e 0.416 em habitat rochoso, 95% CI [0.350, 0.495]. A razão areia/rochoso foi 0.510, 95% CI [0.381, 0.682], teste de Wald z = -4.54, P < .001. O GLM de Poisson incluiu o logaritmo da duração como offset.

Antes de generalizar: verificar sobredispersão, agrupamento e a suposição de contagem proporcional ao esforço. O offset não confirma essas condições.

Artigo da Aula

O’Brien et al. (2025)

Light and Dark Biofilm Adaptation Impacts Larval Settlement in Diverse Coral Species

Escolha do Modelo

20 minutos · grupos de 3 ou 5

Procurem Statistical analysis of settlement assays.

1. Qual é a resposta e qual é o denominador?
2. Porque faz sentido uma família binomial?
3. Que observações partilham o mesmo tanque?
4. Porque é que um GLMM pode ser necessário mesmo depois de escolher a família certa?

Agora: Tutorial

20 minutos · grupos de 3 ou 5

GLMs no R

Façam “emmeans em Logística”.

Separem probabilidades previstas de odds ratios.

Se Sobrar Tempo

GLMs no R

Façam “Manatins”.

Separem associação, teste do efeito, dispersão e limites causais.

Decidir a Unidade

10 minutos · grupos de 3 ou 5

Para cada caso, escrevam o que representa uma linha e a informação que não pode ser omitida:

• larvas fixadas em cada placa;
• peixes contados em transectos com durações diferentes.

Pergunta Extra

Mesmo depois de lidar com sobredispersão, que problema temporal pode continuar?

Se Precisarem de Teoria

Modelos Lineares Generalizados

A Seguir

Finalmente: o que fazemos quando várias linhas pertencem à mesma unidade?