Objetivos

Esta apresentação tem como propósito:

1. Mostrar a evolução metodológica Regressão clássica → GLM → GAM → GAMLSS**, destacando motivações e limitações.
2. Formalizar o GAMLSS: distribuição, parâmetros, ligações e preditores aditivos.
3. Descrever o princípio de estimação por verossimilhança penalizada.
4. Apresentar seleção de modelo em GAMLSS.
5. Discutir diagnóstico: resíduos quantílicos randomizados e worm plot.
6. Apresentar uma breve aplicação em R.

Linha histórica: limitações e avanços

• Regressão clássica (Gauss/Legendre): erro Normal, homoscedasticidade, média linear.
• GLM (Nelder & Wedderburn, 1972): respostas não Normais via família exponencial e função de ligação (Nelder & Wedderburn, 1972).
• GAM (Hastie & Tibshirani, 1986): relações não lineares na média via preditores aditivos suavizados (Hastie & Tibshirani, 1986).
• GAMLSS (Rigby & Stasinopoulos, 2005): modelagem distributiva: parâmetros de localização, escala e forma podem variar com covariáveis (Rigby & Stasinopoulos, 2005).

Regressão clássica

Modelo linear clássico:

\[ Y = X\beta + \varepsilon,\qquad \varepsilon \sim \mathcal{N}(0,\sigma^2 I). \]

Pressupostos centrais:

• Linearidade: \(E(Y\mid X)=X\beta\).
• Homoscedasticidade: \(\mathrm{Var}(Y\mid X)=\sigma^2\), constante.
• Normalidade e independência dos erros.

Em aplicações, essas hipóteses podem falhar: respostas positivas e assimétricas, contagens, proporções, e variância dependente do nível médio.

Motivação para GLM

Os GLMs surgem para lidar com:

1. Distribuição da resposta: \(Y\) pode não ser adequadamente Normal.
2. Variância ligada à média: em várias famílias, a variância condicional cresce com a média.

A meta é manter interpretabilidade com flexibilidade para a distribuição da resposta.

GLM: Modelos Lineares Generalizados

Assume-se \(Y\) na família exponencial:

\[ f(y\mid\theta,\phi)=\exp\left(\frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right). \]

Ligação para a média:

\[ g(\mu_i)=\eta_i=x_i^\top\beta,\qquad \mu_i=E(Y_i\mid X). \]

Variância:

\[ \mathrm{Var}(Y_i\mid X)=\phi\,V(\mu_i). \]

Interpretação: o GLM flexibiliza a distribuição de \(Y\), mas impõe que \(g(\mu)\) seja linear em covariáveis.

Dados de exemplo: rent

Usaremos a base de dados disponível em R: rent (aluguel em Munique).

O diagrama de dispersão sugere uma associação positiva entre a área do imóvel e o valor do aluguel. Também é possível observar que a variabilidade do aluguel aumenta conforme a área cresce, indicando heteroscedasticidade.

Regressão clássica: diagnóstico

• Gráfico de resíduos vs. ajustados: padrão claro de aumento da dispersão dos resíduos conforme os valores ajustados crescem - indício de heteroscedasticidade;

• QQ-plot normal: os pontos se afastam da reta principalmente nas extremidades, indicando que os resíduos não seguem uma distribuição nornam e presença de assimetria e caudas mais pesadas.

Modelo Linear Generalizado (utilizando Gamma)

Para \(Y>0\), um GLM Gamma com função de ligação logarítmica é comum:

\[ Y_i \sim \mathrm{Gamma}(\mu_i,\phi),\qquad \log(\mu_i)=x_i^\top\beta. \]

• gráfico deviance vs. ajustados: ainda é perceptível a mudança na dispersão dos resíduos ao longo dos valores ajustados, a variabilidade parece ser maior em dentre ajustes mais altos. Isso indica que pode haver heteroscedasticidade que não está sendo completamente capturada por um único parâmetro global de dispersão.

