Cookies, NIR spectrometry

Load the dataset.

32 cookies are tested and the spectra are composed of 700 (correlated) features.

url="https://perso.univ-rennes1.fr/valerie.monbet/doc/cours/Biscuits.csv"
url = "~/Dropbox/ENSEIGNEMENT/RADO/Biscuits.csv"
biscuits=read.csv(url,sep=";")
colnames(biscuits)[1] = "fat"
# Extraction de la colonne fat
fat=biscuits[,1]
# Extraction des variables explicatives
X=biscuits[,-1]

Plot one spectrum as example

# Plot spectra
plot(as.numeric(X[2,]),typ="l")
title("Un exemple de spectre")

In the following plot, color is varying according to the fat percent

fatn=fat/max(fat)
colors= rainbow(length(fat)*2)[rev(1:length(fat))]
i.col = rank(fat)
plot(as.numeric(X[1,]),col=colors[i.col[1]],typ="l",ylab = "Absorbances",ylim = c(0,2.5))
for (i in 2:length(fat)){
  lines(as.numeric(X[i,]),col=colors[i.col[i]],typ="l",ylab = "Absorbances")
}
title("Spectres NIR")

train_test_split <- function (X,y,test_size=.25,random_state=NULL){
  # Extraction of a train and a test datasets.
  #
  # Inputs : X, y : data
  #          test_size : a fraction (<1) of the total set or a number of samples (integer) ; default 0.25
  #          random_state : if equal to an interger, it fixes the random seed ; defaut NULL
  # Outputs : X_train, X_test, y_train, y_test
  #
  n = nrow(X)
  if (test_size>1){test_size=test_size/n}
  if (!is.null(random_state)){set.seed(random_state)}
  itest=sample(1:n,round(n*test_size))
  itrain=setdiff(1:n,itest)
  Xtrain=X[itrain,]
  Xtest=X[itest,]
  ytrain=y[itrain]
  ytest=y[itest]
  return(list(X_train=Xtrain,X_test=Xtest,y_train=ytrain,y_test=ytest,i_test=itest))
}

A regular linear model can not be fit without variable selection.

The lm function does not return an error but fit a model based on the first variables available in the dataset. And the rmse evaluated without cross validation is 0.

D = data.frame(fat,scale(X))
mod = lm(fat~.,data = D)
rmse_mod = sqrt(mean((fat - predict(mod,D))^2))

## Warning in predict.lm(mod, D): prediction from a rank-deficient fit may be
## misleading

print(paste("RMSE (no cross validation): ", round(rmse_mod,2)))

## [1] "RMSE (no cross validation):  0"

Now, using croos validation, we see that the generalization error is pretty high (even larger as the one of the naive prediction!).

y = fat
B = 30 # number of samples for the cross validation
mse_lr = mse_naive = matrix(0,B,1)
for (b in 1:B){
  # sampling for train, test sets
  samples = train_test_split(X,y,test_size=0.25,random_state=1)
  X_train=samples$X_train
  X_test=samples$X_test
  y_train=samples$y_train
  y_test=samples$y_test

  # standardization
  mk=apply(X_train,2,mean)
  sk=apply(X_train,2,sd)
  X_train = scale(X_train,center=mk,scale=sk)
  X_test = scale(X_test,center=mk,scale=sk)
  data_train = data.frame(X_train,y_train)
  data_test = data.frame(X_test,y_test)
  
  # learning linear model
  lr = lm(y_train~.,data=data_train)
  y_pred = predict(lr,data_test) # prediction for the test set
  mse_lr[b] = mean((y_test-y_pred)^2) # compute MSE

  y_pred_naive = mean(y_train)
  mse_naive = mean((y_test-y_pred_naive)^2)
}

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(lr, data_test): prediction from a rank-deficient fit may
## be misleading

print(paste("RMSE regression  (no CV): ", round(rmse_mod,2)))

## [1] "RMSE regression  (no CV):  0"

print(paste("RMSE regression (CV): ", round(sqrt(mean(mse_lr)),2)))

## [1] "RMSE regression (CV):  4.42"

print(paste("RMSE naive prediction: ", round(sqrt(mean(mse_naive)),2)))

## [1] "RMSE naive prediction:  1.76"

Feature selection using the correlation with the variable to be predicted

Less us keep the 10 "best" variables (the number 10 is chosen arbitrary).

The scatter plot shows that they are highly correlated to eah other.

p = dim(X)[2]
rho = NULL
for (j in 1:p){
  rho[j] = cor(X[,j],fat)
}

rho_10 = sort(rho)[p-10]
keep = which(rho>=rho_10) # 10 best variables

pairs(X[,keep],pch = 20)

p = dim(X)[2]
rho = NULL
for (j in 1:p){
  rho[j] = cor(X[,j],fat)
}

rho_10 = sort(rho)[p-10]
keep = which(rho>=rho_10) 

