0

Load a data set comprising spectral intensities of 60 samples of gasoline at 401 wavelengths, and their octane ratings in pls package. In my example:

library(pls)

#Data set

data(gasoline)

'data.frame':   60 obs. of  2 variables:
 $ octane: num  85.3 85.2 88.5 83.4 87.9 ...
 $ NIR   : 'AsIs' num [1:60, 1:401] -0.0502 -0.0442 -0.0469 -0.0467 -0.0509 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr  "1" "2" "3" "4" ...
  .. ..$ : chr  "900 nm" "902 nm" "904 nm" "906 nm"

I've like to plot the 3D representation of NIR in function of octane, my holy grail plot is:

https://it.mathworks.com/help/stats/examples/partial-least-squares-regression-and-principal-components-regression.html

But this output plot was create in matlab and not in R. The original code is (https://it.mathworks.com/help/stats/examples/partial-least-squares-regression-and-principal-components-regression.html):

load spectra
whos NIR octane

[dummy,h] = sort(octane);
oldorder = get(gcf,'DefaultAxesColorOrder');
set(gcf,'DefaultAxesColorOrder',jet(60));
plot3(repmat(1:401,60,1)',repmat(octane(h),1,401)',NIR(h,:)');
set(gcf,'DefaultAxesColorOrder',oldorder);
xlabel('Wavelength Index'); ylabel('Octane'); axis('tight');
grid on

Please any ideas or similar function/packages in R? Thanks in advance!

user2554330
  • 37,248
  • 4
  • 43
  • 90
Leprechault
  • 1,531
  • 12
  • 28
  • You can do a plot like that in `rgl` or `plot3D`. If you don't see how, please post some data and we'll show you. – user2554330 Apr 20 '19 at 14:55
  • I try below, but not so good presentation like MATLAB: `library(rgl)` `library(pls)` `data(gasoline)` `str(gasoline)` `gasoline$ID<-seq(1:length(gasoline[,1]))` `open3d()` `x <- gasoline$octane` `y <- gasoline$NIR` `z <- gasoline$D` `plot3d(x, y, z, , type = "l", col = rainbow(1000))` – Leprechault Apr 20 '19 at 17:40
  • 1
    this could also be done in base graphics by using `persp()` to draw an empty plot, *saving the returned object*, and using the `$trans3d()` function it contains to add the lines (see the `Value:` section in `?persp`) – Ben Bolker Apr 21 '19 at 01:03

1 Answers1

2

This gives something like the Matlab plot:

library(rgl) 
library(pls) 
data(gasoline) 
str(gasoline) 
gasoline$ID <- seq(1:length(gasoline[,1])) 
open3d() 
x <- gasoline$octane 
y <- gasoline$NIR 
z <- gasoline$ID 
plot3d(x, 1:ncol(y), y, type = "n", xlab = "", ylab = "", zlab = "")
cols <- rainbow(1000)[(x - min(x))/(max(x) - min(x))*999 + 1]
for (i in seq_along(x))
  lines3d(x[i], 1:ncol(y), y[i,], col = cols[i])

screenshot

There are several differences:

  • The handedness of the coordinate system. (Octane increases towards the viewer in the rgl plot.) You can change it in rgl using par3d("userMatrix" = par3d("userMatrix") %*% diag(c(-1, 1,1,1))) if you really want the Matlab style.
  • I didn't draw axis labels. The default ones in rgl are kind of ugly, and I was too lazy to use mtext3d to draw nicer ones.
  • I didn't bother with background grids. I don't like them, but if you really want them, you can get them with grid3d.
  • Matlab doesn't show perspective. If that's what you want, use par3d(FOV = 0).
  • The aspect ratios are different. I think Matlab is using something like aspect3d(1, 1, 0.75).
  • The colors. rgl uses the R color system, so if you can find a palette you like better than rainbow(1000), use that.
  • The price.
user2554330
  • 37,248
  • 4
  • 43
  • 90