Dear Pep
To choose between different models you can use anova() (i.e.,
anova.lmrob()) which runs robust tests; cf. example below. Since
anova.lmrob() can not handle poly(x, ...) directly, you must run some
manual preparations on the data set (cf. below).
Fitting higher order polynomial can be very tricky. If there are
outliers at one or both end of the range, it is hard to tell whether one
should fit a higher order polynomial or better a low order polynomial
with some outliers at the end of the range. - The goals of your analysis
may give a hint which solution is more preferred. The data themselves
are often not decisive.
If there is no physical law supporting a high order polynomial I prefer
to fit local regression smoother like lowess() or loess.
To plot the models you must install the current version of the package
'robustbase' (cf. below) where predict.lmrob() is implemented.
All the best, Andreas
--------------------- Example ----------------------------------
ADat <- data.frame(x=1:19)
set.seed(3011)
ADat$y <- 4 + 0.02*(ADat$x - 9)^3 + rnorm(nrow(ADat),0,0.7)
ADat$y[c(5,12)] <- c(15,-5)
## Data Preparation:
X <- model.matrix(~ poly (x,6), data=ADat)
ADat2 <- data.frame(y=ADat$y, X[,-1])
names(ADat2)[-1] <- paste("P", 1:6, sep="")
## Robust Fitting
library(robustbase)
pl1 <- lmrob(y ~ P1, data=ADat2)
pl2 <- lmrob(y ~ P1 + P2, data=ADat2)
pl3 <- lmrob(y ~ P1 + P2 + P3, data=ADat2)
pl4 <- lmrob(y ~ . - P5 - P6, data=ADat2)
pl5 <- lmrob(y ~ . - P6, data=ADat2)
pl6 <- lmrob(y ~ ., data=ADat2)
summary(pl6)
## Model Selection
anova(pl6, pl3, test="Wald")
## or alternatively
anova(pl6, pl3, test="Deviance")## Pr(>chisq) = 0.3652
anova(pl6, pl2, test="Deviance")## Pr(>chisq) = < 2.2e-16
## Displaying Model Fits
## For predict: Please load 'robustbase' from R-forge
## see https://r-forge.r-project.org/R/?group_id=59
plot(y ~x, data=ADat)
lines(ADat$x, predict(pl1), col=2)
lines(ADat$x, predict(pl2), col=3)
lines(ADat$x, predict(pl3), col=4)
lines(ADat$x, predict(pl4), col=5)
----------------------------------------------------------------------------------
------- Urspr?ngliche Nachricht --------
Von: Pep Serra <josep.serra at uab.cat>
Gesendet: Mon, 29 Nov 2010 15:52:08 +0100
An: r-sig-robust at r-project.org
Betreff: [RsR] choosing and plotting a model
Dear all,
Trying to choose between different curve fittings within different
robust models (what we would get with cftools in MATLAB, but I am
definetly R-addicted).
Maybe it is a trivial questions:
lm1<-lmrob(y ~ x, data=a)
pl2<-lmrob(y ~ poly (x,2), data=a)
pl3<-lmrob(y ~ poly (x,3), data=a)
pl4<-lmrob(y ~ poly (x,4), data=a)
pl5<-lmrob(y ~ poly (x,5), data=a)
pl6<-lmrob(y ~ poly (x,6), data=a)
I want to choose between different models according to the best fit, but
the problem is that there is no "best fit" in robust (or R2). How could
I do that?
In addition, I would love to plot these models (scatterplot and than
later add the lines to see how the models look like)
I tried, for instance:
lines (x,predict (pl2),col="red")
but I get many curves (I think) or a filled polygon. Is there any way to
get a simple line out of the model?