• QQ-plot: apresenta desvios nas caudas, indicando que o modelo ainda não descreve perfeitamente o comportamento dos dados.

• Conclusão: o GLM (Gamma) foi um avanço importante, já que lida bem com positividade e assimetria, mas ainda é rígiso por exigie linearidade da função de ligação, e tratar a dispersão de forma global.

Motivação para GAM

O GAM é motivado quando:

• Mesmo com distribuição e função de ligação adequados, o GLM ainda impõe uma relação linear entre o preditor linear e as covariáveis, o que pode ser rígido;

• A ideia do GAM é manter distribuição da variável resposta e função de ligação, mas substituir o preditor linear por funções suaves estimadas com penalização, permitindo captar não linearidades.

GAM: Modelos Aditivos Generalizados

No GAM:

\[ g(\mu_i)=\eta_i=\beta_0+\sum_{j=1}^p f_j(x_{ij}), \]

onde \(f_j\) são funções suaves.

• continua assumindo uma distribuição para a distribuição condicional da variável resposta;

• continua usando função de ligação para relacionar a média ao preditor;

• diferença: ao invés do preditor ser linear, assume preditor aditivo;

• as funções \(f_j\) são representadas por bases spline e estimadas com penalização, permitindo capturar não linearidade;

• a dispersão permanece global.

Exemplo em R: GAM (Gamma)

• Efeito suave: o eixo vertical representa a contribuição do termo suave para o preditor \(g(\mu)\), e não o aluguel na escala original. Aqui, a curva estimada é crescente e praticamente linear ao longo do intervalo, com pequenas curvaturas. Isso sugere que o GAM não identificou uma não linearidade forte;

• Deviance vs ajustados:, embora os resíduos estejam em torno de zero, a variabilidade dos resíduos ainda não é uniforme ao longo dos valores ajustados. Isso indica inadequação na forma como a dispersão é tratada, já que o GAM mantém a dispersão como um componente global.

Limitações do GAM e motivação para GAMLSS

O GAM resolve uma limitação importante do GLM, visto que permite que a relação entre covariáveis e a média condicional seja não linear, utilizando funções suaves. Porém, algumas limitações permanecem, que são a motivação para o GAMLSS:

1. A variância pode não ser constante e nem explicada apenas como uma função fixa da média. Ou seja, para determinados valores das covariáveis, a resposta é tem mais variabilidade do que para outros. No GAM, a dispersão é tratada de forma global;

2. A forma da distribuição pode mudar: a distribuição pode ficar mais assimétrica em ou caudas mais pesadas em certas faixas de \(X\). Os GAMs não modelam parâmetros de forma, já que eles atuam principalmente sobre a média;

3. Há cenários em que a família usada no GLM/GAM é limitada. Exemplos: inflação de zeros, truncamento, caudas muito pesadas ou distribuições que não pertencem à família exponencial.

GAMLSS: definição e ideia central

Assume-se:

\[ Y_i \mid x_i \sim D(\boldsymbol{\theta}_i),\qquad \boldsymbol{\theta}_i=(\mu_i,\sigma_i,\nu_i,\tau_i), \]

onde \(D\) não precisa pertencer à família exponencial.

• \(\mu\): parâmetro de localização,

• \(\sigma\): escala (dispersão),

• \(\nu\): assimetria,

• \(\tau\): curtose/caudas,

dependendo da parametrização de \(D\).

Para cada parâmetro \(\theta_k\), define-se uma função de ligação \(g_k\) e um preditor aditivo:

\[ g_k(\theta_{ki})=\eta_{ki}=\beta_{0k}+\sum_{j=1}^{p_k} f_{jk}(x_{ij}), \qquad k\in\{\mu,\sigma,\nu,\tau\}. \]

O GAMLSS é permite que diferentes aspectos da distribuição condicional variem com covariáveis (Rigby & Stasinopoulos, 2005; Stasinopoulos & Rigby, 2007).

GAMLSS como modelo semiparamétrico

