Definições da sessão
#-----------------------------------------------------------------------
# Prof. Dr. Walmes M. Zeviani
# leg.ufpr.br/~walmes · github.com/walmes
# walmes@ufpr.br · @walmeszeviani
# Laboratory of Statistics and Geoinformation (LEG)
# Department of Statistics · Federal University of Paraná
# 2021-mar-04 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(tidyverse)
library(nlme)
Importação e preparo
#-----------------------------------------------------------------------
# Importação e preparo.
# Importação.
da <- read.table(file = "BD.txt", header = TRUE, sep = "\t")
str(da)
## 'data.frame': 4971 obs. of 2 variables:
## $ gro: int 1 1 1 1 1 1 1 1 1 1 ...
## $ dbh: num 362.9 82.8 89.1 82.8 30 ...
# Preparo.
da <- da %>%
mutate(gro = factor(gro),
dbh = dbh/10) %>%
drop_na()
# %>%
# filter(dbh >= 10)
# Gráfico de densidade.
ggplot(data = da,
mapping = aes(x = dbh)) +
facet_wrap(facets = ~gro) +
# scale_x_log10() +
geom_density()
# Gráfico de distribuição acumulada.
ggplot(data = da,
mapping = aes(x = dbh)) +
facet_wrap(facets = ~gro) +
# scale_x_log10() +
stat_ecdf()
Ajuste da distribuição Weibull 2 parâmetros por máxima verossimilhança
#-----------------------------------------------------------------------
# Ajuste das curvas por máxima verossimilhança.
# Cria o log.
da <- da %>%
mutate(ldbh = log10(dbh))
# Ajuste.
ll_fit <- da %>%
group_by(gro) %>%
do({
fit <- MASS::fitdistr(na.omit(.$dbh), densfun = "weibull")
fit$estimate %>%
as.list() %>%
as.data.frame()
}) %>%
ungroup()
ll_fit
## # A tibble: 4 x 3
## gro shape scale
## <fct> <dbl> <dbl>
## 1 1 0.880 3.74
## 2 2 1.01 4.10
## 3 3 0.983 3.75
## 4 4 1.08 4.54
# Valor máximo para fazer a predição.
max_val <- max(da$dbh, na.rm = TRUE)
# Cria valores para plotar a distribuição ajustada.
fit_curves <- ll_fit %>%
group_by(gro) %>%
do({
x <- seq(0, max_val, length.out = 201)
fy <- dweibull(x, shape = .$shape, scale = .$scale)
Fy <- pweibull(x, shape = .$shape, scale = .$scale)
data.frame(dbh = x, fy = fy, Fy = Fy)
}) %>%
ungroup()
head(fit_curves)
## # A tibble: 6 x 4
## gro dbh fy Fy
## <fct> <dbl> <dbl> <dbl>
## 1 1 0 Inf 0
## 2 1 0.324 0.281 0.110
## 3 1 0.647 0.235 0.192
## 4 1 0.971 0.204 0.263
## 5 1 1.29 0.180 0.325
## 6 1 1.62 0.161 0.380
# Ajuste na escala acumulada.
ggplot(data = da,
mapping = aes(x = dbh)) +
facet_wrap(facets = ~gro) +
# scale_x_log10() +
stat_ecdf() +
geom_line(data = fit_curves,
mapping = aes(x = dbh, y = Fy),
color = "red",
inherit.aes = TRUE)
# Ajuste na escala da densidade.
ggplot(data = da,
mapping = aes(x = dbh)) +
facet_wrap(facets = ~gro) +
geom_density() +
geom_line(data = fit_curves,
mapping = aes(x = dbh, y = fy),
color = "red",
inherit.aes = TRUE)
Obtenção da ECDF para ajuste da Weibull via mínimos quadrados
#-----------------------------------------------------------------------
# Fazer a distribuição de valores de diâmetro DENTRO de cada categoria
# para gerar o par (x, Freq_acum(x)) e com isso ajustar o modelo
# Weibull.
# Número de valores distintos e tamanho de amostra.
da %>%
group_by(gro) %>%
summarise(n_support = length(unique(dbh)),
n = n())
## # A tibble: 4 x 3
## gro n_support n
## <fct> <int> <int>
## 1 1 823 972
## 2 2 1087 1371
## 3 3 1119 1397
## 4 4 854 1028
# Obtém a ECDF.
tb_freq <- da %>%
group_by(gro) %>%
arrange(dbh) %>%
mutate(frel = 1/n()) %>%
group_split() %>%
map(function(tb) {
ecdf_fun <- ecdf(x = c(tb$dbh))
u <- unique(tb$dbh)
data.frame(dbh = u, ecdf = ecdf_fun(u))
}) %>%
bind_rows(.id = "gro") %>%
drop_na()
# Cria o log.
tb_freq <- tb_freq %>%
mutate(ldbh = log10(dbh))
# Confere.
tb_freq %>%
group_by(gro) %>%
summarise_at("ecdf", "max")
## # A tibble: 4 x 2
## gro ecdf
## <chr> <dbl>
## 1 1 1
## 2 2 1
## 3 3 1
## 4 4 1
# NOTE: tudo o que foi feito foi para ter a ECDF!
ggplot(data = tb_freq,
mapping = aes(x = dbh, y = ecdf)) +
facet_wrap(facets = ~gro) +
scale_x_log10() +
geom_point() +
stat_ecdf(data = da,
mapping = aes(x = dbh),
color = "green",
inherit.aes = FALSE)
Ajuste do modelo não linear de Weibull com efeitos aleatórios
library(rpanel)
source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/rp.nls.R")
m0 <- rp.nls(model = ecdf ~ (1 - exp(-((ldbh - a)/b)^c)) ,
data = transform(tb_freq, ldbh = log10(dbh)),
start = list(a = c(-10, 10),
b = c(0, 5),
c = c(0.01, 10)),
xlab = "Diâmetros (cm)",
ylab = "Frequência Relativa Acumulada")
#
#-----------------------------------------------------------------------
# Avaliação dos valores iniciais.
# Valores iniciais.
start <- c(a = -1.9, b = 2.5, c = 6)
# Testa os valores iniciais.
ggplot(data = tb_freq,
mapping = aes(x = ldbh, y = ecdf)) +
facet_wrap(facets = ~gro) +
# scale_x_log10() +
geom_point() +
stat_function(color = "red",
fun = function(x) {
with(as.list(start), {
(1 - exp(-((x - a)/b)^c))
})
})
#-----------------------------------------------------------------------
# Aqui tá sendo feito o ajuste de uma distribuição de probabilidades
# usando as frequências relativas acumuladas em função dos valores
# únicos no suporte e usando mínimos quadrados para ajustar o modelo não
# linear que é a acumulada do modelo de distribuição de
# probabilidades. A abordagem por máxima verossimilhança parece mais
# interessante por várias razões.
# Ajuste do modelo não linear de efeito misto em que todos os
# coeficientes da ECDF de Weibull variam (de efeito aleatório) conforme
# os grupos.
n0 <- nlme(model = ecdf ~ (1 - exp(-((ldbh - a)/b)^c)),
data = tb_freq,
fixed = a + b + c ~ 1,
random = a + b + c ~ 1 | gro,
start = start)
summary(n0)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: ecdf ~ (1 - exp(-((ldbh - a)/b)^c))
## Data: tb_freq
## AIC BIC logLik
## -26891.67 -26829.02 13455.83
##
## Random effects:
## Formula: list(a ~ 1, b ~ 1, c ~ 1)
## Level: gro
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## a 0.082151256 a b
## b 0.077366188 -0.846
## c 0.326555546 0.702 -0.213
## Residual 0.007494148
##
## Fixed effects: a + b + c ~ 1
## Value Std.Error DF t-value p-value
## a -0.7604386 0.04131663 3877 -18.40515 0
## b 1.2827320 0.03894435 3877 32.93756 0
## c 3.0040540 0.16378444 3877 18.34151 0
## Correlation:
## a b
## b -0.847
## c 0.688 -0.203
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2901534 -0.6848663 -0.0772755 0.5894960 3.7256666
##
## Number of Observations: 3883
## Number of Groups: 4
fixef(n0)
## a b c
## -0.7604386 1.2827320 3.0040540
ranef(n0)
## a b c
## 1 -0.131776645 0.09261088 -0.43700012
## 2 -0.004738264 0.01903766 0.07278539
## 3 0.080141505 -0.12200276 -0.10420081
## 4 0.056373404 0.01035422 0.46841554
# Ajuste do modelo não linear de efeito misto em que apenas os
# coeficientes de forma e escala apresentam efeito aleatório, sendo o
# coeficiente de locação (a) é fixo.
n1 <- nlme(model = ecdf ~ (1 - exp(-((ldbh - a)/b)^c)),
data = tb_freq,
fixed = a + b + c ~ 1,
# random = a + b ~ 1 | gro,
random = a + c ~ 1 | gro,
# random = b + c ~ 1 | gro,
start = fixef(n0))
anova(n0, n1)
## Model df AIC BIC logLik Test L.Ratio p-value
## n0 1 10 -26891.67 -26829.02 13455.83
## n1 2 7 -26574.97 -26531.12 13294.49 1 vs 2 322.6936 <.0001
summary(n1)
## Nonlinear mixed-effects model fit by maximum likelihood
## Model: ecdf ~ (1 - exp(-((ldbh - a)/b)^c))
## Data: tb_freq
## AIC BIC logLik
## -26574.97 -26531.12 13294.49
##
## Random effects:
## Formula: list(a ~ 1, c ~ 1)
## Level: gro
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## a 0.044725958 a
## c 0.410741139 0.626
## Residual 0.007813012
##
## Fixed effects: a + b + c ~ 1
## Value Std.Error DF t-value p-value
## a -0.768068 0.02283984 3877 -33.62843 0
## b 1.290104 0.00467111 3877 276.18806 0
## c 3.042440 0.20586229 3877 14.77901 0
## Correlation:
## a b
## b -0.201
## c 0.599 0.062
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.48120946 -0.62887776 -0.05188943 0.54638105 4.03215670
##
## Number of Observations: 3883
## Number of Groups: 4
fixef(n1)
## a b c
## -0.768068 1.290104 3.042440
ranef(n1)
## a c
## 1 -0.04128617 -0.663469834
## 2 0.01437813 0.002846493
## 3 -0.04006093 0.240018002
## 4 0.06696896 0.420605339
coef(n1)
## a b c
## 1 -0.8093542 1.290104 2.378970
## 2 -0.7536899 1.290104 3.045287
## 3 -0.8081290 1.290104 3.282458
## 4 -0.7010991 1.290104 3.463046
#-----------------------------------------------------------------------
# Fazer a predição das curvas.
extr_val <- extendrange(tb_freq$ldbh)
grid <- crossing(gro = unique(tb_freq$gro),
ldbh = seq(extr_val[1],
extr_val[2],
length.out = 201))
grid$y <- predict(n1, newdata = grid)
# Resultado do ajuste.
ggplot(data = tb_freq,
mapping = aes(x = ldbh, y = ecdf)) +
facet_wrap(facets = ~gro) +
geom_point() +
geom_line(data = grid,
mapping = aes(x = ldbh, y = y),
inherit.aes = FALSE,
color = "magenta")
# ATTENTION: problemas com a abordagem de mínimos quadrados sobre a
# ECDF.
#
# O tamanho de amostra não é o número de árvores medidas mas sim o
# número de valores únicos de diâmetro. Com isso, dependendo da escala
# do aparelho de medida, pode haver uma diferença bem grande entre as
# quantidades tamanho de amostra e número de valores no suporte. Por
# isso a abordagem de verossimilhança é mais indicada.
#