[R] Add values of rlm coefficients to xyplot

[David Winsemius]

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