Dados

Foram utilizados dados oriundos de inventário florestal em povoamentos de Pinus taeda entre 4 e 26 anos de idade na região de Campo do Tenente, Paraná.

library(readxl)
library(ggplot2)

# Parâmetros
nbs <- 1000 #qtd de simulações
id_start <- 0 #idade inicial
id_end <- 20 #idade final

# Importar dados
dados <- read_excel("dados_processados.xlsx")

# Visualizar dados
plot(HDOM~IDADE,dados,xlab='Idade (anos)',ylab='Altura Dominante (m)',xlim=c(0,30))

Ajuste do modelo de Chapman-Richards

# Ajustar modelo
chapman_richards <- nls(HDOM ~ b0*(1-exp(-b1*IDADE))^b2,
                        start=list(b0=30,b1=0.1,b2=1),
                        data = dados)

# Coeficientes
b0_x <- coef(chapman_richards)[1]
b1_x <- coef(chapman_richards)[2]
b2_x <- coef(chapman_richards)[3]

# Erro padrão dos coeficientes
b0_se <- sqrt(diag(vcov(chapman_richards)))[1]
b1_se <- sqrt(diag(vcov(chapman_richards)))[2]
b2_se <- sqrt(diag(vcov(chapman_richards)))[3]

# Amostragem de coeficientes
b0_dist <- rnorm(nbs,b0_x,b0_se)
b1_dist <- rnorm(nbs,b1_x,b1_se)
b2_dist <- rnorm(nbs,b2_x,b2_se)

Distribuição probabilística dos coeficientes do modelo ajustado

# Plotar distribuições
par(mfrow=c(2,2))
hist(b0_dist, main = 'Distribuição b0');abline(v=b0_x, col='red', lwd=3, lty=2)
hist(b1_dist, main = 'Distribuição b1');abline(v=b1_x, col='red', lwd=3, lty=2)
hist(b2_dist, main = 'Distribuição b2');abline(v=b2_x, col='red', lwd=3, lty=2)
boxplot(data.frame(b0=1-(b0_dist/mean(b0_dist)),
                   b1=1-(b1_dist/mean(b1_dist)),
                   b2=1-(b2_dist/mean(b2_dist))),
        outline = F, ylab='Variação % em torno da média')

Simulação de Monte Carlo

A simulação de Monte Carlo foi realizada considerando a predição da altura dominante e a projeção a partir de dados conhecidos pelo método da diferença algébrica

# Criar data.frame com os coeficientes amostrados
coef_mc <- data.frame(b0 = b0_dist,
                      b1 = b1_dist,
                      b2= b2_dist)

# Transformar em lista
coef_mc_list <- split(coef_mc,seq(nrow(coef_mc)))

# Função de predição
cr_pred <- function(x){
  predicts <- rnorm(nrow(dados),x$b0*(1-exp(-x$b1*id))^x$b2,sigma(chapman_richards)/mean(dados$HDOM)*x$b0*(1-exp(-x$b1*id))^x$b2)
  return(predicts)
}

# Função de projeção
cr_proj <- function(x){
  projections <- dados$HDOM*((1-exp(-x$b1*id))/(1-exp(-x$b1*dados$IDADE)))^x$b2
  return(projections)
}

# Objetos necessários para as simulações
n <- 0
mc_mean <- vector()
model_var <- vector()
design_var <- vector()
total_var <- vector()
ic_inf <- vector()
ic_sup <- vector()

# Inicia as simulações das predições
for(id in id_start:id_end){
  n <- n+1
  mc_proj15 <- lapply(coef_mc_list, cr_pred)
  mc_proj15_means <- unlist(lapply(mc_proj15, mean))
  mc_mean[n] <- mean(mc_proj15_means)
  model_var[n] <- sum((mc_proj15_means-mc_mean[n])^2)/nbs
  design_var[n] <- sum(unlist(lapply(mc_proj15, var)))/nbs
  total_var[n] <- model_var[n]+design_var[n]
  ic_inf[n] <- mc_mean[n]-1.96*sqrt(total_var[n])
  ic_sup[n] <- mc_mean[n]+1.96*sqrt(total_var[n])
}

# Gera a tabela final
tabela_pred <- data.frame(idade=id_start:id_end,
                           model_var=model_var,
                           design_var=design_var,
                           total_var=total_var,
                           ic_inf=ic_inf,
                           ic_sup=ic_sup,
                           predicao=mc_mean)

