library(gamlss)

## Loading required package: splines

## Loading required package: gamlss.data

## 
## Attaching package: 'gamlss.data'

## The following object is masked from 'package:datasets':
## 
##     sleep

## Loading required package: gamlss.dist

## Loading required package: nlme

## Loading required package: parallel

##  **********   GAMLSS Version 5.5-0  **********

## For more on GAMLSS look at https://www.gamlss.com/

## Type gamlssNews() to see new features/changes/bug fixes.

library(readxl)

dados <- read_excel("dados_gamlss.xlsx")
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     3

d=dap; h= altura da arvore; hi=altura relativa da medida do diametro;

#di=diametro

estimação de modelos segmentados

$ y= d/di; 𝑥 = ℎi/ℎ$

attach(dados)

MODELO SIMPLES DE DUAS VARIÁVEIS: Y_LM E X

tabela <- data.frame(hi/h, di)

#gráfico da relação

plot(tabela,
     xlab = "Altura relativa (hi/h)",
     ylab = "Diâmetro relativo (di/d)",
     main = "Perfil de afilamento")

#ajusatndo modelos com GAMLSS

polinomio 3º grau

library(gamlss)

# transformar hi/h -> x
x <- hi/h

distribuição normal

m.p3_1 <- gamlss(di ~ x + I(x^2) + I(x^3),
               family = NO)

## GAMLSS-RS iteration 1: Global Deviance = 20562.24 
## GAMLSS-RS iteration 2: Global Deviance = 20562.24

os residuos não indicaram bom ajuste

escolher nova distribuição para os dados

Melhor distribuição

dist_m.p3_1 <- chooseDist(m.p3_1, type="realAll") ; dist_m.p3_1

## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in RS(): Algorithm RS has not yet converged

## minimum GAIC(k= 2 ) family: BCPEo 
## minimum GAIC(k= 3.84 ) family: BCPEo 
## minimum GAIC(k= 8.29 ) family: BCPEo

##                 2     3.84     8.29
## NO       20572.24 20581.44 20603.69
## GU       20920.92 20930.12 20952.37
## RG       21085.35 21094.55 21116.80
## LO       20554.54 20563.74 20585.99
## NET      20754.21 20763.41 20785.66
## TF       20544.13 20555.17 20581.87
## TF2      20544.13 20555.17 20581.87
## PE       20542.78 20553.82 20580.52
## PE2      20542.78 20553.82 20580.52
## SN1      20574.24 20585.28 20611.98
## SN2      20572.05 20583.09 20609.79
## exGAUS   20574.16 20585.20 20611.90
## SHASH    20551.14 20564.02 20595.17
## SHASHo   20553.15 20566.03 20597.18
## SHASHo2  20544.95 20557.83 20588.98
## EGB2     20555.59 20568.47 20599.62
## JSU      20544.37 20557.25 20588.40
## JSUo     20604.71 20617.59 20648.74
## SEP1     20545.06 20557.94 20589.09
## SEP2     20544.97 20557.85 20589.00
## SEP3     20544.67 20557.55 20588.70
## SEP4     20544.14 20557.02 20588.17
## ST1      20546.83 20559.71 20590.86
## ST2      20546.46 20559.34 20590.49
## ST3      20545.37 20558.25 20589.40
## ST4      20544.58 20557.46 20588.61
## ST5      20548.54 20561.42 20592.57
## SST      20545.31 20558.19 20589.34
## GT       20544.55 20557.43 20588.58
## EXP      26927.76 26935.12 26952.92
## GA       19518.72 19527.92 19550.17
## IG       19930.63 19939.83 19962.08
## LOGNO    19746.32 19755.52 19777.77
## LOGNO2   19746.32 19755.52 19777.77
## WEI      19278.67 19287.87 19310.12
## WEI2     20037.26 20046.46 20068.71
## WEI3     19278.67 19287.87 19310.12
## IGAMMA   20074.55 20083.75 20106.00
## PARETO2  27297.35 27306.55 27328.80
## PARETO2o 27297.46 27306.66 27328.91
## GP       27297.35 27306.55 27328.80
## BCCG     19416.64 19427.68 19454.38
## BCCGo    19323.90 19334.94 19361.64
## GG       19279.71 19290.75 19317.45
## GIG      19520.72 19531.76 19558.46
## LNO      19746.32 19755.52 19777.77
## BCTo     19325.90 19338.78 19369.93
## BCT      19418.64 19431.52 19462.67
## BCPEo    19233.55 19246.43 19277.58
## BCPE     19333.53 19346.41 19377.56
## GB2      19396.65 19409.53 19440.68

