I was going through all kinds of torture to figure this out. An elegant solution, indeed.
Thank you very much. Shawn Way, PE Engineering Manager [EMAIL PROTECTED] -----Original Message----- From: Liaw, Andy [mailto:[EMAIL PROTECTED] Sent: Friday, April 16, 2004 9:15 AM To: 'Shawn Way' Subject: RE: [R] residuals Shawn, This is what I'd do: lof <- function(x, y) { fit1 <- anova(lm(y ~ x)) fit2 <- summary(aov(y ~ factor(x)))[[1]] SSpe <- fit2[["Sum Sq"]][2] DFpe <- fit2[["Df"]][2] SSresid <- fit1[["Sum Sq"]][2] DFresid <- fit1[["Df"]][2] SSlof <- SSresid - SSpe DFlof <- DFresid - DFpe Fstat <- (SSlof / DFlof) / (SSpe / DFpe) pval <- pf(Fstat, DFlof, DFpe, lower.tail=FALSE) out <- matrix(c(SSlof, SSpe, DFlof, DFpe, Fstat, NA, pval, NA), nrow=2, dimnames=list(c("Lack-of-fit", "Pure Error"), c("Sum Sq", "Df", "F value", "Pr(>F)"))) class(out) <- c("anova", class(out)) out } On your data, I get: > lof(data$ref, data$actual) Sum Sq Df F value Pr(>F) Lack-of-fit 0.13441 1 0.7984 0.4374 Pure Error 0.50505 3 Best, Andy > From: Shawn Way > > > I'm trying to determine the lack of fit for regression on the > following: > > data <- data.frame(ref=c(0,50,100,0,50,100), > actual=c(.01,50.9,100.2,.02,49.9,100.1), > level=gl(3,1)) > fit <- lm(actual~ref,data) > fit.aov <- aov(actual~ref+Error(level),data) > > According to the information I have, the lack of fit for this > regression is > the f-ratio between the residuals of the level contribution and the > residuals within. I'm trying to get the information from the > fit to make > this ratio of the residual mean squares. > > Any thoughts? > > ______________________________________________ > [EMAIL PROTECTED] mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > > ---------------------------------------------------------------------------- -- Notice: This e-mail message, together with any attachments,...{{dropped}} ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
