In luca-scr/msir: Model-Based Sliced Inverse Regression

library(knitr)
opts_chunk$set(fig.align = "center", 
               out.width = "70%",
               fig.width = 5, fig.height = 4,
               dev.args=list(pointsize=10),
               par = TRUE, # needed for setting hook 
               collapse = TRUE, # collapse input & ouput code in chunks
               warning = FALSE)
knit_hooks$set(par = function(before, options, envir)
  { if(before && options$fig.show != "none") 
       par(mar=c(4.1,4.1,1.1,1.1), mgp=c(3,1,0), tcl=-0.5)
})

# setupKnitr()
# knit_hooks$set(rgl = hook_webgl)

set.seed(1) # for exact reproducibility

Introduction

Model-based Sliced Inverse Regression (MSIR) is a dimension reduction method based on Gaussian finite mixture models which provides an extension to sliced inverse regression (SIR).

The basis of the MSIR subspace is estimated by modeling the inverse distribution within slice using Gaussian finite mixtures with number of components and covariance matrix parameterization selected by BIC or defined by the user.

The msir package implements the methodology described in Scrucca (2011).

This vignette is written in R Markdown using the knitr package for production.

library(msir)

Example: 1-dimensional nonlinear response curve

n <- 200
p <- 5
b <- as.matrix(c(1,-1,rep(0,p-2)))
x <- matrix(rnorm(n*p), nrow = n, ncol = p)
y <- exp(0.5 * x%*%b) + 0.1*rnorm(n)
MSIR <- msir(x, y)
summary(MSIR)
plot(MSIR, type = "evalues")
plot(MSIR, type = "coefficients", which = 1)
plot(MSIR, type = "2Dplot")

Example: 1-dimensional symmetric response curve

n <- 200
p <- 5
b <- as.matrix(c(1,-1,rep(0,p-2)))
x <- matrix(rnorm(n*p), nrow = n, ncol = p)
y <- (0.5 * x%*%b)^2 + 0.1*rnorm(n)
MSIR <- msir(x, y)
summary(MSIR)
plot(MSIR, type = "evalues")
plot(MSIR, type = "coefficients", which = 1)
plot(MSIR, type = "2Dplot")

Example: 2-dimensional response curve

n <- 300
p <- 5
b1 <- c(1, 1, 1, rep(0, p-3))
b2 <- c(1,-1,-1, rep(0, p-3))
b <- cbind(b1,b2)
x <- matrix(rnorm(n*p), nrow = n, ncol = p)
y <- x %*% b1 + (x %*% b1)^3 + 4*(x %*% b2)^2 + rnorm(n)
MSIR <- msir(x, y)
summary(MSIR)
plot(MSIR, type = "evalues")
plot(MSIR, type = "coefficients", which = 1:2)
plot(MSIR, which = 1:2)
plot(MSIR, which = 1, type = "2Dplot", span = 0.7)
plot(MSIR, which = 2, type = "2Dplot", span = 0.7)

To obtain rotating 3D spinplot use:

plot(MSIR, type = "spinplot")
rgl::rglwidget(width=500, height=450)

plot(MSIR, type = "spinplot", span = 0.75)
rgl::rglwidget(width=500, height=450)

General usage of `spinplot()` function

x1 <- rnorm(100)
x2 <- rnorm(100)
y  <- 2*x1 + x2^2 + 0.5*rnorm(100)

spinplot(x1, y, x2)
rgl::rglwidget(width=500, height=450)

spinplot(x1, y, x2, scaling="aaa")
rgl::rglwidget(width=500, height=450)

spinplot(x1, y, x2, rem.lin.trend = "TRUE")
rgl::rglwidget(width=500, height=450)

spinplot(x1, y, x2, fit.smooth = TRUE)
rgl::rglwidget(width=500, height=450)

spinplot(x1, y, x2, fit.ols = TRUE)
rgl::rglwidget(width=500, height=450)

x <- iris[,1:3]
y <- iris[,5]

spinplot(x)
rgl::rglwidget(width=500, height=450)

spinplot(x, markby = y)
rgl::rglwidget(width=500, height=450)

spinplot(x, markby = y, pch = c(0,3,1), 
         col.points = c("lightcyan", "yellow", "lightgreen"), 
         background = "black")
rgl::rglwidget(width=500, height=450)