# Gera gráficos das predições
pred_uncertainty <- ggplot()+
  geom_line(aes(y=total_var,x=idade, linetype = 'Variância total'),tabela_pred, size=1)+
  geom_line(aes(y=design_var,x=idade, linetype = 'Variância do desenho'),tabela_pred, size=1)+
  geom_line(aes(y=model_var,x=idade, linetype = 'Variância do modelo'),tabela_pred, size=1)+
  xlab('Idade (anos)')+
  ylab('Variância')+
  ggtitle('Incerteza associada ao procedimento de predição')+
  scale_x_continuous(expand=c(0,0))+
  scale_y_continuous(expand=c(0,0),limits=c(0,6))+
  scale_linetype_manual(name=NULL,values=c(3,2,1))+
  theme_bw()

pred_curva <- ggplot()+
  geom_ribbon(aes(x=idade,ymin=ic_inf,ymax=ic_sup, fill = 'Intervalo de confiança'), linetype = 2, color = 'black', tabela_pred)+
  geom_line(aes(y=predicao,x=idade, linetype = 'Curva de predição'),tabela_pred, size=1)+
  xlab('Idade (anos)')+
  ylab('Variância')+
  scale_x_continuous(expand=c(0,0))+
  scale_y_continuous(expand=c(0,0))+
  scale_linetype_manual(name=NULL,values=c(1))+
  scale_fill_manual(name=NULL,values='gray')+
  theme_bw()


# Objetos necessários para as simulações
n <- 0
mc_mean <- vector()
model_var <- vector()
design_var <- vector()
total_var <- vector()
ic_inf <- vector()
ic_sup <- vector()

# Inicia as simulações das predições
for(id in id_start:id_end){
  n <- n+1
  mc_proj15 <- lapply(coef_mc_list, cr_proj)
  mc_proj15_means <- unlist(lapply(mc_proj15, mean))
  mc_mean[n] <- mean(mc_proj15_means)
  model_var[n] <- sum((mc_proj15_means-mc_mean[n])^2)/nbs
  design_var[n] <- sum(unlist(lapply(mc_proj15, var)))/nbs
  total_var[n] <- model_var[n]+design_var[n]
  ic_inf[n] <- mc_mean[n]-1.96*sqrt(total_var[n])
  ic_sup[n] <- mc_mean[n]+1.96*sqrt(total_var[n])
}

# Gera a tabela final
tabela_proj <- data.frame(idade=id_start:id_end,
                           model_var=model_var,
                           design_var=design_var,
                           total_var=total_var,
                           ic_inf=ic_inf,
                           ic_sup=ic_sup,
                           projecao=mc_mean)

# Gera gráficos das projeções
library(ggplot2)
proj_uncertainty <- ggplot()+
  geom_line(aes(y=total_var,x=idade, linetype = 'Variância total'),tabela_proj, size=1)+
  geom_line(aes(y=design_var,x=idade, linetype = 'Variância do desenho'),tabela_proj, size=1)+
  geom_line(aes(y=model_var,x=idade, linetype = 'Variância do modelo'),tabela_proj, size=1)+
  xlab('Idade (anos)')+
  ylab('Variância')+
  ggtitle('Incerteza associada ao procedimento de projeção')+
  scale_x_continuous(expand=c(0,0))+
  scale_y_continuous(expand=c(0,0),limits=c(0,6))+
  scale_linetype_manual(name=NULL,values=c(3,2,1))+
  theme_bw()

proj_curva <- ggplot()+
  geom_ribbon(aes(x=idade,ymin=ic_inf,ymax=ic_sup, fill = 'Intervalo de confiança'), linetype = 2, color = 'black', tabela_proj)+
  geom_line(aes(y=projecao,x=idade, linetype = 'Curva de projeção'),tabela_proj, size=1)+
  xlab('Idade (anos)')+
  ylab('Variância')+
  scale_x_continuous(expand=c(0,0))+
  scale_y_continuous(expand=c(0,0))+
  scale_linetype_manual(name=NULL,values=c(1))+
  scale_fill_manual(name=NULL,values='gray')+
  theme_bw()

# Plotar todos os gráficos
gridExtra::grid.arrange(pred_uncertainty,pred_curva,proj_uncertainty,proj_curva)

Composição da variância aos 15 anos

# Componentes da variância total aos 15 anos
tabela <- round(tabela_pred[16,c(2,3,4)],3)
tabela[2,] <- round(tabela_proj[16,c(2,3,4)],3)
colnames(tabela) <- c('Variância do modelo','Variância do desenho', 'Variância total')
rownames(tabela) <- c('Predição','Projeção')
DT::datatable(tabela,options = list(dom='t'))