best_dist_m.p3_1 <- dist_m.p3_1[c("EXP","GA","IG","LOGNO","WEI","WEI2","BCCG","GIG","BCT","LO","NO","exGAUS","PE","TF2"),]
min(best_dist_m.p3_1[,1]) # p/ k=2.00  #melhor WEI

## [1] 19278.67

min(best_dist_m.p3_1[,2]) # p/ k=3.84  #melhor WEI

## [1] 19287.87

min(best_dist_m.p3_1[,3]) # p/ k=6.74  #melhor WEI

## [1] 19310.12

Ajustando com distribuição selecionada bcpe0

m.pm_2 <- gamlss(di ~ x + I(x^2) + I(x^3),
                 family = BCPEo()
            )

## GAMLSS-RS iteration 1: Global Deviance = 19348.37 
## GAMLSS-RS iteration 2: Global Deviance = 19225.51 
## GAMLSS-RS iteration 3: Global Deviance = 19221.09 
## GAMLSS-RS iteration 4: Global Deviance = 19219.95 
## GAMLSS-RS iteration 5: Global Deviance = 19219.65 
## GAMLSS-RS iteration 6: Global Deviance = 19219.57 
## GAMLSS-RS iteration 7: Global Deviance = 19219.56 
## GAMLSS-RS iteration 8: Global Deviance = 19219.55 
## GAMLSS-RS iteration 9: Global Deviance = 19219.55 
## GAMLSS-RS iteration 10: Global Deviance = 19219.55

gráfico e resumo

summary(m.pm_2)

## ******************************************************************
## Family:  c("BCPEo", "Box-Cox Power Exponential-orig.") 
## 
## Call:  gamlss(formula = di ~ x + I(x^2) + I(x^3), family = BCPEo()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.993487   0.009575  312.64   <2e-16 ***
## x           -2.800681   0.109510  -25.57   <2e-16 ***
## I(x^2)       4.701200   0.320250   14.68   <2e-16 ***
## I(x^3)      -4.333045   0.255315  -16.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.3856     0.0114  -121.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.15523    0.04682   24.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.07998    0.04076   26.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3982 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  3975 
##                       at cycle:  10 
##  
## Global Deviance:     19219.55 
##             AIC:     19233.55 
##             SBC:     19277.58 
## ******************************************************************

plot(m.pm_2)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.0008648483 
##                        variance   =  1.000517 
##                coef. of skewness  =  -0.0007431948 
##                coef. of kurtosis  =  2.978641 
## Filliben correlation coefficient  =  0.9995483 
## ******************************************************************

wp(m.pm_2)

## Warning in wp(m.pm_2): Some points are missed out 
## increase the y limits using ylim.all

AIC

AIC(m.p3_1, m.pm_2) # modelo pm_2= polinomial com wei -aic -12730,490

##        df      AIC
## m.pm_2  7 19233.55
## m.p3_1  5 20572.24

modelo de Prodan / polinomio de 4? grau

mod_pol_4 <- gamlss(di ~ x + I(x^2) + I(x^3) + I(x^4),
              family = NO)

## GAMLSS-RS iteration 1: Global Deviance = 20467.65 
## GAMLSS-RS iteration 2: Global Deviance = 20467.65