O GAMLSS é descrito como semiparamétrico porque:

• paramétrico: impõe uma forma probabilística \(D\) para \(Y\mid X\);
• não paramétrico: permite que componentes dos preditores sejam funções suaves \(f(\cdot)\), estimadas via bases spline e penalização.

Isso flexibiliza ao mesmo tempo que controla a complexidade (via suavização e seleção de modelo).

Componentes do modelo: \(\mathcal{M}=\{D,G,T,\mathcal{L}\}\)

Organização útil:

\[ \mathcal{M}=\{D, G, T, \mathcal{L}\}, \]

• \(D\): distribuição de \(Y\mid X\).
• \(G\): ligações \(g_\mu,g_\sigma,g_\nu,g_\tau\).
• \(T\): termos nos preditores (covariáveis, interações, funções suaves) para \(\mu,\sigma,\nu,\tau\).
• \(\mathcal{L}\): parâmetros de suavização (intensidade da penalização) em cada preditor.

Ajustar um GAMLSS exige especificar (ou selecionar) esses quatro componentes.

O componente \(D\): escolha da distribuição da resposta

A seleção de \(D\) deve considerar:

• coerência com o suporte (exemplo: \(Y>0\), \(Y\in(0,1)\), contagens, zeros);
• considerações de assimetria e curtose;
• comparação por GAIC/AIC/BIC entre as famílias adequadas;
• validação por resíduos quantílicos e worm plot.

O componente \(G\): funções de ligação

Em GAMLSS, cada parâmetro \(\theta_k \in \{\mu,\sigma,\nu,\tau\}\) é modelado via \[ g_k(\theta_{k,i}) = \eta_{k,i}. \]

As funções de ligação \(g_k\) garantem:

1. domínio válido do parâmetro (exemplo: \(\sigma>0\Rightarrow \log(\sigma)\));
2. estabilidade numérica e interpretabilidade.

O componente \(T\): quais parâmetros variam com covariáveis

No GAMLSS, é necessário escolher um preditor para cada parâmetro \(\theta_k\in\{\mu,\sigma,\nu,\tau\}\). Assim, precisamos responder a pergunta “quais parâmetros dependem de covariáveis e quais permanecem constantes?”.

O procedimento recomendado é começar por \(\mu\), em seguida \(\sigma\) quando a variabilidade muda conforme os valores de \(X\). Por final \(\nu\) e \(\tau\), já que incluir covariáveis em parâmetros de forma aumenta muito a complexidade do de modelo e pode dificultar interpretação e convergência.

\[ \mu \;\rightarrow\; \sigma \;\rightarrow\; \nu \;\rightarrow\; \tau. \]

Assim, o componente \(T\) deve ser escolhido de forma parcimoniosa, utilizando os critérios GAIC, AIC, BIC, e diagnóstico.

Estimação: verossimilhança penalizada

Quando os preditores incluem termos suaves (splines), o ajuste do GAMLSS é formulado como a maximização de uma log-verossimilhança penalizada, que equilibra qualidade de ajuste e complexidade das funções estimadas:

\[ \ell_p(\boldsymbol{\beta}) = \ell(\boldsymbol{\beta}) -\frac{1}{2}\sum_{k\in\{\mu,\sigma,\nu,\tau\}} \sum_{j\in\mathcal{S}_k} \lambda_{kj}\, \boldsymbol{\beta}_{kj}^{\top}\mathbf{P}_{kj}\boldsymbol{\beta}_{kj}. \]

• \(\ell(\boldsymbol{\beta})=\sum_{i=1}^{n}\log f(y_i\mid \boldsymbol{\theta}_i)\) é a log-verossimilhança, com
  \(\boldsymbol{\theta}_i=(\mu_i,\sigma_i,\nu_i,\tau_i)\) dependente de covariáveis via as ligações \(g_k(\theta_{k,i})=\eta_{k,i}\).
