Seleção de modelos - Alometria

• A princípio decidimos verificar se havia efeito do grupo (leg e não-leg). Essa diferença existe e esta demonstrada no script abaixo. Mas, depois mudamos de idéia (ver resultados e explicações a seguir)

Script seleção de modelos GRUPO

• No final, decidimos agrupar leg e não leg (apesar do efeito de grupo encontrado) porque para comparar temos que ter uma área igual no entorno de leg e não-leg!!!!
• Os três modelos contrastados foram:
  1. Raio copa não varia em função do DAP
  2. Raio copa varia linearmente em função do DAP
  3. Raio copa varia exponencialmente em função do DAP
• Analisando os valores de AIC, o modelo mais plausível foi o 3, e a equação encontrada nele foi: 0.047*(DAP)^0.8.
• Script dessa seleção:

Seleção de modelos DAP x COPA

• Amostra estratificada e balanceada
• Dados de DAP e altura
• Quatro medidas de raio copa, em cruz a partir do tronco, tirei média depois.

Dados brutos finais

• Coletei no campo dados de diâmetro maior e diâmetro menor da copa (cm)
• Calculei a área considerando uma elipse → pi * (DMaior/2) * (DMenor/2)
• Dúvida: chego no raio considerando que essa área é de uma circunferência?
• Acho que to errando nas transformações, pois gero raios de copa gigantes…

Dados brutos piloto

Escolhi as árvores estratificadas, entretanto, em campo algumas não leg estavam mortas, outras emaranhadas com outra copa impossibilitando a visualização e algumas não foram encontradas. Algumas árvores estavam em parcelas que ainda não foram recenseadas e o acesso estava muito difícil. Acima de 450 de dap, a grande maioria é balizia e do outro grupo há alguns Calophyllum e Ocotea.

Testando efeito do grupo

1. Modelo 1 = copa ~ dap
2. Modelo 2 = copa ~ dap + grupo
3. Modelo 3 = copa ~ dap * grupo

Comparando com anova e AIC, há interação…

1. Modelo linear para legs
   • Regressão simples: Raio.copa.leg = 0.49 + 0.01*(Dap.leg)… significativa, com r² = 0.9.
2. Modelos para não legs
   • Regressão simples: Raio.copa.nleg = 1.3 + 0.004* Dap.nleg… significativa com r² = 0.5
   • M-M: Raio.copa.nleg = 2.92*Dap.nleg/(56.6 + Dap.nleg)

Os AICs foram similares, fiquei com a regressão simples. Além disso, o erro da parâmetro b na MM foi enorme (26!).

Arquivo RData com modelos de alometria julho 09

Script 1ª análise(leg + nleg) julho09

Script 2ª análise (efeito do grupo) julho09

plot(rm.copa~dap, col=c("red", "blue")[as.factor(copa$Grupo)],pch=19, xlab= "DAP (cm)", ylab= "raio da copa (m)")
## rm.copa é o raio em metros da copa e o dap é o dap em cm. Copa$grupo é a coluna de dados com o  fator leg x nleg 

 locator(1)->locxy
## aqui vc clica  no ponto xy no gráfico onde deseja colocar a legenda e guarda no objeto locxy.

legend(locxy, legend=c("leguminosae", "outras familias"), col=c("red", "blue"), pch=c(19,19))
## inclui a legenda na posição locxy

Como eram poucos dados plotei todos leg e não leg juntas… \

\

Comparei 5 modelos: \ 1)O mais simples :

copa.vero1<-function(m=2.542,s=1.96)
	{
     -sum(stats::dnorm(alo$Raio.copa.m, m, s, log=TRUE))
	}
copa1<-mle(copa.vero1)

2) Com a função original:

copa.vero2<-function(a=0,b=1,s=1.96)
	{
	m = a+b*((0.1*alo$Dap.mm)^0.5/2)
	-sum(stats::dnorm(x=alo$Raio.copa.m,mean=m,sd=s,log=T))
	}

copa2←mle(copa.vero2) 3) Com outra função, sem elevar a 0.5:

copa.vero3<-function(a=0,b=1,s=1.96)
	{
	m = a+b*((0.1*alo$Dap.mm)/2)
	-sum(stats::dnorm(x=alo$Raio.copa.m,mean=m,sd=s,log=T))
	}

copa3←mle(copa.vero3) 4) Elevando ao quadrado sem dividir:

copa.vero4<-function(a=0,b=1,s=1.96)
	{
	m = a+b*(0.1*alo$Dap.mm)^2
	-sum(stats::dnorm(x=alo$Raio.copa.m,mean=m,sd=s,log=T))
	}