fatn=fat/max(fat)
colors= rainbow(length(fat)*2)[rev(1:length(fat))]
i.col = rank(fat)
plot(as.numeric(X[1,]),col=colors[i.col[1]],typ="l",ylab = "Absorbances",ylim = c(0,2.5),
     main = "Correlation selection")
matlines(t(matrix(c(keep,keep),length(keep),2)),matrix(c(0,2.5),2,length(keep)),col="black")
for (i in 2:length(fat)){
  lines(as.numeric(X[i,]),col=colors[i.col[i]],typ="l",ylab = "Absorbances")
}

Remark that the feature selection as to be done on the train dataset.

B = 30 # number of samples for the cross validation
mse_lr_cor = mse_naive = matrix(0,B,1)
for (b in 1:B){

  # sampling for train, test sets
  samples = train_test_split(X,y,test_size=0.25,random_state=1)
  X_train=samples$X_train
  X_test=samples$X_test
  y_train=samples$y_train
  y_test=samples$y_test

  rho = NULL
  for (j in 1:p){
    rho[j] = cor(X_train[,j],y_train)
  }
  rho_10 = sort(rho)[p-10]
  keep = which(rho>=rho_10) # 10 best variables

  # standardization
  mk=apply(X_train,2,mean)
  sk=apply(X_train,2,sd)
  X_train = scale(X_train,center=mk,scale=sk)
  X_test = scale(X_test,center=mk,scale=sk)
  data_train = data.frame(y_train,X_train[,keep])
  data_test = data.frame(y_test,X_test[,keep])
  
  # learning linear model
  colnames(data_train)[1] = "fat"
  colnames(data_test)[1] = "fat"
  lr = lm(fat~.,data=data_train)
  y_pred = predict(lr,data_test) # prediction for the test set
  mse_lr_cor[b] = mean((y_test-y_pred)^2) # compute MSE

  y_pred_naive = mean(y_train)
  mse_naive = mean((y_test-y_pred_naive)^2)
}
print(paste("RMSE linear regression  (no CV): ", round(rmse_mod,2)))

## [1] "RMSE linear regression  (no CV):  0"

print(paste("RMSE linear regression  with correlation selection (CV): ", round(sqrt(mean(mse_lr_cor)),2)))

## [1] "RMSE linear regression  with correlation selection (CV):  1.17"

print(paste("Error of naive prediction: ", round(sqrt(mean(mse_naive)),2)))

## [1] "Error of naive prediction:  1.76"

Another trick to reduce the dimension is to run a PCA and use the firsts component as input in the linear model. The numer of component to be included in the model is chosen arbitrarely equal to 7. This number may affect the quality of the model.

Here again, the PCA has to be run on the train dataset.

n = dim(X)[1]
B = 30 # number of samples for the cross validation
mse_lr = mse_pca = NULL
for (b in 1:B){
  samples = train_test_split(X,fat,test_size=.25,random_state=1)
  X_train=samples$X_train
  X_test=samples$X_test
  y_train=samples$y_train
  y_test=samples$y_test
  
  mk=apply(X_train,2,mean)
  sk=apply(X_train,2,sd)
  X_train = scale(X_train,center=mk,scale=sk)
  X_test = scale(X_test,center=mk,scale=sk)
  
  data_train = data.frame(y_train,X_train)
  data_test = data.frame(y_test,X_test)
  colnames(data_train)[1] = colnames(data_test)[1] = "fat"
  mod = lm(fat~.,data=data_train)
  y_pred = predict(mod,data_test)  
  mse_lr[b] = mean((y_test-y_pred)^2)
  
  pca = PCA(X,ncp=7,graph=FALSE,ind.sup=samples$i_test)  
  data_train_pca = data.frame(y_train,pca$ind$coord)
  data_test_pca = data.frame(y_test,pca$ind.sup$coord)  
  colnames(data_train_pca)[1] = colnames(data_test_pca)[1] = "fat"
  lr_pca = lm(fat~.,data=data_train_pca) 
  y_pred_pca = predict(lr_pca,data_test_pca)
  mse_pca[b] = mean((y_test-y_pred_pca)^2)
}

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

## Warning in predict.lm(mod, data_test): prediction from a rank-deficient fit may
## be misleading

print(paste("RMSE of linear regression:", round(sqrt(mean(mse_lr)),2)))

## [1] "RMSE of linear regression: 4.42"

print(paste("RMSE of linear regression on pca components:", round(sqrt(mean(mse_pca)),2)))

## [1] "RMSE of linear regression on pca components: 0.18"

Fit Ridge and Lasso regression

alpha_values = c(1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1,5) # valeurs lambda


