| [[http://cmq.esalq.usp.br|{{:publico:tree-of-life.gif?100|CMQ: Centro de Métodos Quantitativos}}]] ^ **__Centro de Métodos Quantitativos__** \\ Departamento de Ciências Florestais \\ Escola Superior de Agricultura "Luiz de Queiroz" \\ **UNIVERSIDADE DE SÃO PAULO** | [[publico:tutoriais:uso-r:start|{{:publico:tutoriais:r-relampago:r-logo.png?90}}]] |
^ |
Curso Básico de R
para Mensuração Florestal:
Modelos Lineares
^ |
| ^ Modelos Dendrométricos | |
\\
====== 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-hipso.csv|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ção ''lm''.
* 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|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 7//cm// (''v7cm''), e
* volume comercial até o diâmetro de 12//cm// (''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é 7//cm//.
# 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é 12//cm// (''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