• \(\mathcal{S}_k\) são os termos suaves presentes no preditor do parâmetro \(k\).
• \(\mathbf{P}_{kj}\) é a matriz de penalização associada ao suavizador: ela quantifica a rugosdade/complexidade do efeito suave.
• \(\lambda_{kj}\ge 0\) é o parâmetro de suavização, que controla o equilíbrio de ajuste/suavidade:
  \(\lambda\) grande \(\Rightarrow\) funções mais suaves;
  \(\lambda\) pequeno \(\Rightarrow\) maior flexibilidade/ondulações.

Estimação iterativa

O GAMLSS possui preditores distintos para cada parâmetro: \[ g_\mu(\mu_i)=\eta_{\mu,i},\quad g_\sigma(\sigma_i)=\eta_{\sigma,i},\quad g_\nu(\nu_i)=\eta_{\nu,i},\quad g_\tau(\tau_i)=\eta_{\tau,i}. \] Assim, a maximização da log-versossimilhança é feita por ciclos de atualização:

O backfitting é um esquema iterativo de ajuste por blocos, onde se atualiza um conjunto de parâmetos mantendo os demais fixos. O ciclo é mostrado abaixo:

1. Fixar os preditores atuais de \((\sigma,\nu,\tau)\) e atualizar o preditor de \(\mu\);
2. Fixar \((\mu,\nu,\tau)\) e atualizar \(\sigma\);
3. Fixar \((\mu,\sigma,\tau)\) e atualizar \(\nu\);
4. Fixar \((\mu,\sigma,\nu)\) e atualizar \(\tau\);
5. Repetir o ciclo até convergência.

A convergência ocorre quando as mudanças em \(\ell_p\) e/ou nos preditores tornam-se muito pequenas, indicando estabilização.

(Rigby & Stasinopoulos, 2005; Stasinopoulos & Rigby, 2007).

Critérios de informação: GAIC, AIC e BIC-like

Na seleção de modelos em GAMLSS, são utilizados critérios de informação que equilibram qualidade de ajuste e complexidade. Um cirtério bastante utilizado é o GAIC (Generalized Akaike Information Criterion):

\[ GAIC(k)= -2\,\ell(\hat{\boldsymbol{\theta}}) + k\cdot df, \]

em que:

• \(\ell(\hat{\boldsymbol{\theta}})\) é a log-verossimilhança maximizada do modelo;

• \(df\) representa os graus de liberdade do modelo, incluido parâmetros e a complexidade dos termos suaves;

• \(k\) controla a força da penalização por complexidade.

• \(k=2 \Rightarrow GAIC(2)\) coincide com o AIC (penalização mais branda, favorecendo modelos mais flexíveis).

• \(k=\log(n) \Rightarrow\) critério BIC-like (penalização mais forte, favorecendo modelos mais parcimoniosos).

Diagnóstico: resíduos quantílicos randomizados (RQR)

• No GAMLSS, resíduos clássicos podem não ser comparáveis para diferentes observações, já que a variância e a forma da distribuição podem variar conforme as covariáveis.

• Para adotar uma padronização, usamos o resíduo quantílico randomizado (RQR), que transforma cada observação para a escala Normal padrão.

Para cada \(i\), calcula-se o valor sob o modelo ajustado:

1. \[ u_i = F_Y\!\left(y_i \mid \hat{\boldsymbol{\theta}}_i\right), \] onde \(F_Y(\cdot\mid \hat{\boldsymbol{\theta}}_i)\) é a acumulada da distribuição assumida para \(Y_i\mid X_i\).

2. \[ r_i = \Phi^{-1}(u_i). \]

Diagnóstico: worm plot

• O worm plot é uma extensão do QQ-plot aplicada aos resíduos quantílicos randomizados. Nele, para cada faixa de uma covariável, é apresentado o desvio entre os quantis observados e os quantis teóricos. Assim, em vez de observar apenas um QQ-plot global, o worm plot permite investigar em quais regiões o modelo é inadequado.

