1

I have loaded the data set Animals2 from the package library (robustbase), and I am interested on working with the logarithm of these data.

library(robustbase)
x<-log(Animals2)
plot(x, main="Plot Animals2", col="darkgreen")

Now I need to add the least-squares line, the LMS line, and the line of Siegel, by using different colours. Finally, I also wanted to show two plots with the estimated 2-dimensional densities. The first plot should be a 3D visualization of the density (I was trying to use command persp) and the second plot with the command image (applied to the density estimate).

I am a bit stuck and would appreciate some help.

nils
  • 25
  • 5
  • `abline()` will let you plot linear model lines in the scatterplot. The `kde2d()` function from the `MASS` package calculates 2-D density estimates and you should then be able to use `persp()` and `image()` to make the plots. – DaveArmstrong May 20 '22 at 10:27
  • @DaveArmstrong I am having trouble with the line of Siegel. I am using the following code but it doesnt work: `install.packages("mblm") library("mblm") mblm(Animals2$body ~ Animals2$brain , repeated=TRUE)` – nils May 20 '22 at 10:36
  • You should be able to do it with `mblm(body ~ brain, data=x, repeated=TRUE)`. Just beware that when you do `plot(x)` it puts `brain` on the y-axis and `body` on the x-axis. In your model, you've got it the other way around. If you want `brain` to be the dependent variable, use `mblm(brain ~ body, ...)` otherwise, you should use `plot(body ~ brain, data=x)` to put `body` rather than `brain` on the y-axis. – DaveArmstrong May 20 '22 at 14:40
  • @DaveArmstrong Okay thank you very much I understand now. And which of the 3 the regression lines would you say is more suitable for the data? – nils May 22 '22 at 16:25
  • The Least Squares line is drawn toward the cluster of three outliers. To me, this makes it less suitable, but the other lines are discounting (down-weighting) those points, so the inferential qualities of those lines may be somewhat dubious. The question of which line is most suitable is one that would be better posed to [CrossValidated](https://stats.stackexchange.com). – DaveArmstrong May 22 '22 at 17:44

1 Answers1

1

Assuming that mblm() does the Siegel model, you could use this:

  library(robustbase)
  library(mblm)
  data(Animals2)
  x <- log(Animals2)
  plot(x, main="Plot Animals2", col="darkgreen")
  abline(lm(brain ~ body, data=x), col="red")
  abline(MASS::lqs(brain ~ body, data=x, method="lms"), col="blue")
  abline(mblm(brain ~ body, data=x, repeated=TRUE), col="black")
  legend("topleft", c("LM", "LMS", "Siegel"), 
         col = c("red", "blue", "black"), 
         lty = c(1,1,1), 
         inset=.01)



d2 <- MASS::kde2d(x$body, x$brain)
persp(d2)

image(d2, xlab="body", ylab="brain")

Created on 2022-05-20 by the reprex package (v2.0.1)

DaveArmstrong
  • 18,377
  • 2
  • 13
  • 25