Ajuste gamlss para a função de afilamento kozak (2004)
Eliel Martins
2026-04-17
Carregar pacotes
library(dplyr)
library(readxl)
library(gamlss)
library(gamlss.nl)Lendo dados
dados <- read_excel("dados_gamlss.xlsx")attach(dados)Visualização
head(dados)## # A tibble: 6 × 7
## cod tree d h hi di i
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 401833MUT 40 6.5 6.5 0.1 9.6 3
## 2 401833MUT 40 6.5 6.5 0.5 8 3
## 3 401833MUT 40 6.5 6.5 1 6.8 3
## 4 401833MUT 40 6.5 6.5 1.3 6.5 3
## 5 401833MUT 40 6.5 6.5 2 5.8 3
## 6 401833MUT 40 6.5 6.5 3 4.8 3dispersão dos dados
plot(di, hi/h)
Definindo a função do modelo Kozak 2004
Kozak_2004 <- function(hi, h, d, alpha0, alpha1, alpha2, beta1, beta2, beta3, beta4, beta5, beta6) {
m <- (1 - (hi/h)^(1/3)) / (1 - ((1/3)/h)^(1/3))
C <- beta1 * (hi/h)^4 + beta2 * (1/exp(d/h)) + beta3 * m^0.1 + beta4 * (1/d) + beta5 * h ^ (1 - (hi/h)^(1/3)) + beta6 * m
di <- (alpha0 * d^(alpha1) * h^(alpha2)) * m^C
return(di)
}ajustando por nls funão não linear
nls_model <- nls(di ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2, beta1, beta2, beta3, beta4, beta5, beta6),
start = list(alpha0 = 0.2, alpha1 = 0.5, alpha2 = 0.3, beta1 = 0.4, beta2 = 0.2, beta3 = 0.5, beta4 = 0.8, beta5 = 0.4, beta6 = 0.6))relatorio do modelo ajustado b2 e b6 não significativos para o modelo
summary(nls_model)##
## Formula: di ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2, beta1, beta2,
## beta3, beta4, beta5, beta6)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## alpha0 1.6644638 0.0236972 70.239 < 2e-16 ***
## alpha1 0.8833136 0.0115060 76.770 < 2e-16 ***
## alpha2 0.0009896 0.0133084 0.074 0.94073
## beta1 0.3915819 0.0184098 21.270 < 2e-16 ***
## beta2 0.0218395 0.0548372 0.398 0.69046
## beta3 0.5546203 0.0177143 31.309 < 2e-16 ***
## beta4 0.4225621 0.1421017 2.974 0.00296 **
## beta5 0.0119973 0.0061979 1.936 0.05298 .
## beta6 0.0219612 0.0380047 0.578 0.56340
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7963 on 3973 degrees of freedom
##
## Number of iterations to convergence: 8
## Achieved convergence tolerance: 3.53e-06salvando coeficientes
ib <- coef(nls_model)
coef(nls_model)## alpha0 alpha1 alpha2 beta1 beta2 beta3
## 1.6644637848 0.8833136057 0.0009896338 0.3915818541 0.0218395077 0.5546202987
## beta4 beta5 beta6
## 0.4225621204 0.0119972721 0.0219611712transformando hi/h em hi_h (facilitar)
hi_h <- hi/hAjuste do Modelo com GAMLSS assumindo uma distribuição Normal
m0 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8))visualizar summary e residuos
plot(m0)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.005684386
## variance = 1.000219
## coef. of skewness = 0.04214068
## coef. of kurtosis = 5.512621
## Filliben correlation coefficient = 0.9904048
## ******************************************************************wp(m0)
summary(m0)## *******************************************************************
## Family: c("NO", "Normal")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~I(sqrt(hi_h)) +
## hi_h, mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01, -1,
## 0.01), control = NL.control(iterlim = 1000, ndigit = 15,
## steptol = 1e-08))
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.664584 0.025206 66.0395 0.000e+00
## alpha1 0.902815 0.011934 75.6523 0.000e+00
## alpha2 -0.020144 0.013880 -1.4513 1.467e-01
## beta1 0.427352 0.022063 19.3692 1.404e-83
## beta2 -0.017668 0.058131 -0.3039 7.612e-01
## beta3 0.557968 0.018700 29.8373 1.284e-195
## beta4 0.433910 0.141696 3.0623 2.197e-03
## beta5 0.007853 0.006646 1.1817 2.373e-01
## beta6 0.057393 0.039784 1.4426 1.491e-01
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.2503 0.0445 5.626 1.847e-08
## I(sqrt(hi_h)) -3.0362 0.1993 -15.234 2.119e-52
## hi_h 3.1760 0.1954 16.257 1.999e-59
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 12
## Residual Deg. of Freedom: 3970
## at cycle: 61
##
## Global Deviance: 9194.985
## AIC: 9218.985
## SBC: 9294.459
## *******************************************************************o modelo de distribuição normal não a heterogeneidade dos dados com caudas pesadas
Summary of the Quantile Residuals
mean = 0.005684386
variance = 1.000219
coef. of skewness = 0.04214068
coef. of kurtosis = 5.512621 Filliben correlation coefficient = 0.9904048 ******************************************************************
Modelo kozack com distribuição logistica LO
m1 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-0.01,0.01),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = LO())visualizar summary e residuos
plot(m1)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.01086794
## variance = 1.007958
## coef. of skewness = 0.03426083
## coef. of kurtosis = 3.327807
## Filliben correlation coefficient = 0.9991151
## ******************************************************************wp(m1)
summary(m1)## *******************************************************************
## Family: c("LO", "Logistic")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~I(sqrt(hi_h)) +
## hi_h, mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01, -0.01,
## 0.01), control = NL.control(iterlim = 1000, ndigit = 15,
## steptol = 1e-08), family = LO())
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.656591 0.023502 70.48803 0.000e+00
## alpha1 0.915145 0.011012 83.10088 0.000e+00
## alpha2 -0.032171 0.012708 -2.53165 1.135e-02
## beta1 0.437059 0.021268 20.55023 7.661e-94
## beta2 -0.002004 0.055145 -0.03635 9.710e-01
## beta3 0.550085 0.017910 30.71385 3.719e-207
## beta4 0.372620 0.131193 2.84025 4.508e-03
## beta5 0.002370 0.006454 0.36719 7.135e-01
## beta6 0.097638 0.037665 2.59226 9.535e-03
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) -0.3198 0.05415 -5.906 3.516e-09
## I(sqrt(hi_h)) -3.4257 0.24298 -14.099 3.874e-45
## hi_h 3.6540 0.23715 15.408 1.453e-53
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 12
## Residual Deg. of Freedom: 3970
## at cycle: 50
##
## Global Deviance: 8982.139
## AIC: 9006.139
## SBC: 9081.614
## *******************************************************************Comparação entre AIC dos Modelos
AIC(m0,m1)## df AIC
## m1 12 9006.139
## m0 12 9218.985Modelos GAMLSS de 3 Parametros nu constante = 1
Distribuição TF2
m2 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-0.01,0.01),
nu.formula = ~1,
nu.start = c(0.1),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = TF2())Resíduos
plot(m2)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.01112568
## variance = 1.000144
## coef. of skewness = 0.03129896
## coef. of kurtosis = 3.0146
## Filliben correlation coefficient = 0.9996696
## ******************************************************************wp(m2)
summary(m2)## *******************************************************************
## Family: c("TF2", "t Family 2")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~I(sqrt(hi_h)) +
## hi_h, mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01, -0.01,
## 0.01), nu.formula = ~1, nu.start = c(0.1), control = NL.control(iterlim = 1000,
## ndigit = 15, steptol = 1e-08), family = TF2())
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.656090 0.023379 70.83637 0.000e+00
## alpha1 0.916030 0.010945 83.69118 0.000e+00
## alpha2 -0.033051 0.012617 -2.61962 8.803e-03
## beta1 0.439486 0.021248 20.68332 4.895e-95
## beta2 0.002804 0.054843 0.05113 9.592e-01
## beta3 0.548191 0.017848 30.71451 3.644e-207
## beta4 0.363518 0.129874 2.79901 5.126e-03
## beta5 0.001556 0.006459 0.24093 8.096e-01
## beta6 0.103524 0.037571 2.75540 5.862e-03
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.3043 0.05702 5.336 9.519e-08
## I(sqrt(hi_h)) -3.5347 0.25271 -13.987 1.868e-44
## hi_h 3.7833 0.24674 15.333 4.598e-53
##
## -------------------------------------------------------------------
## Nu link function: logshiftto2
## Nu Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 1.423 0.1407 10.11 4.784e-24
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 13
## Residual Deg. of Freedom: 3969
## at cycle: 154
##
## Global Deviance: 8967.472
## AIC: 8993.472
## SBC: 9075.236
## *******************************************************************Distribuição PE
m3 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
nu.formula = ~ 1,
nu.start = c(0.1),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = PE())Residuos e summary do modelo
plot(m3)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.02194838
## variance = 0.9975771
## coef. of skewness = 0.01776158
## coef. of kurtosis = 3.338591
## Filliben correlation coefficient = 0.9987243
## ******************************************************************wp(m3)
summary(m3)## *******************************************************************
## Family: c("PE", "Power Exponential")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~I(sqrt(hi_h)) +
## hi_h, mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01, -1,
## 0.01), nu.formula = ~1, nu.start = c(0.1), control = NL.control(iterlim = 1000,
## ndigit = 15, steptol = 1e-08), family = PE())
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.6475141 NaN NaN NaN
## alpha1 0.9277981 0.006322 146.7656 0.000e+00
## alpha2 -0.0440492 0.007875 -5.5938 2.221e-08
## beta1 0.4487972 0.013926 32.2278 7.205e-228
## beta2 -0.0221470 0.017929 -1.2353 2.167e-01
## beta3 0.5527938 0.002923 189.1365 0.000e+00
## beta4 0.3352918 0.096315 3.4812 4.991e-04
## beta5 -0.0002936 0.001416 -0.2073 8.358e-01
## beta6 0.1242208 0.004778 25.9979 5.231e-149
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.2832 0.05545 5.107 3.282e-07
## I(sqrt(hi_h)) -3.3801 0.24863 -13.595 4.296e-42
## hi_h 3.6016 0.24385 14.770 2.299e-49
##
## -------------------------------------------------------------------
## Nu link function: log
## Nu Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.2816 0.03122 9.019 1.892e-19
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 13
## Residual Deg. of Freedom: 3969
## at cycle: 158
##
## Global Deviance: 9004.349
## AIC: 9030.349
## SBC: 9112.113
## *******************************************************************Aic dos modelos de 3 parametros
AIC(m2,m3)## df AIC
## m2 13 8993.472
## m3 13 9030.349#Modelos GAMlSS com 4 parâmetros nu e tau constantes
Distribuição BCT
m4 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
nu.formula = ~ 1,
tau.formula = ~1,
nu.start = c(0.1),
tau.start = c(0.1),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = BCT())Resíduos e coeficientes
plot(m4)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.007393218
## variance = 0.9995241
## coef. of skewness = 0.118508
## coef. of kurtosis = 3.01934
## Filliben correlation coefficient = 0.999348
## ******************************************************************wp(m4)
summary(m4)## *******************************************************************
## Family: c("BCT", "Box-Cox t")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~I(sqrt(hi_h)) +
## hi_h, mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01, -1,
## 0.01), nu.formula = ~1, tau.formula = ~1, nu.start = c(0.1),
## tau.start = c(0.1), control = NL.control(iterlim = 1000,
## ndigit = 15, steptol = 1e-08), family = BCT())
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.673111 0.019186 87.2054 0.000e+00
## alpha1 0.911683 0.010925 83.4485 0.000e+00
## alpha2 -0.031885 0.012842 -2.4830 1.303e-02
## beta1 0.410302 0.018411 22.2855 5.107e-110
## beta2 0.087220 0.055607 1.5685 1.168e-01
## beta3 0.534280 0.018513 28.8599 3.808e-183
## beta4 0.199451 0.106014 1.8814 5.992e-02
## beta5 -0.001259 0.006296 -0.1999 8.415e-01
## beta6 0.124451 0.033619 3.7018 2.141e-04
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) -2.815 0.05696 -49.42 0.000e+00
## I(sqrt(hi_h)) -3.403 0.25359 -13.42 4.686e-41
## hi_h 5.247 0.24746 21.20 9.202e-100
##
## -------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 1.031 0.1857 5.554 2.797e-08
##
## -------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 2.035 0.1138 17.88 1.76e-71
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 14
## Residual Deg. of Freedom: 3968
## at cycle: 86
##
## Global Deviance: 8797.073
## AIC: 8825.073
## SBC: 8913.126
## *******************************************************************Distribuição JSU
m5 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
nu.formula = ~ 1,
tau.formula = ~1,
nu.start = c(0.1),
tau.start = c(0.1),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = JSU())Gráfico de resíduos
plot(m5)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.007691703
## variance = 1.000365
## coef. of skewness = -0.004822054
## coef. of kurtosis = 3.044932
## Filliben correlation coefficient = 0.9996314
## ******************************************************************wp(m5)
Distribuição SEP2
m6 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
nu.formula = ~ 1,
tau.formula = ~1,
nu.start = c(0.1),
tau.start = c(0.1),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = SEP2())Gráfico de resíduos
plot(m6)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.001237366
## variance = 0.9968429
## coef. of skewness = -0.1465531
## coef. of kurtosis = 3.36367
## Filliben correlation coefficient = 0.9981808
## ******************************************************************wp(m6)
summary(m6)## *******************************************************************
## Family: c("SEP2", "Skew Exponential Power type 2")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~I(sqrt(hi_h)) +
## hi_h, mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01, -1,
## 0.01), nu.formula = ~1, tau.formula = ~1, nu.start = c(0.1),
## tau.start = c(0.1), control = NL.control(iterlim = 1000,
## ndigit = 15, steptol = 1e-08), family = SEP2())
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.587e+00 0.0003882 4088.78008 0.000e+00
## alpha1 9.488e-01 NaN NaN NaN
## alpha2 -5.851e-02 NaN NaN NaN
## beta1 4.953e-01 0.0076797 64.49851 0.000e+00
## beta2 -3.808e-02 0.0091104 -4.17937 2.923e-05
## beta3 5.763e-01 0.0042294 136.24960 0.000e+00
## beta4 5.082e-01 0.0232000 21.90699 2.228e-106
## beta5 2.915e-05 0.0020904 0.01395 9.889e-01
## beta6 6.913e-02 0.0103893 6.65438 2.845e-11
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.1762 0.05508 3.199 1.378e-03
## I(sqrt(hi_h)) -3.4442 0.24273 -14.190 1.059e-45
## hi_h 3.6984 0.23828 15.522 2.480e-54
##
## -------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.4388 0.0214 20.5 2.112e-93
##
## -------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.269 0.02984 9.016 1.951e-19
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 14
## Residual Deg. of Freedom: 3968
## at cycle: 429
##
## Global Deviance: 8966.827
## AIC: 8994.827
## SBC: 9082.881
## *******************************************************************Distribuição SST
m7 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
nu.formula = ~ 1,
tau.formula = ~1,
nu.start = c(0.1),
tau.start = c(0.1),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = SST())Gráfico de resíduos
plot(m7)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.006137366
## variance = 1.000581
## coef. of skewness = -0.01353557
## coef. of kurtosis = 3.014731
## Filliben correlation coefficient = 0.999706
## ******************************************************************wp(m7)
summary(m7)## *******************************************************************
## Family: c("SST", "SST")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~I(sqrt(hi_h)) +
## hi_h, mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01, -1,
## 0.01), nu.formula = ~1, tau.formula = ~1, nu.start = c(0.1),
## tau.start = c(0.1), control = NL.control(iterlim = 1000,
## ndigit = 15, steptol = 1e-08), family = SST())
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.6649271 0.024188 68.83190 0.000e+00
## alpha1 0.9163405 0.010870 84.29754 0.000e+00
## alpha2 -0.0352987 0.012631 -2.79462 5.196e-03
## beta1 0.4367216 0.021362 20.44379 6.823e-93
## beta2 -0.0025682 0.054860 -0.04681 9.627e-01
## beta3 0.5525817 0.018096 30.53537 8.844e-205
## beta4 0.3531561 0.130111 2.71427 6.642e-03
## beta5 0.0003907 0.006494 0.06017 9.520e-01
## beta6 0.1060993 0.037584 2.82296 4.758e-03
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.3009 0.05704 5.275 1.329e-07
## I(sqrt(hi_h)) -3.5341 0.25279 -13.980 2.057e-44
## hi_h 3.7920 0.24700 15.352 3.421e-53
##
## -------------------------------------------------------------------
## Nu link function: log
## Nu Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.03427 0.02404 1.425 0.154
##
## -------------------------------------------------------------------
## Tau link function: logshiftto2
## Tau Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 1.421 0.1408 10.09 6.023e-24
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 14
## Residual Deg. of Freedom: 3968
## at cycle: 166
##
## Global Deviance: 8965.443
## AIC: 8993.443
## SBC: 9081.496
## *******************************************************************Aic dos modelos de 4 parametro
AIC(m4, m5, m6, m7)## df AIC
## m4 14 8825.073
## m7 14 8993.443
## m6 14 8994.827
## m5 14 8996.624nu e tau em função de hi/h m8
Distribuição SST
m8 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
nu.formula = ~ hi_h,
tau.formula = ~1,
nu.start = c(0.1,0.1),
tau.start = c(0.1),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = SST())Gráfico resíduos
plot(m8)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.0008474134
## variance = 1.000469
## coef. of skewness = 0.01901981
## coef. of kurtosis = 2.993278
## Filliben correlation coefficient = 0.9996119
## ******************************************************************wp(m8)
Distribuição JSU
m9 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ I(sqrt(hi_h))+ hi_h,
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1,0.01),
nu.formula = ~ hi_h,
tau.formula = ~ 1,
nu.start = c(0.01,0.01),
tau.start = c(0.01),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = JSU())Gráfico
plot(m9)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.005982124
## variance = 0.9981927
## coef. of skewness = 0.005124982
## coef. of kurtosis = 3.016811
## Filliben correlation coefficient = 0.9995401
## ******************************************************************wp(m9)
Distribuição BCT , nu e tau hi_h
m10 <- nlgamlss(
di,
mu.formula = ~ Kozak_2004(hi, h, d, alpha0, alpha1, alpha2,
beta1, beta2, beta3, beta4, beta5, beta6),
sigma.formula = ~ pb(hi_h),
mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
beta1 = ib[4], beta2 = ib[5], beta3 = ib[6],
beta4 = ib[7], beta5 = ib[8], beta6 = ib[9]),
sigma.start = c(0.01,-1),
nu.formula = ~ hi_h,
tau.formula = ~ hi_h,
nu.start = c(0.01,0.01),
tau.start = c(0.01,0.01),
control = NL.control(iterlim = 1000, ndigit = 15, steptol = 1e-8),
family = BCT())Resíduos
plot(m10)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = 0.01252643
## variance = 0.9938559
## coef. of skewness = 0.05248493
## coef. of kurtosis = 3.066212
## Filliben correlation coefficient = 0.9994093
## ******************************************************************wp(m10)
summary(m10)## *******************************************************************
## Family: c("BCT", "Box-Cox t")
##
## Call: nlgamlss(di, mu.formula = ~Kozak_2004(hi, h, d, alpha0, alpha1,
## alpha2, beta1, beta2, beta3, beta4, beta5, beta6), sigma.formula = ~pb(hi_h),
## mu.start = c(alpha0 = ib[1], alpha1 = ib[2], alpha2 = ib[3],
## beta1 = ib[4], beta2 = ib[5], beta3 = ib[6], beta4 = ib[7],
## beta5 = ib[8], beta6 = ib[9]), sigma.start = c(0.01,
## -1), nu.formula = ~hi_h, tau.formula = ~hi_h, nu.start = c(0.01,
## 0.01), tau.start = c(0.01, 0.01), control = NL.control(iterlim = 1000,
## ndigit = 15, steptol = 1e-08), family = BCT())
##
## Fitting method: "JL()"
##
## -------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t-value p-value
## alpha0 1.677747 0.017421 96.303 0.000e+00
## alpha1 0.906268 0.009866 91.854 0.000e+00
## alpha2 -0.027843 0.011440 -2.434 1.494e-02
## beta1 0.362120 0.017168 21.093 9.252e-99
## beta2 0.154235 0.053757 2.869 4.116e-03
## beta3 0.528194 0.017703 29.836 1.333e-195
## beta4 0.160861 0.103697 1.551 1.208e-01
## beta5 -0.006437 0.005628 -1.144 2.527e-01
## beta6 0.129411 0.030786 4.204 2.627e-05
##
## -------------------------------------------------------------------
## Sigma link function: log
## Migma Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) -3.694 0.03737 -98.85 0.000e+00
## pb(hi_h) 2.506 0.06797 36.87 1.623e-297
##
## -------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) -0.3458 0.615 -0.5624 0.573865
## hi_h 2.6858 0.991 2.7102 0.006725
##
## -------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t-value p-value
## (Intercept) 0.6909 0.1146 6.031 1.633e-09
## hi_h 6.6834 0.9544 7.003 2.505e-12
##
## -------------------------------------------------------------------
## No. of observations in the fit: 3982
## Degrees of Freedom for the fit: 15
## Residual Deg. of Freedom: 3967
## at cycle: 66
##
## Global Deviance: 8880.503
## AIC: 8910.503
## SBC: 9004.846
## *******************************************************************