[R-br] ANOVA em regressão não linear!

Thiago de Paula Protásio depaulaprotasio em yahoo.com.br
Segunda Janeiro 16 19:27:37 BRST 2012


 Walmes e demais membros da lista;


Vejam até onde consegui desenvolver com base no exemplo proposto na página
das "Rídiculas".
Tentei aplicar a função do Walmes para obter a anova do modelo não linear e
o R², mas não obtive sucesso.
Desde já agradeço a ajuda de todos!



A <- 6.25551
> B <- -6.85223
> C <- -0.24253


> model8 <- nls(log(Carbtotal)~A+B*exp(C*Idade), data=citriodora,
+           start=list(A=A, B=B, C=C))
> summary(model8)

Formula: log(Carbtotal) ~ A + B * exp(C * Idade)

Parameters:
  Estimate Std. Error t value Pr(>|t|)
A  6.25552    0.26527   23.58 1.14e-12 ***
B -6.85224    0.27415  -25.00 5.14e-13 ***
C -0.24253    0.02415  -10.04 8.88e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2841 on 14 degrees of freedom

Number of iterations to convergence: 1
Achieved convergence tolerance: 2.429e-06

> F <- attr(model8$m$fitted(), "gradient")  ## extraindo a matriz gradiente
avaliada nas estimativas dos parâmetros

> m0 <- lm(log(Carbtotal)~-1+F, data=citriodora)# passando a matriz
gradiente para a lm(), importante remover intercepto
> anova(m0) # partição da soma de quadrados total


Analysis of Variance Table

Response: log(Carbtotal)
          Df Sum Sq Mean Sq F value    Pr(>F)
F          3 271.26  90.422  1119.9 < 2.2e-16 ***
Residuals 14   1.13   0.081
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(model8, lm(log(Carbtotal)~1, citriodora)) # SQ do modelo não linear
corrigido para intercepto
Analysis of Variance Table

Model 1: log(Carbtotal) ~ A + B * exp(C * Idade)
Model 2: log(Carbtotal) ~ 1
  Res.Df Res.Sum Sq Df  Sum Sq F value    Pr(>F)
1     14      1.130
2     16     82.335 -2 -81.204  502.88 9.192e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


> R2 <- function(nls.obj){
+   da <- eval(nls.obj$data)
+   resp.name <- all.vars(summary(nls.obj)$formula)[1]
+   form <- paste(resp.name, "~1", sep="")
+   m0 <- lm(form, da)
+   an <- anova(nls.obj, m0)
+   sqn <- deviance(nls.obj)
+   sqe <- deviance(m0)
+   r2 <- 1-(sqn/sqe)
+   aov <- data.frame(fv=c("regression","residuals"),
+                     gl=c(-an$Df[2],an$Res.Df[1]),
+                     sq=c(-an$Sum[2],an$Res.Sum[1]))
+   aov$qm <- aov$sq/aov$gl
+   aov$F <- c(aov$qm[1]/aov$qm[2], NA)
+   aov$"Pr(>F)" <- c(1-pf(aov$F[1], df1=aov$gl[1], df2=aov$gl[2]), NA)
+   names(aov) <- c(" ","Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")
+   return(list(anova=aov, R2=r2))
+ }


> R2(model8)


Erro em anovalist.nls(object, ...) :
  'anova' is only defined for sequences of "nls" objects
Além disso: Warning message:
In anovalist.nls(object, ...) :
  models with response "Carbtotal" removed because response differs from
model 1


-- 
Thiago de Paula Protásio
Acadêmico de Engenharia Florestal
Universidade Federal de Lavras
Departamento de Ciências Florestais
Ciência e Tecnologia da Madeira
(035) 9183 - 2246
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120116/fa10546c/attachment.html>


Mais detalhes sobre a lista de discussão R-br