B = 30 # number of samples for the cross validation
mse_ridge = mse_lasso = matrix(0,B,length(alpha_values))
for (b in 1:B){
  samples = train_test_split(X,fat,test_size=.25,random_state=1)
  X_train=samples$X_train
  X_test=samples$X_test
  y_train=samples$y_train
  y_test=samples$y_test  
  
  mk=apply(X_train,2,mean)
  sk=apply(X_train,2,sd)
  X_train = scale(X_train,center=mk,scale=sk)
  X_test = scale(X_test,center=mk,scale=sk)
  
  for (i in 1:length(alpha_values)){
     ridge = glmnet(as.matrix(X_train),y_train,alpha=0,lambda=alpha_values[i])
     y_pred = predict(ridge,as.matrix(X_test))
     mse_ridge[b,i] = mean((y_pred-y_test)^2)
     lasso = glmnet(as.matrix(X_train),y_train,alpha=1,lambda=alpha_values[i])
     y_pred = predict(lasso,as.matrix(X_test))
     mse_lasso[b,i] = mean((y_pred-y_test)^2)
  }
}
     
print(paste("RMSE of Ridge regression (CV):", round(sqrt(min(apply(mse_ridge,2,mean))),2)))

## [1] "RMSE of Ridge regression (CV): 0.24"

print(paste("RMSE of Lasso regression  (CV): ", round(sqrt(min(apply(mse_lasso,2,mean))),2)))

## [1] "RMSE of Lasso regression  (CV):  0.25"

print(paste("RMSE of linear regression on pca components  (CV):", round(sqrt(mean(mse_pca)),2)))

## [1] "RMSE of linear regression on pca components  (CV): 0.18"

print(paste("RMSE linear regression  with correlation selection (CV): ", round(sqrt(mean(mse_lr_cor)),2)))

## [1] "RMSE linear regression  with correlation selection (CV):  1.17"

Interpretability of the models

It is usually interesting to go back to the feature selection. For spectrometry data, it is natural to plot the spectra and highlight the locations of the selected variables.

For that, let us use all the data.

p = dim(X)[2]
rho = NULL
for (j in 1:p){
  rho[j] = cor(X[,j],fat)
}

rho_10 = sort(rho)[p-10]
keep = which(rho>=rho_10) 

fatn=fat/max(fat)
colors= rainbow(length(fat)*2)[rev(1:length(fat))]
i.col = rank(fat)
plot(as.numeric(X[1,]),col=colors[i.col[1]],typ="l",ylab = "Absorbances",ylim = c(0,2.5),
     main = "Correlation selection")
matlines(t(matrix(c(keep,keep),length(keep),2)),matrix(c(0,2.5),2,length(keep)),col="black")
for (i in 2:length(fat)){
  lines(as.numeric(X[i,]),col=colors[i.col[i]],typ="l",ylab = "Absorbances")
}

Selection driven by the PCA.

pca = PCA(scale(X),ncp=7,graph=FALSE) 
res = dimdesc(pca)
keep =  which(res[[1]][[1]][,2]<.05 | res[[2]][[1]][,2]<.05 )

## Warning in res[[1]][[1]][, 2] < 0.05 | res[[2]][[1]][, 2] < 0.05: la taille d'un
## objet plus long n'est pas multiple de la taille d'un objet plus court

fatn=fat/max(fat)
colors= rainbow(length(fat)*2)[rev(1:length(fat))]
i.col = rank(fat)
plot(as.numeric(X[1,]),col=colors[i.col[1]],typ="l",ylab = "Absorbances",ylim = c(0,2.5),
     main = "PCA selection")
matlines(t(matrix(c(keep,keep),length(keep),2)),matrix(c(0,2.5),2,length(keep)),col="black")
for (i in 2:length(fat)){
  lines(as.numeric(X[i,]),col=colors[i.col[i]],typ="l",ylab = "Absorbances")
}

Selection driven by the PCA.

pca = PCA(scale(X),ncp=7,graph=FALSE) 
res = dimdesc(pca)

par(mfrow=c(2,1))
plot(res[[1]][[1]][,1],pch=20,ylab= "Cor(PC1,X)")
plot(res[[2]][[1]][,1],pch=20,ylab= "Cor(PC2,X)")

Selection of the Lasso regression

cv = cv.glmnet(scale(X),fat)
mod = glmnet(scale(X),fat,lambda=cv$lambda.1se)
keep =  mod$beta@i

fatn=fat/max(fat)
colors= rainbow(length(fat)*2)[rev(1:length(fat))]
i.col = rank(fat)
plot(as.numeric(X[1,]),col=colors[i.col[1]],typ="l",ylab = "Absorbances",ylim = c(0,2.5),
     main = "Lasso selection")
matlines(t(matrix(c(keep,keep),length(keep),2)),matrix(c(0,2.5),2,length(keep)),col="black")
for (i in 2:length(fat)){
  lines(as.numeric(X[i,]),col=colors[i.col[i]],typ="l",ylab = "Absorbances")
}

Variables selection - Cookies

Cookies, NIR spectrometry

Feature selection using the correlation with the variable to be predicted

Fit Ridge and Lasso regression

Interpretability of the models