Ferramentas do usuário

Ferramentas do site


Barra lateral

CMQ
Centro de Métodos Quantitativos


USP ESALQ
Depto. de Ciências Florestais
ESALQ
UNIVERSIDADE de SÃO PAULO
Av. Pádua Dias, 11
Caixa Postal 09
13418-900 - Piracicaba - SP
BRASIL
biometria:tutoriais:r-basico-mensuracao:modelos
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

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:

  1. ajustar modelos dendrométricos utilizando regressão linear;
  2. verificar a qualidade das predições dos modelos ajustados e selecionar o modelo mais adequado;
  3. gerar predições a partir do modelo selecionado para dados de árvores individuais;
  4. 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çã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:

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

biometria/tutoriais/r-basico-mensuracao/modelos.txt · Última modificação: 2015/08/10 20:48 (edição externa)