====== Modelo Linear Misto - lme ====== ====== Análise de Dados ====== 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 ^ ^ Model ^ df ^ AIC ^ BIC ^ logLik ^ Test ^ L.Ratio ^ p-value ^ ^lme.all | 1 | 16 | 265.9761 | 318.9874 | -116.9881 | ^lm.decomp | 2 | 14 | 263.3403 | 309.7252 | -117.6702 | 1 vs 2 | 1.364197 | 0.5056 | ## 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) ==== k para modelo não linear ==== 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 ==== Gráfico base ==== 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) {{graf1.jpg?700|}} ==== GRÁFICO COM OS TRATAMENTOS POR FLORESTA ==== coplot(massa.rem/100~dias| amb+fauna, data=decomp[decomp$flor== unique(decomp$flor)[1],]) **Floresta Estacional Semidecidual** {{graf.estac.jpg?700|}} coplot(massa.rem/100~dias| amb+fauna, data=decomp[decomp$flor== unique(decomp$flor)[2],]) // // **Floresta Altântica - PECB** {{graf.atlantica.jpg?700|}} coplot(massa.rem/100~dias| amb+fauna, data=decomp[decomp$flor== unique(decomp$flor)[3],]) **Floresta de Restinga- PEIC** {{graf.rest.jpg?700|}} ==== Graficos pelos tratamento por ambiente==== - 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 ** {{graf.abg.jpg?700|}} ** Ambiente de decomposição de raizes ** {{graf.belg.jpg?700|}} coplot(massa.rem/100~dias| flor+fauna, data=decomp[decomp$amb== unique(decomp$amb)[1],]) ==== CALCULAR O K PARA cada nivel de tratamento ==== **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) ===== Resultado do Modelo Misto ===== 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) ^ ^numDF^ denDF ^ F-value ^ p-value^ |(Intercept)| 1 | 202 |13267.795| <.0001| |dias | 1 | 202 |73.801| <.0001| |as.factor(flor)| 2 | 202 |7.804| 0.0005| |as.factor(amb)| 1 | 202 |6.602| 0.0109| |as.factor(fauna)| 1 | 202 |32.659| <.0001| |as.factor(flor):as.factor(amb)| 2 | 202 |1.523| 0.2205| |as.factor(flor):as.factor(fauna) | 2 | 202 |3.477 |0.0328| |as.factor(amb):as.factor(fauna) | 1 | 202 |6.399 |0.0122| |as.factor(flor):as.factor(amb):as.factor(fauna) | 2 | 202 |1.497 |0.2264|