Tabela de conteúdos
Objetivo
O objetivo desse curso é apresentar os procedimentos necessários para se construir Modelos Dendrométricos. Modelos dendrométricos são modelos empíricos utilizados para obter predições de atributos de árvores individuais que são de mensuração difícil (altura) ou de mensuração destrutiva (volume, biomassa).
Ao final desse curso você deverá ser capaz de:
- ajustar modelos dendrométricos utilizando regressão linear;
- verificar a qualidade das predições dos modelos ajustados e selecionar o modelo mais adequado;
- gerar predições a partir do modelo selecionado para dados de árvores individuais;
- agregar dados de árvores individuais para o nível de parcela.
Modelos de Relação Hipsométrica
Inicialmente, é necessário fazer a leitura dos dados que utilizaremos com exemplo para construção de modelos de relação hipsométrica no arquivo exemplo-hispo.csv.
hipso = read.csv("exemplo-hipso.csv",header=T,as.is=T) head(hipso)
Esses dados se referem a florestas plantads de Eucalyptus grandis. Somente as árvores com alturas medidas são apresentadas.
A relação altura-diâmetro para todo conjunto de dados mostra grande variação, embora tem uma relação média próxima da linear:
scatter.smooth( hipso$dap, hipso$ht, col="red")
Verifiquemos o número de árvores por parcela:
table( hipso$regiao, hipso$parcela )
O grande número de observações em algumas “parcelas” e regiões (Pilar do Sul, por exemplo) mostra que a variável parcela
não representa unicamente as parcelas. É necessário construir uma nova variável:
hipso$parc = paste(hipso$regiao, hipso$inv, hipso$faz, hipso$proj, hipso$talhao, hipso$parcela, sep="-") head(hipso) length( table(hipso$parc) ) hist( table(hipso$parc) )
Nota-se que com exceção de uma parcela, todas possuem mais de 10 árvores para ajuste da relação hipsométrica.
Análise Gráfica da Relação Hipsométrica
Para explorar melhor a relação altura-diâmetro num número grande de parcelas podemos utilizar o pacote gráfico lattice
com a função xyplot
:
require(lattice) xyplot( x = ht ~ dap, data=hipso)
A função xyplot
dois argumentos:
- Primeiro:
x
é uma fórmula semelhante a que se utiliza para definir um modelo linear na funçãolm
. - Segundo:
data
é o data frame contendo as variáveis.
A construção de gráficos através de fórmulas aumenta muito a flexibilidade na construção de gráficos.
Consideremos agora, a possibilidade de construir uma relação hipsométrica para cada região. Será que obteríamos bons modelos?
xyplot( ht ~ dap | regiao, data=hipso)
Note que nessa fórmula a barra vertical (|
) entra com a idéia de condição, isto é, construir um gráfico altura-diâmetro para cada região.
Ainda há uma grande variabilidade por região. Se a relação hipsométrica for construída para cada parcela, será que obteríamos bons modelos?
xyplot( ht ~ dap | parc, data=hipso)
Podemos tentar uma forma mais parsimoniosa de construção por região:
REG = unique( hipso$regiao ) REG xyplot( ht ~ dap | parc, data=hipso[ hipso$regiao == REG[1], ]) xyplot( ht ~ dap | parc, data=hipso[ hipso$regiao == REG[2], ]) xyplot( ht ~ dap | parc, data=hipso[ hipso$regiao == REG[3], ]) xyplot( ht ~ dap | parc, data=hipso[ hipso$regiao == REG[4], ])
Construção de Modelos de Relação Hipsométrica
A cada gráfico podemos fazer corresponder um modelo de relação hipsométrica.
Um modelo geral de relação hipsométrica é facilmente ajustado:
modgeral.lin = lm( ht ~ dap, data = hipso ) plot(modgeral.lin) summary(modgeral.lin)
O que os gráficos e o resumo falam da qualidade desse modelo?
Podemos nos perguntar se outros modelos, além do linear, apresentariam um desempenho melhor:
#Modelo Logarítmico modgeral.log = lm( log(ht) ~ log(dap), data = hipso ) plot(modgeral.log) summary(modgeral.log) #Modelo Schumacher modgeral.schu = lm( log(ht) ~ I(1/dap), data = hipso ) plot(modgeral.schu) summary(modgeral.schu) #Modelo Parabolóide modgeral.parab = lm( ht ~ dap + I(dap^2), data = hipso ) plot(modgeral.parab) summary(modgeral.parab)
Podemos também construir uma relação hipsométrica para cada região:
#Modelo Linear modreg.lin = lm( ht ~ dap * factor(regiao), data = hipso ) plot(modreg.lin) summary(modreg.lin) #Modelo Parabolóide modreg.parab = lm( ht ~ (dap + I(dap^2)) * factor(regiao), data = hipso ) plot(modreg.parab) summary(modreg.parab)
Houve melhora nos modelos, quando se modela no nível de região?
Podemos também construir uma relação hipsométrica para cada parcela:
#Modelo Linear modpar.lin = lm( ht ~ dap * factor(parc), data = hipso ) plot(modpar.lin) summary(modpar.lin) #Modelo Parabolóide modpar.parab = lm( ht ~ (dap + I(dap^2)) * factor(parc), data = hipso ) plot(modpar.parab) summary(modpar.parab)
Houve melhora nos modelos, quando se modela no nível de parcela?
Podemos fazer uma comparação de todos os modelos em termos de Coeficiente de Determinação (R2):
summary(modgeral.lin)$r.squared summary(modgeral.parab)$r.squared summary(modreg.lin)$r.squared summary(modreg.parab)$r.squared summary(modpar.lin)$r.squared summary(modpar.parab)$r.squared
Ou então uma comparação geral através do teste F, a partir do modelo mais simples até o modelo mais complexo:
anova(modgeral.lin, modreg.lin, modpar.lin) anova(modgeral.parab, modreg.parab, modpar.parab)
A comparação dos modelos também pode ser realizada em termos de AIC (Akaike Information Criterion):
AIC(modgeral.lin,modgeral.parab, modreg.lin,modreg.parab, modpar.lin,modpar.parab)
Ou de uma forma mais sofisticada, para construir uma tabela mais explicativa:
aic.tab = cbind(AIC(modgeral.lin, modreg.lin, modpar.lin), AIC(modgeral.parab, modreg.parab, modpar.parab)) rownames(aic.tab) = c("Geral","Regiao","Parcela") colnames(aic.tab) = paste(colnames(aic.tab), sort(rep(c("Linear","Parabolico"),2)) ) aic.tab
O que essa série de comparações nos diz?
Se verificarmos os resíduos dos modelos, veremos que os indicadores de desempenho dos modelos são coerentes:
#Verificando por Região boxplot( residuals(modgeral.parab) ~ regiao, data=hipso , main="Geral", ylim=c(-15,10)) abline(h=0,col="red") boxplot( residuals(modreg.parab) ~ regiao, data=hipso , main="Regiao", ylim=c(-15,10)) abline(h=0,col="red") boxplot( residuals(modpar.parab) ~ regiao, data=hipso , main="Parcela", ylim=c(-15,10)) abline(h=0,col="red") #Verificando por Talhão bwplot( residuals(modgeral.parab) ~ factor(talhao) | regiao , data=hipso , main="Geral", xlab="Talhão") bwplot( residuals(modreg.parab) ~ factor(talhao) | regiao , data=hipso , main="Região", xlab="Talhão") bwplot( residuals(modpar.parab) ~ factor(talhao) | regiao , data=hipso , main="Parcela", xlab="Talhão")
Exercício
Trabalhando nesses mesmos dados, verifique se a inclusão da variável idade
num modelo geral ou regional gera desempenho semelhante à modelagem por parcela.
Exercício
Trabalhando nesses mesmos dados, verifique se a inclusão da variável média das alturas dominates num modelo geral ou regional gera desempenho semelhante à modelagem por parcela. A média das alturas dominantes necessita ser calculada.
Modelos de Equação de Volume
Como exemplo para construção de equações de volume, utilizaremos os dados do arquivo exemplo-volume-cax.csv:
volcax = read.csv("exemplo-volume-cax.csv",header=T,as.is=T) head(volcax)
Nesses dados há duas informações de volume:
* volume comercial até o diâmetro de 7cm (v7cm
), e
* volume comercial até o diâmetro de 12cm (v12cm
), e
Há três variáveis nominais que identificam o caixetal onde foram cubadas as árvores:
regiao
: três regiões A (litoral Sul - SP), B (litoral Sul - RJ) e C (Vale do Ribeira);municip
: municípios dos caixetais;local
: código numérico de cada caixetal.
Comecemos com um modelo geral de equação de volume para o volume comercial até 7cm.
# Modelo de Spurr scatter.smooth( volcax$dap^2*volcax$ht, volcax$v7cm, col="grey20") v7.spurr = lm( v7cm ~ I(dap^2*ht), data=volcax ) plot(v7.spurr) summary(v7.spurr) scatter.smooth( volcax$dap^2*volcax$ht, volcax$v7cm, col="grey20") abline(coef(v7.spurr), col="green") # Modelo de Schumacher-Hall pairs( log(volcax[, c("dap","ht","v7cm")]) ) splom( log(volcax[, c("dap","ht","v7cm")]), panel = function(x,y){ panel.splom(x,y); panel.loess(x,y,col="red")} ) v7.sh = lm( log(v7cm) ~ log(dap) + log(ht), data=volcax ) plot(v7.sh) summary(v7.sh)
Os gráficos diagnósticos revelam que o modelo Spurr tem problema de heteroscedasticidade, enquanto que o modelo Schumacher-Hall, não tem esse problema e possui distribuição do resíduo mais próxima à distribuição Normal.
O coeficiente de determinação também sugere que o modelo Schumacher-Hall é melhor.
Equações de Volume com Variáveis Qualitativas
Uma questão importante é se podemos obter uma equação de volume melhor introduzirmos variávies qualitativas no modelo. Tomemos como exemplo o modelo Schumacher-Hall:
v7.sh.reg1 = lm( log(v7cm) ~ log(dap) + log(ht) + factor(regiao), data=volcax ) plot(v7.sh.reg1) summary(v7.sh.reg1) v7.sh.reg2 = lm( log(v7cm) ~ (log(dap) + log(ht)) * factor(regiao), data=volcax ) plot(v7.sh.reg2) summary(v7.sh.reg2) AIC(v7.sh, v7.sh.reg1, v7.sh.reg2) anova(v7.sh, v7.sh.reg1, v7.sh.reg2)
Qual dos modelos é mais apropriado?
Outro aspecto importante é saber se podemos dispensar a medição da altura se identificarmos a região do caixetal:
v7.sh.dr = lm( log(v7cm) ~ log(dap) * factor(regiao), data=volcax ) plot(v7.sh.dr) summary(v7.sh.dr) AIC(v7.sh.dr, v7.sh.reg1)
Qual a conclusão?
Exercício
Utilizando os dados de volume das árvores de caixeta, construa um modelo adequado para o volume comercial até 12cm (v12cm
).
Exercício
Verifique se o ajuste da equação de volume por caixetal (local
) torna dispensável a medição da altura.
Autor
João Luís Ferreira Batista
Laboratório de Biometria Ecológica
Centro de Métodos Quantitativos
Departamento de Ciências Florestais
Escola Superior de Agricultura "Luiz de Queiroz"
UNIVERSIDADE DE SÃO PAULO