summary(mod_pol_4)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = di ~ x + I(x^2) + I(x^3) + I(x^4),      family = NO) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.5818     0.1679 122.562   <2e-16 ***
## x            -71.2103     3.1773 -22.412   <2e-16 ***
## I(x^2)       222.3415    16.1219  13.791   <2e-16 ***
## I(x^3)      -333.4121    29.2956 -11.381   <2e-16 ***
## I(x^4)       169.4391    17.3181   9.784   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.15108    0.01121   102.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3982 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  3976 
##                       at cycle:  2 
##  
## Global Deviance:     20467.65 
##             AIC:     20479.65 
##             SBC:     20517.39 
## ******************************************************************

plot(mod_pol_4) # resíduos assimétricos. É possível obter um ajuste melhor?

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -7.107643e-17 
##                        variance   =  1.000251 
##                coef. of skewness  =  -0.09472375 
##                coef. of kurtosis  =  3.513196 
## Filliben correlation coefficient  =  0.9983102 
## ******************************************************************

wp(mod_pol_4) # Melhorar a modelagem!

## Warning in wp(mod_pol_4): Some points are missed out 
## increase the y limits using ylim.all

escolher a distribuição gamlss

dist_mod_pol_4 <- chooseDist(mod_pol_4, type="realAll") ; dist_mod_pol_4

## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in RS(): Algorithm RS has not yet converged

## minimum GAIC(k= 2 ) family: BCPEo 
## minimum GAIC(k= 3.84 ) family: BCPEo 
## minimum GAIC(k= 8.29 ) family: BCPEo

##                 2     3.84     8.29
## NO       20479.65 20490.69 20517.39
## GU       20768.92 20779.96 20806.66
## RG       21034.87 21045.91 21072.61
## LO       20452.51 20463.55 20490.25
## NET      20653.94 20664.98 20691.68
## TF       20446.39 20459.27 20490.42
## TF2      20446.39 20459.27 20490.42
## PE       20442.78 20455.66 20486.81
## PE2      20442.78 20455.66 20486.81
## SN1      20481.65 20494.53 20525.68
## SN2      20476.70 20489.58 20520.73
## exGAUS   20484.07 20496.95 20528.10
## SHASH    20447.57 20462.29 20497.89
## SHASHo   20450.44 20465.16 20500.76
## SHASHo2  20443.01 20457.73 20493.33
## EGB2     20453.30 20468.02 20503.62
## JSU      20443.77 20458.49 20494.09
## JSUo     20502.69 20517.41 20553.01
## SEP1     20444.94 20459.66 20495.26
## SEP2     20444.63 20459.35 20494.95
## SEP3     20443.50 20458.22 20493.82
## SEP4     20440.84 20455.56 20491.16
## ST1      20448.23 20462.95 20498.55
## ST2      20448.69 20463.41 20499.01
## ST3      20446.56 20461.28 20496.88
## ST4      20443.04 20457.76 20493.36
## ST5      20451.01 20465.73 20501.33
## SST      20446.38 20461.10 20496.70
## GT       20445.44 20460.16 20495.76
## EXP      26926.29 26935.49 26957.74
## GA       19471.16 19482.20 19508.90
## IG       19879.45 19890.49 19917.19
## LOGNO    19695.73 19706.77 19733.47
## LOGNO2   19695.73 19706.77 19733.47
## WEI      19242.86 19253.90 19280.60
## WEI2     20011.90 20022.94 20049.64
## WEI3     19242.86 19253.90 19280.60
## IGAMMA   20022.20 20033.24 20059.94
## PARETO2  27295.91 27306.95 27333.65
## PARETO2o 27297.53 27308.57 27335.27
## GP       27295.91 27306.95 27333.65
## BCCG     19288.51 19301.39 19332.54
## BCCGo    19281.34 19294.22 19325.37
## GG       19242.68 19255.56 19286.71
## GIG      19473.17 19486.05 19517.20
## LNO      19695.73 19706.77 19733.47
## BCTo     19283.34 19298.06 19333.66
## BCT      19290.51 19305.23 19340.83
## BCPEo    19203.40 19218.12 19253.72
## BCPE     19213.47 19228.19 19263.79
## GB2      19355.33 19370.05 19405.65