• Muito útil para o GAMLSS, já que \(\mu(x)\), \(\sigma(x)\), \(\nu(x)\), \(\tau(x)\) podem variar com \(x\). Ou seja, é possível que o ajuste seja bom em uma parte do domínio e ruim em outra.

• Ajuste adequado: os pontos distribuídos aleatoriamente em torno de zero, dentro das bandas de referência.

(Buuren & Fredriks, 2001)

Aplicação em R: GAMLSS

(1) Seleção inicial de \(D\):

• O fitDist() realiza uma avaliação inicial de algumas distribuições e seleciona a família que minimiza o critério de informação definido por \(k\). Aqui, como \(k=2\), o critério é o AIC.

• A família selecionada foi BCCG (Box–Cox–Cole–Green), muito utilizada em cenários com resposta positiva e assimetria O parâmetro \(\nu\) permite capturar a forma da distribuição.

• O Global Deviance corresponde a \(-2\ell(\hat{\theta})\), o AIC e o BIC avaliam ajuste e complexidade (aqui \(df=3\), pois são estimados \(\mu,\sigma,\nu\)).

## 
## Family:  c("BCCG", "Box-Cox-Cole-Green") 
## Fitting method: "nlminb" 
## 
## Call:  gamlssML(formula = y, family = DIST[i]) 
## 
## Mu Coefficients:
## [1]  749.3
## Sigma Coefficients:
## [1]  -0.752
## Nu Coefficients:
## [1]  0.2531
## 
##  Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   1966 
## Global Deviance:     28607.7 
##             AIC:     28613.7 
##             SBC:     28630.4

(2) Modelo inicial (BCCG): \(\mu(Fl)\), e \(\sigma\), \(\nu\) constante

Aqui, ajustamos um GAMLSS usando a distribuição BCCG para o aluguel.
A ideia é começar com um modelo simples. Modelamos \(\mu\) como uma função suave de \(Fl\) (área), usando dispersão e assimetria constantes.

As linhas de Global Deviance mostram o processo iterativo de ajuste.
Como o valor vai estabilizando há indicação de que o algoritmo convergiu.

No próximo passo, vamos permitir que \(\sigma\) dependa de \(Fl\) e ver se o critério melhora.

## GAMLSS-RS iteration 1: Global Deviance = 28138.67 
## GAMLSS-RS iteration 2: Global Deviance = 28080.4 
## GAMLSS-RS iteration 3: Global Deviance = 28078.86 
## GAMLSS-RS iteration 4: Global Deviance = 28078.71 
## GAMLSS-RS iteration 5: Global Deviance = 28078.69 
## GAMLSS-RS iteration 6: Global Deviance = 28078.69 
## GAMLSS-RS iteration 7: Global Deviance = 28078.69

## [1] 28088.53

(3) Avaliando \(\sigma(Fl)\)

Agora vamos manter a mesma forma para a média \(\mu(Fl)\), mas permitir que a dispersão (escala) também dependa de \(Fl\). Aqui, \(\sigma\) passa a ser uma função suave da área. Isso é útil quando os dados mostram que a variabilidade aumenta ou diminui conforme \(Fl\).

Mesmo com o segundo modelo sendo mais complexo (maior \(df\), porque agora há um spline em \(\sigma\)), o AIC caiu bastante:

## GAMLSS-RS iteration 1: Global Deviance = 28092.9 
## GAMLSS-RS iteration 2: Global Deviance = 28030.24 
## GAMLSS-RS iteration 3: Global Deviance = 28028.11 
## GAMLSS-RS iteration 4: Global Deviance = 28027.64 
## GAMLSS-RS iteration 5: Global Deviance = 28027.53 
## GAMLSS-RS iteration 6: Global Deviance = 28027.5 
## GAMLSS-RS iteration 7: Global Deviance = 28027.5 
## GAMLSS-RS iteration 8: Global Deviance = 28027.5 
## GAMLSS-RS iteration 9: Global Deviance = 28027.5

