On Nov 9, 2010, at 8:01 AM, PtitBleu wrote:
Hello again,
Thanks to this forum, I found a lapply command to apply "rlm" to
each Device
of the Name column of the df1 data frame (df1 data frame is given in
the
previous post).
I'm now looking for a nice way to create a data frame with the name
of the
device, the intercept and the slope. I now use a for loop (see
below) but
I'm sure there is a more elegant way to retrieve these coefficients,
isn't
it ?
Thanks in advance for your help to improve this script.
Ptit Bleu.
# Library including rlm
library(MASS)
# To apply rlm to each Device of df1$Name
bbb<-lapply(split(df1, df1$Name), function(X)(rlm(col2 ~ col3,
data=X)})
# To get Intercept and Slope for each device
dfrlm1<-data.frame()
for (u in 1:length(names(bbb))) {
dfrlm1<-rbind(dfrlm1, data.frame(as.character(names(bbb)[u]),
bbb[[u]][[1]][[1]], bbb[[u]][[1]][[2]]))
}
names(dfrlm1)<-c("Name", "Intercept", "Slope")
That looks a bit painful. It might be easier to simple wrap the rlm
call in the lapply(split(...)) step with coef() and then you would
have just the coefficients in a neat bundle. Generally that list can
be made into a data.frame or matrix with do.call("rbind", list-
object), possibly needing as.data.frame to finish the bundling process.
?coef # which is a function that is a method for many regression fit
classes
methods(coef)
> class(rlm(stack.loss ~ ., stackloss))
[1] "rlm" "lm"
Then, of course, there are the plyr and data.table packages that often
make the process even beautiful.
--
David
# To plot the graphs and display the rlm coefficients
x11(15,12)
xyplot(df1$col2 ~ df1$col3 | df1$Name,
panel = function(x, y,...) {
panel.abline(h=seq(20,30,1), col="gray")
panel.abline(v=seq(20,70,5), col="gray")
panel.xyplot(x, y, type="p", col="red", pch=20,...)
panel.abline(rlm(y ~ x), col="blue")
panel.text(40,21, paste("y =",
round(dfrlm1$Intercept[panel.number()],3),"+ (",
round(dfrlm1$Slope[panel.number()],3),") *x"))
}, scales=list(cex=1.2),
xlab=list("X", cex=1.4), ylab=list("Y", cex=1.4),
xlim=c(20,70), ylim=c(20,30))
--
David Winsemius, MD
West Hartford, CT
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