best_dist_mod_pol_4 <- dist_mod_pol_4[c("EXP","GA","IG","LOGNO","WEI","WEI2","BCCG","GIG","BCT","LO","NO","exGAUS","PE","TF2"),] ; best_dist_mod_pol_4

##               2     3.84     8.29
## EXP    26926.29 26935.49 26957.74
## GA     19471.16 19482.20 19508.90
## IG     19879.45 19890.49 19917.19
## LOGNO  19695.73 19706.77 19733.47
## WEI    19242.86 19253.90 19280.60
## WEI2   20011.90 20022.94 20049.64
## BCCG   19288.51 19301.39 19332.54
## GIG    19473.17 19486.05 19517.20
## BCT    19290.51 19305.23 19340.83
## LO     20452.51 20463.55 20490.25
## NO     20479.65 20490.69 20517.39
## exGAUS 20484.07 20496.95 20528.10
## PE     20442.78 20455.66 20486.81
## TF2    20446.39 20459.27 20490.42

best_dist_mod_pol_4 <-  na.omit(best_dist_mod_pol_4) ; best_dist_mod_pol_4

##               2     3.84     8.29
## EXP    26926.29 26935.49 26957.74
## GA     19471.16 19482.20 19508.90
## IG     19879.45 19890.49 19917.19
## LOGNO  19695.73 19706.77 19733.47
## WEI    19242.86 19253.90 19280.60
## WEI2   20011.90 20022.94 20049.64
## BCCG   19288.51 19301.39 19332.54
## GIG    19473.17 19486.05 19517.20
## BCT    19290.51 19305.23 19340.83
## LO     20452.51 20463.55 20490.25
## NO     20479.65 20490.69 20517.39
## exGAUS 20484.07 20496.95 20528.10
## PE     20442.78 20455.66 20486.81
## TF2    20446.39 20459.27 20490.42

min(best_dist_mod_pol_4[,1]) # p/ k=2.00  #melhor dist bcpeo

## [1] 19242.86

min(best_dist_mod_pol_4[,2]) # p/ k=3.84  #melhor dist bcpeo

## [1] 19253.9

min(best_dist_mod_pol_4[,3]) # p/ k=6.74  #melhor dist bcpeo

## [1] 19280.6

mod_pol_4_1 <- gamlss(di ~ x + I(x^2) + I(x^3) + I(x^4),
                      family = BCPEo ()) # familia BCPEo

## GAMLSS-RS iteration 1: Global Deviance = 19318.2 
## GAMLSS-RS iteration 2: Global Deviance = 19192.52 
## GAMLSS-RS iteration 3: Global Deviance = 19188.7 
## GAMLSS-RS iteration 4: Global Deviance = 19187.72 
## GAMLSS-RS iteration 5: Global Deviance = 19187.47 
## GAMLSS-RS iteration 6: Global Deviance = 19187.42 
## GAMLSS-RS iteration 7: Global Deviance = 19187.4 
## GAMLSS-RS iteration 8: Global Deviance = 19187.4 
## GAMLSS-RS iteration 9: Global Deviance = 19187.4

summary(mod_pol_4_1) #

## ******************************************************************
## Family:  c("BCPEo", "Box-Cox Power Exponential-orig.") 
## 
## Call:  gamlss(formula = di ~ x + I(x^2) + I(x^3) + I(x^4),  
##     family = BCPEo()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.02670    0.01156 261.938  < 2e-16 ***
## x            -3.82659    0.21724 -17.614  < 2e-16 ***
## I(x^2)       10.59711    1.11757   9.482  < 2e-16 ***
## I(x^3)      -15.46808    2.04727  -7.555 5.15e-14 ***
## I(x^4)        6.62161    1.21836   5.435 5.81e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.39008    0.01142  -121.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.14425    0.04747   24.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.04808    0.04016    26.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3982 
## Degrees of Freedom for the fit:  8
##       Residual Deg. of Freedom:  3974 
##                       at cycle:  9 
##  
## Global Deviance:     19187.4 
##             AIC:     19203.4 
##             SBC:     19253.72 
## ******************************************************************