##           df      AIC
## m2 12.003086 28051.50
## m1  4.919832 28088.53

(4) Avaliando \(\nu(Fl)\)

Aqui, mantemos a estrutura já escolhida para \(\mu(Fl)\) e \(\sigma(Fl)\) e testamos se vale a pena avaliar o parâmetro \(\nu\) da distribuição BCCG. Aqui, \(\nu(Fl)\) passa a ser uma função suave de \(Fl\): isso significa que estamos testando se a forma associada ao parâmetro \(\nu\)) muda com a área \(Fl\).

O modelo 2 continua o menor AIC, assim não há ganho em deixar \(\nu\) variar com \(Fl\). Além disso, note que o modelo 3 é um pouco mais complexo, e essa complexidade não foi compensada por melhora no ajuste.

## GAMLSS-RS iteration 1: Global Deviance = 28091.46 
## GAMLSS-RS iteration 2: Global Deviance = 28028.84 
## GAMLSS-RS iteration 3: Global Deviance = 28027.09 
## GAMLSS-RS iteration 4: Global Deviance = 28026.75 
## GAMLSS-RS iteration 5: Global Deviance = 28026.66 
## GAMLSS-RS iteration 6: Global Deviance = 28026.64 
## GAMLSS-RS iteration 7: Global Deviance = 28026.63 
## GAMLSS-RS iteration 8: Global Deviance = 28026.63 
## GAMLSS-RS iteration 9: Global Deviance = 28026.62

##           df      AIC
## m2 12.003086 28051.50
## m3 12.686807 28052.00
## m1  4.919832 28088.53

Resíduos quantílicos vs covariável

O gráfico mostra os resíduos quantílicos randomizados em função da covariável \(Fl\). A ideia aqui é verificar se ainda existe algum padrão nos resíduos ao longo de \(Fl\).

• Os resíduos estão distribuídos em torno de 0 ao longo de todo o intervalo de \(Fl\);

• O modelo conseguiu capturar a variação de dispersão (especialmente quando \(\sigma\) é modelada como função de \(Fl\));

• Existem alguns valores extremos, mas eles parecem pontuais, não concentrados em uma região específica de \(Fl\).

Diagnóstico do modelo escolhido

Idealmente, os pontos devem estar distribuídos ao redor de zero dentro das bandas.

Neste gráfico, temos um padrão:

• proximidade de zero na região central, mas

• apresenta desvios nas caudas.

Indica que, apesar do bom ajuste global, pode ocorrer:

• leve sub ou superestimação de extremos;

• avaliar parâmetros de forma, por exemplo, \(\nu\) ou \(\tau\)).

Conclusão: o diagnóstico aponta boa adequação global, com indícios de possíveis discrepâncias pequenas nas caudas. A decisão de complexificar o modelo deve considerar se esses desvios nas caudas são relevantes para o objetivo do estudo (por exemplo, previsão de quantis extremos).

Conclusões

• O GAMLSS permite modelar vários parâmetros da distribuição, não apenas a média: \[ (\mu,\sigma,\nu,\tau). \]

• Isso é útil em casos de heterocedasticidade, assimetria ao longo de \(X\), caudas mais leves/pesadas dependendo de \(X\).

• A seleção do modelo deve ser parcimoniosa;

• Validação é indispensável: resíduos quantílicos randomizados; worm plot para avaliar inadequações locais (centro vs caudas).

• o GAMLSS oferece grande flexibilidade, mas exige escolher bem \(D\), evitar complexidade desnecessária e confirmar adequação com diagnóstico.

Alguns materiais

• Ômega Talks com Ricardo Nakamura: Ômega Talks | Ricardo Nakamura · Além da regressão à média: GAMLSS - teoria e aplicações

• Palestra Fernanda de Bastiani: X Bienal da Sociedade Brasileira de Matemática - Palestra Plenária 4: Fernanda de Bastiani

• Site do grupo de pesquisa do GAMLSS: http://www.gamlss.com

Referências