copa4←mle(copa.vero4)

5) Com média e desvio variando pela mesma função:

copa.vero5<-function(a=1.43,b=0.002,s1=1,s2=2)
	{
	m = a+b*((0.1*alo$Dap.mm)/2)
	s = s1+s2*((0.1*alo$Dap.mm)/2)
	-sum(stats::dnorm(x=alo$Raio.copa.m,mean=m,sd=s,log=T))
	}

copa5←mle(copa.vero5)

* Como a curiosidade é grande rodei alguns nulos alterando a função r.copa para o modelo 3 e o resultado foi:

• Para Abarbr, 259 médias ficaram inferiores ao valor observado. Com a função anterior esse valor era 0.
• Balipe continuo com 0, como possui bastante dbh grande acredito que a função possa estar errando aí…
• Ormoar teve 8 valores inferiores, com a função anterior era 1 (sem muita diferença então)

Continuo não entendendo porque os modelos tem parâmetros como “0.01*” e “/2”. Não fazem o minimo sentido. Esses são os parâmetros que vc. quer estimar… eles estão incluídos no “a” e “b” da regressão linear…

Acho que devemos nos limitar nesse momento apenas ao lm. Minhas sugestões talvez não tenham ficado claras.

1. passo: Confrontar o modelo linear com o fator “leg” e “nleg”, com o modelo linear se esse fator, se não me engano vc. chamou “Grupo”.

• mod1=lm(copa~dap*grupo)
• mod2=lm(copa~dap)

use: summary(); anova() para ver os modelos separadamente e depois anova ou AICtab para confrontá-los. Se no modelo um há interação significa que leg e nleg devem ser tratados diferentemente na hora de estimar a copa. Caso o modelo só com o dap seja melhor, então podemos fazer um modelo único.

2.passo: decidido qual o modelo base, devemos estimar os parâmetros (leg/nleg junto ou isolado dependendo do passo anterior)

• modelo linear simples e com desvio padrão variando

 copa~a+ b* dap

3. estimado a e b voltamos ao modelo nulo para ajustá-lo

• mesmo para um piloto os dados estão mal coletados. Não há a distribuição de diâmetros estratificada como combinamos e tb. não há correspondência leg x nleg, como vc. afirmou! Para evitar essas surpresas deve olhar os dados ainda quando está em campo para fechar lacunas de coleta. Com esses dados não há como chegar a um modelo legal de predição, mas ao menos podemos fazer um ensaio para esse relatório. Assim os códigos ficam prontos para serem usados quando vc. coletar mais dados.

• a partir desse relatórios não estamos mais na fase de pilotos e tentivas. Agora é vida real e vc. precisa se conscientizar que os prazos estão contra vc…. infelizmente não há mais tempo para ficarmos contemplando e pirando em modelos. Também não há mais tempo para vc. ficar quebrando a cabeça com códigos. QQ travada que der deve me procurar imediatamente para resolvermos.

• até agora há mais esforço meu que seu no trabalho intelectual de definição dos modelos nulos. Eu gosto da brincadeira, mas o brinquedo não é meu! A partir de agora vc. assume o brinquedo e eu apenas facilito sua vida com os códigos e ajudo apontando caminhos. A decisão de até onde sua monografia irá chegar é sua.

• nesse sentido estabeleço o que julto ser o mínimo para ser atingido (ao final da tese):
  1. alometria da copa x diâmetro para melhorar o modelo nulo (n min: 100 arvores legs e 100 nleg), estratificado e balanceado
  2. testar o modelo nulo de índices para todas as legs (usando os dados do recenso!)
  3. testar para jovens com DAP<1cm leg x nleg (levantar os jovens em campo), n mínimo de parcelas: 30 abaixo de legs e 30 abaixo de nlegs

• até o dia da entrega do seu relatório vc. deve atualizar o que fez diariamente no wiki. O que vamos fazer é:
  1. seleção do modelo de predição de raio de copa como os dados que tem.
  2. ajustar o modelo nulo com essa alometria e rodar para todas as legs.
  3. as análises devem ser feitas essa semana para vc. poder escrever o relatório na semana que vem.

• depois do seu relatório vc. deve atualizar mensalmente o wiki e prever um boa e longa estada em campo para:

1. levantar os dados de alometria copa x dap
2. levantar dados de jovens dap< 1cm abaixo de leg e nleg
3. agilizar a entrada de dados da parcela permanente para estar apta a usar os dados ao final do seu trabalho

• durante esse tempo deve focar em leitura específica sobre o seu trabalho e mais ainda sobre no artigo que irá apresentar na sua qualificação. Quando terminar o relatório decidimos qual será o artigo…