plot(mod_pol_4_1)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.001205362 
##                        variance   =  1.001037 
##                coef. of skewness  =  0.001121911 
##                coef. of kurtosis  =  3.000216 
## Filliben correlation coefficient  =  0.9994539 
## ******************************************************************

wp(mod_pol_4_1)

## Warning in wp(mod_pol_4_1): Some points are missed out 
## increase the y limits using ylim.all

# wei

mod_pol4_wei <- gamlss(di ~ x + I(x^2) + I(x^3) + I(x^4),
                      family = WEI ()) # familia wei

## GAMLSS-RS iteration 1: Global Deviance = 19384.17 
## GAMLSS-RS iteration 2: Global Deviance = 19233.78 
## GAMLSS-RS iteration 3: Global Deviance = 19230.89 
## GAMLSS-RS iteration 4: Global Deviance = 19230.86 
## GAMLSS-RS iteration 5: Global Deviance = 19230.86

aic e resíduos distribuição wei

summary(mod_pol4_wei) #

## ******************************************************************
## Family:  c("WEI", "Weibull") 
## 
## Call:  gamlss(formula = di ~ x + I(x^2) + I(x^3) + I(x^4),  
##     family = WEI()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.11633    0.01171 266.155  < 2e-16 ***
## x            -3.94204    0.22126 -17.816  < 2e-16 ***
## I(x^2)       11.27626    1.12291  10.042  < 2e-16 ***
## I(x^3)      -16.64565    2.03170  -8.193 3.40e-16 ***
## I(x^4)        7.23549    1.19439   6.058 1.51e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.5031     0.0124   121.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3982 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  3976 
##                       at cycle:  5 
##  
## Global Deviance:     19230.86 
##             AIC:     19242.86 
##             SBC:     19280.59 
## ******************************************************************

plot(mod_pol4_wei)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.0008760225 
##                        variance   =  0.9907635 
##                coef. of skewness  =  0.07994959 
##                coef. of kurtosis  =  2.645348 
## Filliben correlation coefficient  =  0.9983432 
## ******************************************************************

wp(mod_pol4_wei)

## Warning in wp(mod_pol4_wei): Some points are missed out 
## increase the y limits using ylim.all

polinomio de 5? grau

mod_pol_5 <- gamlss(di ~ x + I(x^2) + I(x^3) + I(x^4)+ I(x^5),
                    family = NO)

## GAMLSS-RS iteration 1: Global Deviance = 20459.3 
## GAMLSS-RS iteration 2: Global Deviance = 20459.3

summary(mod_pol_5)

## ******************************************************************
## Family:  c("NO", "Normal") 
## 
## Call:  gamlss(formula = di ~ x + I(x^2) + I(x^3) + I(x^4) +  
##     I(x^5), family = NO) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  identity
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.8483     0.1925 108.293  < 2e-16 ***
## x            -82.9567     5.2301 -15.861  < 2e-16 ***
## I(x^2)       328.6212    40.8767   8.039 1.18e-15 ***
## I(x^3)      -680.3780   126.0024  -5.400 7.06e-08 ***
## I(x^4)       635.5655   165.4542   3.841 0.000124 ***
## I(x^5)      -219.7360    77.5317  -2.834 0.004618 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.15003    0.01121   102.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3982 
## Degrees of Freedom for the fit:  7
##       Residual Deg. of Freedom:  3975 
##                       at cycle:  2 
##  
## Global Deviance:     20459.3 
##             AIC:     20473.3 
##             SBC:     20517.33 
## ******************************************************************

plot(mod_pol_5) # res????duos assimétricos. É poss????vel obter um ajuste melhor?

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -3.891138e-16 
##                        variance   =  1.000251 
##                coef. of skewness  =  -0.1039208 
##                coef. of kurtosis  =  3.518281 
## Filliben correlation coefficient  =  0.9982397 
## ******************************************************************

