Ajustando funções de afilamento com GAMLSS
Eliel Martins
2025-12-10
Bases de dados
Desenvolvida com dados de inventários florestais em povoamentos nativos e plantados.
Pacotes necessários
library(readxl)
library(gamlss)Leitura dos dados
dados<- read_xlsx("Eucalipto_CampoVerde.xlsx")attach(dados) # padronizar o nome de entradahead(dados) # mostrar primiros dadosNesta etapa, foram construídas variáveis com o objetivo de padronizar as medições ao longo do fuste
e viabilizar diferentes abordagens de modelagem (linear e não linear).
y_lm <- di.cm/dap.cm # variável resposta
x <- hi.m/ht.m #
y_nl <- di.cmAjustando modelos:
Modelos de regressão polinomial e modelos não-lineares parametrizados são alternativas flexíveis para capturar curvaturas, pontos de inflexão e assintotas.
Aqui, serão testados cinco modelos com preditores lineares nos parâmetros: polinômios de 3º, 4º e 5º graus, e um modelo tradicionais da área florestal — Hradetzky (1976).
A escolha entre as distribuições foi realizada com base em critérios de informação, particularmente o GAIC (Generalized Akaike Information Criterion), considerando diferentes valores do parâmetro de penalização k (k = 2, k = 3,84 e k = 6,74). O critério GAIC é definido como: GAIC=−2log(L)+k×df
em que L representa a verossimilhança do modelo e df o número de graus de liberdade estimados. Valores menores de GAIC indicam melhor equilíbrio entre qualidade de ajuste e complexidade do modelo.
O pacote gamlss.dist fornece acesso a mais de 100 distribuições de resposta explícitas voltadas para aplicações no modelo GAMLSS e, além disso, facilita a geração de distribuições de variáveis de resposta “implícitas” definidas pelo usuário.
Será realizado ajuste dos modelos com preditores lineares via gamlss utilizando as distribuições TF2, JSU e SST com, além da verificação dos resíduos e avaliação do desempenho dos modelos ajustados.
Polinômio de 3º Grau (Cúbico)
Captura até um ponto de inflexão Modelo paramétrico linear nos coeficientes
\[ \frac{d(h)}{D} = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 \]
onde \[ x = \frac{h}{H}\] Regressão Linear
m.p3_lm <- lm(y_lm ~ x + I(x^2) + I(x^3)) summary(m.p3_lm)Qualidade do ajuste: R² = 0.9743
R² ajustado = 0.9742
Erro padrão residual = 0.0513
Todos os termos são altamente significativos (p < 2e-16).
análise dos resíduos
summary(m.p3_lm)
plot(m.p3_lm)
wp(m.p3_lm) # nao mostra bom ajuste
Ajustando com gamlss familia NO
Captura Média e Variância
\[Y∼N(μ(x),σ(x))\]
## polinomio de 3º grau
m.p3 <- gamlss(y_lm ~ x + I(x^2) + I(x^3),
sigma.formula = ~ x + I(x^2) + I(x^3),
family = NO)summary(m.p3) #Assimetria = 0.719 (moderada); Curtose = 4.41 (cauda mais pesada que Normal)
plot(m.p3)
wp(m.p3) # não indicou bom ajuste
Equação do modelo ajustado por distribuição normal
Média (μ): μ = 1,1085 - 1,2208·x + 0,9822·x² - 0,9121·x³ Desvio Padrão: σ = exp(-2,4553 - 4,4405·x + 8,6931·x² - 5,3509·x³)
testando as distribuições selecionadas TF2, JSU, SST
dist_m.p3 <- chooseDist(m.p3, type="realAll") ; dist_m.p3 # definindo as distribuições
best_dist_m.p3 <- dist_m.p3[c("TF2","JSU","SST"),] # selecionar dentre as três
min(best_dist_m.p3[,1]) # p/ k=2.00 #melhor dist ? JSU
min(best_dist_m.p3[,2]) # p/ k=3.84 #melhor dist ? JSU
min(best_dist_m.p3[,3]) # p/ k=6.74 #melhor dist ? JSUbest_dist_m.p3 # a melhor AIC dentre as 3 selecionadas: JSUAjustando com JSU - distribuição de quatro parâmetros
m.p3_1 <- gamlss(y_lm ~ x + I(x^2) + I(x^3),
sigma.formula = ~ x + I(x^2) + I(x^3),
family = JSU)
summary(m.p3_1)
plot(m.p3_1)
wp(m.p3_1) # resíduos no intervalo de confianças
ν = 1.4259 (p = 0.006) τ = 0.8139 (p < 0.001)
Comparação entre os modelos:
| Modelo | Global Deviance | AIC |
|---|---|---|
| Normal (NO) | -2729.324 | -2713.324 |
| Johnson SU | -2801.79 | -2781.79 |
A distribuição Johnson SU (JSU), mostrou-se superior à distribuição Normal (NO) para todos os critérios de informação avaliados. A especificação do modelo é dada por:
μ = 1,105 - 1,217·x + 0,998·x² - 0,927·x³ log(σ) = -2,677 - 3,148·x + 6,991·x² - 4,620·x³ ν = 1,426 (assimetria constante) τ = 2,256 (curtose constante)
A significância estatística dos parâmetros ν (p = 0,006) e τ (p < 0,001) confirma a presença de assimetria e curtose não-normais nos dados, justificando a escolha da distribuição JSU. A análise dos resíduos quantílicos indica excelente ajuste, com estatísticas muito próximas dos valores teóricos ideais (média = 0,002; variância = 1,001; assimetria = 0,024; curtose = 3,000).
modelo: μ = 1,105402 - 1,216621·x + 0,998129·x² - 0,927394·x³ Parâmetro de Escala (σ - Desvio Padrão)
log(σ) = -2,67717 - 3,14763·x + 6,99114·x² - 4,61954·x³ σ = exp(-2,67717 - 3,14763·x + 6,99114·x² - 4,61954·x³) Parâmetro de Assimetria (ν - Nu)
ν = 1,4259 (constante) Parâmetro de Curtose (τ - Tau)
log(τ) = 0,8139 τ = exp(0,8139) = 2,256 (constante)
Modelo de Prodan - polinomio de 4º grau
\[ μ(x)=β 0 +β 1 x+β 2 x^2 +β 3 x^3 +β 4 x^4 \]
Ajustando com distribuição NO
mod_pol_4 <- gamlss(y_lm ~ x + I(x^2) + I(x^3) + I(x^4),
sigma.formula = ~ x + I(x^2) + I(x^3) + I(x^4),
family = NO)verificando resíduos e sifnificância
summary(mod_pol_4)
plot(mod_pol_4) # resíduos assimétricos
wp(mod_pol_4) # Melhorar a modelagem
Coeficiente de Assimetria = 0.773, indica:Distribuição assimétrica à direita Presença de mais resíduos positivos extremos
Curtose = 4.45, Valor> 3 - leptocurtose, caudas mais pesadas que a Normal
Filliben Correlation = 0.9837: indica pequenos desvios da normalidade, especialmente nas caudas.
Selecionando as distribuições
dist_mod_pol_4 <- chooseDist(mod_pol_4, type="realAll") ; dist_mod_pol_4
best_dist_mod_pol_4 <- dist_mod_pol_4[c("TF2","JSU","SST"),] ; best_dist_mod_pol_4
best_dist_mod_pol_4 <- na.omit(best_dist_mod_pol_4) ; best_dist_mod_pol_4
min(best_dist_mod_pol_4[,1]) # p/ k=2.00 #melhor dist ? SST
min(best_dist_mod_pol_4[,2]) # p/ k=3.84 #melhor dist ? SST
min(best_dist_mod_pol_4[,3]) # p/ k=6.74 #melhor dist ? SST2 3.84 6.74 TF2 -2778.796 -2758.556 -2726.656 JSU -2827.089 -2805.009 -2770.209 SST -2828.296 -2806.216 -2771.416
A distribuição SST forneceu o menor AIC para o modelo
\[ Y∼SST(μ(x),σ(x),ν,τ) \]
mod_pol_4_1 <- gamlss(y_lm ~ x + I(x^2) + I(x^3) + I(x^4),
sigma.formula = ~ x + I(x^2) + I(x^3) + I(x^4),
family = SST
)
summary(mod_pol_4_1)
plot(mod_pol_4_1) # bom ajuste 
wp(mod_pol_4_1) # dentro do intervalo de confiança
\[ μ(x)=β 0 +β 1 x+β 2 x^2 +β3x^3 +β 4x^4 \] \[log(σ(x))=γ 0 +γ 1 x+γ 2 x^2 +γ 3 x^3 +γ4 x^4 \] \[log(ν)=δ 0 \] \[ τ=constante (implícita)\]
Modelo
\[μ(x) =1.1268−1.68995x+3.22563x^2 −4.49433x^3 +1.82285x^4\] \[log(σ(x))=−2.7849−0.1013x−9.5413x^2 +24.7893x^3 −16.4381x^4\] \[log(ν)=0.39333\] \[ τ =≈constante (não //estimada) \]
Polinômio de 5º grau:
\[μ(x)=β 0 +β 1 x+β 2 x^2 +β 3 x^3 +β 4 x^4 +β 5 x^5\]
Ajustando com distbuição - family: Normal
mod_pol_5 <- gamlss(y_lm ~ x + I(x^2) + I(x^3) + I(x^4)+ I(x^5),
sigma.formula = ~ x + I(x^2) + I(x^3) + I(x^4)+ I(x^5),
family = NO)
summary(mod_pol_5) # não capturou sigma
plot(mod_pol_5) # resíduos assimétricos
wp(mod_pol_5) # precisa melhorar o ajuste dos resíduos
x⁴ e x⁵ não contribuem para a média (sem sifnificância)
σ: heterocedasticidade não foi capturada (somente β0 significativo)
selecionando a distribuição com GAMLSS
dist_mod_pol_5 <- chooseDist(mod_pol_5, type="realAll") ; dist_mod_pol_5
best_dist_mod_pol_5 <- dist_mod_pol_5[c("TF2", "JSU","SST"),] ; best_dist_mod_pol_5
min(best_dist_mod_pol_5[,1]) # p/ k=2.00 #melhor dist SST
min(best_dist_mod_pol_5[,2]) # p/ k=3.84 #melhor dist SST
min(best_dist_mod_pol_5[,3]) # p/ k=6.74 #melhor dist SSTbest_dist_mod_pol_5SST apresentou mellhor distribuição
2 3.84 6.74 TF2 -2778.640 -2754.720 -2717.020 JSU -2826.835 -2801.075 -2760.475 SST -2827.995 -2802.235 -2761.635
Ajustando com a distribuição SST
mod_pol_5_1 <- gamlss(y_lm ~ x + I(x^2) + I(x^3) + I(x^4)+ I(x^5),
sigma.formula = ~ x + I(x^2) + I(x^3) + I(x^4) + I(x^5),
family = SST)
summary(mod_pol_5_1)
plot(mod_pol_5_1)
wp(mod_pol_5_1) # adequado
μ(x)= coeficientes estatisticamente significativos (p < 0.001, exceto x⁵ que ainda é significativo p = 0.000115).
Sinais alternados (- + - + -) indicam múltiplas mudanças de concavidade ao longo da altura: Permite capturar detalhes finos do afilamento do tronco. o polinômio de 5º grau é capaz de modelar com precisão pequenas irregularidades na forma do fuste que não seriam capturadas com polinômios de grau menor
log(σ(x))= apenas o termo x⁵ é significativo (p = 0.0386), x e x⁴ têm significância marginal (p < 0.1).
Isso indica que a variabilidade muda ao longo da altura, mas de forma complexa e não linear.
Assimetria (ν) = exp(0.39461) ≈ 1.484, significativa (p < 0.001)
Indica assimetria positiva, ou seja, a cauda direita é mais longa
Modelo de Hradetsky
O modelo de Hradetsky é uma generalização polinomial contínua, mas em vez de usar apenas potências inteiras (x², x³…), ele usa uma série de potências fracionárias ou pequenas sobre a variável independente x.
A ideia \[ y=f(x)+ϵ\] com: \[ f(x)= i ∑ βix^p i \]
o modelo tenta capturar padrões não-lineares complexos em x usando várias potências fracionárias e inteiras.
Ajustando com a distribuição Normal NO
pot1 <- seq(0.005, 0.01, by = 0.001) ; pot1 # pequenos expoentes <0.01
pot2 <- seq(0.01, 0.1, by = 0.01) ; pot2 # pequenos expoentes 0.01-0.1
pot3 <- seq(0.1, 1, by = 0.1) ; pot3 # 0.1-1
pot4 <- seq(1, 30, by = 1) ; pot4 # grandes expoentes 1-30
m.hrad <- gamlss(y_lm ~ I(x^0.005) + I(x^pot1[2]) + I(x^pot1[3])+ I(x^pot1[4]) + I(x^pot1[5]) + I(x^pot1[6]) +
I(x^pot2[1]) + I(x^pot2[2]) + I(x^pot2[3])+ I(x^pot2[4]) + I(x^pot2[5]) + I(x^pot2[6]) + I(x^pot2[7]) + I(x^pot2[8])+ I(x^pot2[9]) + I(x^pot2[10]) +
I(x^pot3[1]) + I(x^pot3[2]) + I(x^pot3[3])+ I(x^pot3[4]) + I(x^pot3[5]) + I(x^pot3[6]) + I(x^pot3[7]) + I(x^pot3[8])+ I(x^pot3[9]) + I(x^pot3[10]) +
I(x^pot4[1]) + I(x^pot4[2]) + I(x^pot4[3])+ I(x^pot4[4]) + I(x^pot4[5]) + I(x^pot4[6]) + I(x^pot4[7]) + I(x^pot4[8])+ I(x^pot4[9]) + I(x^pot4[10]) +
I(x^pot4[11]) + I(x^pot4[12]) + I(x^pot4[13])+ I(x^pot4[14]) + I(x^pot4[15]) + I(x^pot4[16]) + I(x^pot4[17]) + I(x^pot4[18])+ I(x^pot4[19]) + I(x^pot4[20]) +
I(x^pot4[21]) + I(x^pot4[22]) + I(x^pot4[23])+ I(x^pot4[24]) + I(x^pot4[25]) + I(x^pot4[26]) + I(x^pot4[27]) + I(x^pot4[28])+ I(x^pot4[29]) + I(x^pot4[30])
#sigma.formula = ~ I(x^0.005) + I(x^pot1[2]) + I(x^pot1[3])+ I(x^pot1[4]) + I(x^pot1[5]) + I(x^pot1[6]) +
#I(x^pot2[1]) + I(x^pot2[2]) + I(x^pot2[3])+ I(x^pot2[4]) + I(x^pot2[5]) + I(x^pot2[6]) + I(x^pot2[7]) + I(x^pot2[8])+ I(x^pot2[9]) + I(x^pot2[10]) +
#I(x^pot3[1]) + I(x^pot3[2]) + I(x^pot3[3])+ I(x^pot3[4]) + I(x^pot3[5]) + I(x^pot3[6]) + I(x^pot3[7]) + I(x^pot3[8])+ I(x^pot3[9]) + I(x^pot3[10]) +
#I(x^pot4[1]) + I(x^pot4[2]) + I(x^pot4[3])+ I(x^pot4[4]) + I(x^pot4[5]) + I(x^pot4[6]) + I(x^pot4[7]) + I(x^pot4[8])+ I(x^pot4[9]) + I(x^pot4[10]) +
#I(x^pot4[11]) + I(x^pot4[12]) + I(x^pot4[13])+ I(x^pot4[14]) + I(x^pot4[15]) + I(x^pot4[16]) + I(x^pot4[17]) + I(x^pot4[18])+ I(x^pot4[19]) + I(x^pot4[20]) +
#I(x^pot4[21]) + I(x^pot4[22]) + I(x^pot4[23])+ I(x^pot4[24]) + I(x^pot4[25]) + I(x^pot4[26]) + I(x^pot4[27]) + I(x^pot4[28])+ I(x^pot4[29]) + I(x^pot4[30]),
)
summary(m.hrad)
plot(m.hrad) # assimetria
wp(m.hrad) # não estimou bem os residuos 
media = 4.448e-12
variancia = 1.001183 → aproximadamente 1, consistente com Normal
v = 1.294, bastante assimetria
t = 7.363, alta curtose, presença de outliers ou caudas pesadas
Filliben correlation = 0.9658, qq-plot ainda razoável, mas menor que modelos mais simples
O modelo não conseguiu capturar bem a forma da distribuição residual.
Há resíduos assimétricos e caudas pesadas, mesmo com tantas potências.
Isso sugere que a distribuição Normal pode não ser adequada para o Hradetsky ou que o modelo está muito instável.
Selecionado a melhor distribuição do modelo com Gamlss
dist_m.hrad <- chooseDist(m.hrad, type="realAll") ; dist_m.hrad
best_dist_m.hrad <- dist_m.hrad[c("TF2","JSU","SST"),] ; best_dist_m.hrad
best_dist_m.hrad <- na.omit(best_dist_m.hrad)
min(best_dist_m.hrad[,1]) # p/ k=2.00 #melhor dist SST
min(best_dist_m.hrad[,2]) # p/ k=3.84 #melhor dist SST
min(best_dist_m.hrad[,3]) # p/ k=6.74 #melhor dist SST best_dist_m.hrad # SSTAjustando com a distribuição SST
m.hrad_1 <- gamlss(y_lm ~ I(x^0.005) + I(x^pot1[2]) + I(x^pot1[3])+ I(x^pot1[4]) + I(x^pot1[5]) + I(x^pot1[6]) +
I(x^pot2[1]) + I(x^pot2[2]) + I(x^pot2[3])+ I(x^pot2[4]) + I(x^pot2[5]) + I(x^pot2[6]) + I(x^pot2[7]) + I(x^pot2[8])+ I(x^pot2[9]) + I(x^pot2[10]) +
I(x^pot3[1]) + I(x^pot3[2]) + I(x^pot3[3])+ I(x^pot3[4]) + I(x^pot3[5]) + I(x^pot3[6]) + I(x^pot3[7]) + I(x^pot3[8])+ I(x^pot3[9]) + I(x^pot3[10]) +
I(x^pot4[1]) + I(x^pot4[2]) + I(x^pot4[3])+ I(x^pot4[4]) + I(x^pot4[5]) + I(x^pot4[6]) + I(x^pot4[7]) + I(x^pot4[8])+ I(x^pot4[9]) + I(x^pot4[10]) +
I(x^pot4[11]) + I(x^pot4[12]) + I(x^pot4[13])+ I(x^pot4[14]) + I(x^pot4[15]) + I(x^pot4[16]) + I(x^pot4[17]) + I(x^pot4[18])+ I(x^pot4[19]) + I(x^pot4[20]) +
I(x^pot4[21]) + I(x^pot4[22]) + I(x^pot4[23])+ I(x^pot4[24]) + I(x^pot4[25]) + I(x^pot4[26]) + I(x^pot4[27]) + I(x^pot4[28])+ I(x^pot4[29]) + I(x^pot4[30]),
#sigma.formula = ~ I(x^0.005) + I(x^pot1[2]) + I(x^pot1[3])+ I(x^pot1[4]) + I(x^pot1[5]) + I(x^pot1[6]) +
#I(x^pot2[1]) + I(x^pot2[2]) + I(x^pot2[3])+ I(x^pot2[4]) + I(x^pot2[5]) + I(x^pot2[6]) + I(x^pot2[7]) + I(x^pot2[8])+ I(x^pot2[9]) + I(x^pot2[10]) +
#I(x^pot3[1]) + I(x^pot3[2]) + I(x^pot3[3])+ I(x^pot3[4]) + I(x^pot3[5]) + I(x^pot3[6]) + I(x^pot3[7]) + I(x^pot3[8])+ I(x^pot3[9]) + I(x^pot3[10]) +
#I(x^pot4[1]) + I(x^pot4[2]) + I(x^pot4[3])+ I(x^pot4[4]) + I(x^pot4[5]) + I(x^pot4[6]) + I(x^pot4[7]) + I(x^pot4[8])+ I(x^pot4[9]) + I(x^pot4[10]) +
#I(x^pot4[11]) + I(x^pot4[12]) + I(x^pot4[13])+ I(x^pot4[14]) + I(x^pot4[15]) + I(x^pot4[16]) + I(x^pot4[17]) + I(x^pot4[18])+ I(x^pot4[19]) + I(x^pot4[20]) +
#I(x^pot4[21]) + I(x^pot4[22]) + I(x^pot4[23])+ I(x^pot4[24]) + I(x^pot4[25]) + I(x^pot4[26]) + I(x^pot4[27]) + I(x^pot4[28])+ I(x^pot4[29]) + I(x^pot4[30]),
family = SST)
summary(m.hrad_1) #apenas intercepto significativo
plot(m.hrad_1)
wp(m.hrad_1) # melhorou, resíduos abaixo do intervalo de confiança
| Parâmetro | Interpretação | Valor estimado | Significância |
|---|---|---|---|
| Sigma | dispersão da Normal | -2.985 | *** |
| Nu | assimetria | 0.429 | *** |
| Tau | caudas pesadas | 0.993 | *** |
significativos estatisticamente
mean = -0.0026: bom
variance = 1.0069: quase 1, consistente
skewness = 0.05: resíduos quase simétricos
kurtosis = 2.94: levemente leptocúrticos, perto de Normal
Filliben correlation = 0.99897: qq-plot praticamente perfeito Os resíduos estão aproximadamente normais, simétricos e sem outliers extremos.
O modelo Hradetsky ajustado com a distribuição SST ajusta a relação complexa entre y e x muito melhor do que polinômios de grau 5 ou Hradetsky com Normal. Ele lida com assimetria e caudas pesadas, resultando em resíduos quase normais e ajuste confiável.