inst/doc/CorReg.R

In CorReg: Linear Regression Based on Linear Structure Between Variables

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, fig.align = "center", fig.width = 7, fig.height = 5)

## -----------------------------------------------------------------------------
library(CorReg)

# dataset generation
base <- mixture_generator(n = 15, p = 10, ratio = 0.4, tp1 = 1, tp2 = 1, tp3 = 1, positive = 0.5,
                          R2Y = 0.8, R2 = 0.9, scale = TRUE, max_compl = 3, lambda = 1)

X_appr <- base$X_appr # learning sample
Y_appr <- base$Y_appr # response variable for the learning sample
Y_test <- base$Y_test # response variable for the validation sample
X_test <- base$X_test # validation sample
TrueZ <- base$Z # True generative structure (binary adjacency matrix)


## -----------------------------------------------------------------------------
TrueZ
# TrueZ[i,j]==1 means that X[,j] linearly depends on X[,i]


## -----------------------------------------------------------------------------
# density estimation for the MCMC (with Gaussian Mixtures)
density <- density_estimation(X = X_appr, nbclustmax = 10, detailed = TRUE, package = "Rmixmod")
Bic_null_vect <- density$BIC_vect # vector of the BIC found (1 value per covariate)


## -----------------------------------------------------------------------------
# MCMC to find the structure
res <- structureFinder(X = X_appr, verbose = 0, reject = 0, Maxiter = 1500, nbini = 30, 
                       candidates = -1, Bic_null_vect = Bic_null_vect, star = TRUE, p1max = 15, clean = TRUE)
hatZ <- res$Z_opt # found structure (adjacency matrix)
hatBic <- res$bic_opt # associated BIC


## -----------------------------------------------------------------------------
# BIC comparison between true and found structure
bicopt_vect <- BicZ(X = X_appr, Z = hatZ, Bic_null_vect = Bic_null_vect)
bicopt_True <- BicZ(X = X_appr, Z = TrueZ, Bic_null_vect = Bic_null_vect)
sum(bicopt_vect)
sum(bicopt_True)


## -----------------------------------------------------------------------------
# Structure comparison
compZ <- compare_struct(trueZ = TrueZ, Zalgo = hatZ) # qualitative comparison
compZ


## -----------------------------------------------------------------------------
# interpretation of found and true structure ordered by increasing R2
readZ(Z = hatZ, crit = "R2", X = X_appr, output = "all", order = 1) # <NA>line : name of subregressed covariate
readZ(Z = TrueZ, crit = "R2", X = X_appr, output = "all", order = 1) # <NA>line : name of subregressed covariate
  

## -----------------------------------------------------------------------------
# Regression coefficients estimation
select = "NULL" # without variable selection (otherwise, choose "lar" for example to use lasso selection)
resY <- correg(X = X_appr, Y = Y_appr, Z = hatZ, compl = TRUE, expl = TRUE, pred = TRUE,
               select = select, K = 10)


## -----------------------------------------------------------------------------
# MSE computation
MSE_complete <- MSE_loc(Y = Y_test, X = X_test, A = resY$compl$A) # classical model on X
MSE_marginal <- MSE_loc(Y=Y_test, X = X_test, A = resY$expl$A) # reduced model without correlations
MSE_plugin <- MSE_loc(Y = Y_test, X = X_test, A = resY$pred$A) # plug-in model
MSE_true <- MSE_loc(Y = Y_test, X = X_test, A = base$A) # True model
 

## -----------------------------------------------------------------------------
# MSE comparison
MSE <- data.frame(MSE_complete, MSE_marginal, MSE_plugin, MSE_true)
MSE # estimated structure
barplot(as.matrix(MSE), main = "MSE on validation dataset", sub = paste("select=", select), col = "blue")
abline(h = MSE_complete, col = "red")