wp(mod_pol_5) # Melhorar a modelagem!

## Warning in wp(mod_pol_5): Some points are missed out 
## increase the y limits using ylim.all

escolhendo a distribuição com gamlss

dist_mod_pol_5 <- chooseDist(mod_pol_5, type="realAll") ; dist_mod_pol_5

## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged
## Warning in RS(): Algorithm RS has not yet converged

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = fv)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in `formals<-`(`*tmp*`, envir = new.env(), value = alist(nu = nu)):
## 'fun' is not a function

## Warning in body(fun): argument is not a function

## Warning in RS(): Algorithm RS has not yet converged

## minimum GAIC(k= 2 ) family: BCPEo 
## minimum GAIC(k= 3.84 ) family: BCPEo 
## minimum GAIC(k= 8.29 ) family: BCPEo

##                 2     3.84     8.29
## NO       20473.30 20486.18 20517.33
## GU       20745.74 20758.62 20789.77
## RG       21036.03 21048.91 21080.06
## LO       20444.52 20457.40 20488.55
## NET      20644.10 20656.98 20688.13
## TF       20439.08 20453.80 20489.40
## TF2      20439.08 20453.80 20489.40
## PE       20435.92 20450.64 20486.24
## PE2      20435.92 20450.64 20486.24
## SN1      20475.31 20490.03 20525.63
## SN2      20469.26 20483.98 20519.58
## exGAUS   20478.31 20493.03 20528.63
## SHASH    20439.53 20456.09 20496.14
## SHASHo   20442.52 20459.08 20499.13
## SHASHo2  20435.35 20451.91 20491.96
## EGB2     20444.96 20461.52 20501.57
## JSU      20435.55 20452.11 20492.16
## JSUo     20493.98 20510.54 20550.59
## SEP1     20438.30 20454.86 20494.91
## SEP2     20438.09 20454.65 20494.70
## SEP3     20435.62 20452.18 20492.23
## SEP4     20433.33 20449.89 20489.94
## ST1      20440.66 20457.22 20497.27
## ST2      20441.36 20457.92 20497.97
## ST3      20438.59 20455.15 20495.20
## ST4      20434.46 20451.02 20491.07
## ST5      20443.43 20459.99 20500.04
## SST      20438.34 20454.90 20494.95
## GT       20438.30 20454.86 20494.91
## EXP      26928.29 26939.33 26966.03
## GA       19473.16 19486.04 19517.19
## IG       19881.31 19894.19 19925.34
## LOGNO    19697.73 19710.61 19741.76
## LOGNO2   19697.73 19710.61 19741.76
## WEI      19244.58 19257.46 19288.61
## WEI2     20013.81 20026.69 20057.84
## WEI3     19244.58 19257.46 19288.61
## IGAMMA   20024.15 20037.03 20068.18
## PARETO2  27297.90 27310.78 27341.93
## PARETO2o 27298.82 27311.70 27342.85
## GP       27297.90 27310.78 27341.93
## BCCG     19284.32 19299.04 19334.64
## BCCGo    19283.27 19297.99 19333.59
## GG       19244.46 19259.18 19294.78
## GIG      19475.16 19489.88 19525.48
## LNO      19697.73 19710.61 19741.76
## BCTo     19285.27 19301.83 19341.88
## BCT      19286.32 19302.88 19342.93
## BCPEo    19204.45 19221.01 19261.06
## BCPE     19207.58 19224.14 19264.19
## GB2      19357.32 19373.88 19413.93

best_dist_mod_pol_5 <- dist_mod_pol_5[c("EXP","GA","IG","LOGNO","WEI","WEI2","BCCG","GIG","BCT","LO","NO","exGAUS","PE","TF2"),] ; best_dist_mod_pol_5

