16 Jan
2012
16 Jan
'12
21:27
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