Como escolher o tipo de parâmetros GAM

Comecei a trabalhar com o GAM em R e adquiri o excelente livro de Simon Wood sobre o tópico ( "Modelos aditivos generalizados, uma introdução ao R" ). Com base em um de seus exemplos, estou analisando o seguinte:

library(mgcv)
data(trees)
ct1<-gam(log(Volume) ~ Height + s(Girth), data=trees)

Eu tenho duas perguntas gerais para este exemplo:

1. Como se decide quando uma variável na estimativa do modelo deve ser paramétrica (como Altura) ou quando deve ser suave (como Girth)? Um mantém uma vantagem sobre o outro e existe uma maneira de determinar qual é o tipo ideal para uma variável? Se alguém tiver alguma literatura sobre esse tópico, ficarei feliz em saber.

2. Digamos que eu queira olhar mais de perto os pesos de ct1: ct1$coefficients. Posso usá-los como o gamprocedimento os produz, ou tenho que transformá-los antes de analisá-los, considerando que estou me adaptando log(Volume)? No caso deste último, acho que teria que usarexp (ct1$coefficients)

Suponho que isso seja melhor explicado na 2ª edição do livro de Simon (que deve sair em alguns dias), pois ele e seus alunos só trabalharam parte da teoria nos últimos anos depois que Simon escreveu seu livro.

O que Marra & Wood (2011) mostrou foi que, se queremos fazer a seleção em um modelo com termos suaves, uma abordagem muito boa é adicionar uma penalidade extra a todos os termos suaves. Essa penalidade adicional trabalha com a penalidade de suavidade para esse termo para controlar a ondulação do termo e se um termo deve estar no modelo.

Portanto, a menos que você tenha uma boa teoria para assumir formas / efeitos suaves ou lineares / paramétricos para as covariáveis, você pode abordar o problema como escolhendo entre todos os modelos (representáveis ​​pela combinação aditiva de combinações lineares das funções básicas) entre um com suaviza cada covariável desde o modelo que contém apenas uma interceptação.

Por exemplo:

library(mgcv)
data(trees)
ct1 <- gam(log(Volume) ~ s(Height) + s(Girth), data=trees, method = "REML", select = TRUE)

> summary(ct1)

Family: gaussian 
Link function: identity 

Formula:
log(Volume) ~ s(Height) + s(Girth)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.27273    0.01492   219.3   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
            edf Ref.df      F  p-value    
s(Height) 0.967      9  3.249 3.51e-06 ***
s(Girth)  2.725      9 75.470  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.975   Deviance explained = 97.8%
-REML = -23.681  Scale est. = 0.0069012  n = 31

Observando o resultado (especificamente na seção Coeficientes paramétricos ), observamos que ambos os termos são altamente significativos. Mas observe os graus efetivos de valor da liberdade para a suavidade de Height; é ~ 1. O que esses testes estão fazendo é explicado em Wood (2013).

Isso sugere para mim que Heightdeve entrar no modelo como um termo paramétrico linear. Podemos avaliar isso plotando o encaixe suave:

> plot(ct1, select = 1, shade = TRUE, scale = 0, seWithMean = TRUE)

que dá:

Isso mostra claramente que a forma selecionada do efeito de Heighté linear.

No entanto, se você não soubesse disso antecipadamente (e não o conhecesse, caso contrário não teria feito a pergunta), agora não poderá reequipar o modelo para esses dados usando apenas um termo linear para Height. Isso lhe causaria problemas reais de inferência na linha. A saída em summary()foi responsável pelo fato de você ter feito essa seleção. Se você reequipar o modelo com um efeito paramétrico linear de Height, a saída não saberia disso e você obteria valores p excessivamente otimistas.

Quanto à questão 2, como já mencionado nos comentários, não, não exponencie os coeficientes desse modelo. Além disso, não mergulhe nos modelos ajustados, pois o conteúdo desses componentes nem sempre é o que você poderia esperar. Use as funções do extrator; neste caso coef().

Mais adiante, quando Simon chegar aos GLMs e GAMs, você o verá modelar esses dados por meio de um Gamma GLM:

ct1 <- gam(Volume ~ Height + s(Girth), data=trees, method = "REML",
           family = Gamma(link = "log"))

Nesse modelo, como o ajuste está sendo feito na escala do preditor linear (na escala logarítmica), os coeficientes podem ser exponenciados para obter algum efeito parcial, mas é melhor usar predict(ct1, ...., type = "response")para recuperar valores / previsões ajustados no escala da resposta (em m ^ 3).

Marra, G. & Wood, SN Seleção prática de variáveis ​​para modelos de aditivos generalizados. Comput. Estado. Data Anal. 55, 2372- 2387 (2011).

Madeira, SN Em valores de p para componentes lisos de um modelo de aditivo generalizado estendido. Biometrika 100, 221-228 (2013).