##               2     3.84     8.29
## EXP    26928.29 26939.33 26966.03
## GA     19473.16 19486.04 19517.19
## IG     19881.31 19894.19 19925.34
## LOGNO  19697.73 19710.61 19741.76
## WEI    19244.58 19257.46 19288.61
## WEI2   20013.81 20026.69 20057.84
## BCCG   19284.32 19299.04 19334.64
## GIG    19475.16 19489.88 19525.48
## BCT    19286.32 19302.88 19342.93
## LO     20444.52 20457.40 20488.55
## NO     20473.30 20486.18 20517.33
## exGAUS 20478.31 20493.03 20528.63
## PE     20435.92 20450.64 20486.24
## TF2    20439.08 20453.80 20489.40

min(best_dist_mod_pol_5[,1]) # p/ k=2.00  #melhor dist ? Bcpeo

## [1] 19244.58

min(best_dist_mod_pol_5[,2]) # p/ k=3.84  #melhor dist ? bcpeo

## [1] 19257.46

min(best_dist_mod_pol_5[,3]) # p/ k=6.74  #melhor dist ? bcpeo

## [1] 19288.61

mod_pol_5_1 <- gamlss(di ~ x + I(x^2) + I(x^3) + I(x^4)+ I(x^5),
                      family = BCPEo())

## GAMLSS-RS iteration 1: Global Deviance = 19317.83 
## GAMLSS-RS iteration 2: Global Deviance = 19191.51 
## GAMLSS-RS iteration 3: Global Deviance = 19187.75 
## GAMLSS-RS iteration 4: Global Deviance = 19186.77 
## GAMLSS-RS iteration 5: Global Deviance = 19186.53 
## GAMLSS-RS iteration 6: Global Deviance = 19186.47 
## GAMLSS-RS iteration 7: Global Deviance = 19186.46 
## GAMLSS-RS iteration 8: Global Deviance = 19186.45 
## GAMLSS-RS iteration 9: Global Deviance = 19186.45

summary(mod_pol_5_1) #

## ******************************************************************
## Family:  c("BCPEo", "Box-Cox Power Exponential-orig.") 
## 
## Call:  gamlss(formula = di ~ x + I(x^2) + I(x^3) + I(x^4) +  
##     I(x^5), family = BCPEo()) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.03256    0.01309 231.619  < 2e-16 ***
## x            -4.08726    0.34567 -11.824  < 2e-16 ***
## I(x^2)       12.95239    2.66667   4.857 1.24e-06 ***
## I(x^3)      -23.15364    8.14727  -2.842  0.00451 ** 
## I(x^4)       16.94824   10.64520   1.592  0.11144    
## I(x^5)       -4.86938    4.97542  -0.979  0.32779    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.39010    0.01142  -121.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  identity 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.14554    0.04742   24.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Tau link function:  log 
## Tau Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.05125    0.04035   26.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  3982 
## Degrees of Freedom for the fit:  9
##       Residual Deg. of Freedom:  3973 
##                       at cycle:  9 
##  
## Global Deviance:     19186.45 
##             AIC:     19204.45 
##             SBC:     19261.06 
## ******************************************************************

plot(mod_pol_5_1)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.001089763 
##                        variance   =  1.001057 
##                coef. of skewness  =  0.0009420844 
##                coef. of kurtosis  =  2.999261 
## Filliben correlation coefficient  =  0.9994755 
## ******************************************************************

wp(mod_pol_5_1)

## Warning in wp(mod_pol_5_1): Some points are missed out 
## increase the y limits using ylim.all

Ajuste de modelos de afilamento com gamlss

eliel

2026-02-26

d=dap; h= altura da arvore; hi=altura relativa da medida do diametro;

estimação de modelos segmentados

$ y= d/di; 𝑥 = ℎi/ℎ$

MODELO SIMPLES DE DUAS VARIÁVEIS: Y_LM E X

polinomio 3º grau

distribuição normal

os residuos não indicaram bom ajuste

escolher nova distribuição para os dados

Melhor distribuição

Ajustando com distribuição selecionada bcpe0

gráfico e resumo

AIC

modelo de Prodan / polinomio de 4? grau

escolher a distribuição gamlss

aic e resíduos distribuição wei

polinomio de 5? grau

escolhendo a distribuição com gamlss