10 de novembro 2009

Estou fazendo as análises com modelos mistos, utilizando o pacote nlms no R. Preciso definir quais os fatores randomicos do modelo. A principio são os dias/blocos

lme(logmassa ~ dias + flor * amb * fauna, random= ~1|bloco)->lme.bloco 

##modelo total levando em conta o fator randomico bloco

lm(lgmassa ~ dias + flor * amb * fauna, data=decomp)->lm.decomp ## modelo simples

anova(lme.all,lm.decomp)

RESULTADO

## não há acrescimo de informação pelo incorporação do bloco. Cuidado! o bloco é um fator que na realidade só diminiui o erro e aumenta o poder do teste.

## as análises deverão levar em consideração o bloco, entretanto os graficos serão feitos com o decaimento do modelo (k)

nls.k1=nls(I(massa.rem/100)~a+exp(k*dias),start=list(a=0.1, k=-.004),data=decomp)

é preciso incluir um valor de iniciação do k, aqui o modelo prevê uma estimativa de intercepto, parece melhor!

summary(nls.k1)

Formula: I(massa.rem/100) ~ a + exp(k * dias)
Parameters:
^^Estimate ^ Std. Error ^t^ value Pr(>|t|)^    
| a | -0.0550089 |  0.0161366 |  -3.409 0.000779|  
| k | -0.0046770 |  0.0002649 | -17.655 |  < 2e-16|  
 ---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
Residual standard error: 0.1134 on 214 degrees of freedom
Number of iterations to convergence: 8 
Achieved convergence tolerance: 5.045e-06 

plot(massa.rem/100~dias, data=decomp)
exp.fit=function(x,a,k){a+exp(k*x)}
curve(exp.fit(x,-0.0550089,-0.0046770), from=0, to=300, col="red",lty=2, add=TRUE)

pode ser feito direto:

curve(a+exp(k*x), from=0, to=300, col="red",lty=2, add=TRUE)

coplot(massa.rem/100~dias| amb+fauna, data=decomp[decomp$flor== unique(decomp$flor)[1],])

Floresta Estacional Semidecidual

coplot(massa.rem/100~dias| amb+fauna, data=decomp[decomp$flor== unique(decomp$flor)[2],])

Floresta Altântica - PECB

coplot(massa.rem/100~dias| amb+fauna, data=decomp[decomp$flor== unique(decomp$flor)[3],])

Floresta de Restinga- PEIC

1. Fator Ambiente de Decomposição

coplot(massa.rem/100~dias| flor+fauna, data=decomp[decomp$amb== unique(decomp$amb)[1],])

Ambiente de decomposição serrapilheira

Ambiente de decomposição de raizes

coplot(massa.rem/100~dias| flor+fauna, data=decomp[decomp$amb== unique(decomp$amb)[1],])

Função do R

k.calc=function(dados=decomp)
{
library(nlme)
unico.flor=unique(dados$flor)
unico.bloco=unique(dados$bloco)
unico.amb=unique(dados$amb)
unico.fauna=unique(dados$fauna)
##resulta=matrix(NA,ncol=10, nrow=24)
floresta =ambiente=exclfauna= blocorand= coefa = coefk= a5 =k5= a95= k95=r2=rep(NA,24)
u=0
   for (f in 1: length(unico.flor))
  {
  flor1=dados[dados$flor==unico.flor[f],]
  for (j in 1: length(unico.amb))
  {
     amb1=flor1[flor1$amb==unico.amb[j],]
      for (h in 1: length(unico.fauna))
      {
      fauna1=amb1[amb1$fauna==unico.fauna[h],]
          for(i in 1:length(unico.bloco))
          {
          u=u+1
          bloco1=fauna1[fauna1$bloco==unico.bloco[i],]
          nls.k1=nls(I(massa.rem/100)~a*exp(k*dias),start=list(a=1, k=-.004),data=bloco1)
          lm.k1=lm(I(massa.rem/100)~1,data=bloco1)
          ssres.nls=sum(resid(nls.k1)^2)
          ssres.lm=sum(resid(lm.k1)^2)
          r2[u]=(ssres.lm-ssres.nls)/ssres.lm
          floresta[u]=unico.flor[f]
          ambiente[u]=unico.amb[j]
          exclfauna[u]=unico.fauna[h]
          blocorand[u]=unico.bloco[i]
          coefa[u]=coef(nls.k1)[1]
          coefk[u]=coef(nls.k1)[2]
          a5=confint(nls.k1)[1]
          k5=confint(nls.k1)[2]
          a95=confint(nls.k1)[3]
          k95=confint(nls.k1)[4]
          }
       }
     }   
  }
 resulta=data.frame(cbind(floresta,ambiente,exclfauna,blocorand,r2,coefa,coefk,a5=as.numeric(a5),k5=as.numeric(k5), a95=as.numeric(a95),k95=as.numeric(k95))) 
return(resulta)
}

Para calcular agora é só usar a função

k.calc(decomp)->decomp.k

Estou tendo problemas com a função, as variáveis numerica estão saindo como fatores (10 novembro de 2008)

Deve antes de iniciar a modelagem baixa o pacote nlme

library(nlme)

Modelando todas os fatores na decomposição:

lme.decomp=lme(lgmassa~dias+as.factor(flor)*as.factor(amb)*as.factor(fauna),random=~dias|bloco, data=decomp)
anova(lme.decomp)