test-Kernel_Ridge_MM.R
In KRMM: Kernel Ridge Mixed Model

### SOURCE LIBRARY KRMM
library(KRMM)


### SIMULATE DATA 
set.seed(123)
p=200
N=100

beta=rnorm(p, mean=0, sd=1.0)
X=matrix(runif(p*N, min=0, max=1), ncol=p, byrow=TRUE)  #X: covariates (i.e. predictors)

f=X%*%beta                                              #f: data generating process (i.e. DGP)
E=rnorm(N, mean=0, sd=0.5)

Y=f+E                                                  #Y: response data

hist(f)
hist(beta)
Nb_train=floor((2/3)*N)


###	CREATE TRAINING AND TARGET SETS FOR RESPONSE AND PREDICTOR VARIABLES

Index_train=sample(1:N, size=Nb_train, replace=FALSE)

### Covariates (i.e. predictors) for training and target sets

Predictors_train=X[Index_train, ]
Response_train=Y[Index_train]

Predictors_target=X[-Index_train, ]
Response_target=Y[-Index_train]
True_value_target=f[-Index_train]    #True value (generated by DGP) we want to predict

length(Response_train)
length(Response_target)
N


### PREDICTION WITH KERNEL RIDGE REGRESSION SOLVED WITHIN THE MIXED MODEL FRAMEWORK  


Linear_KRR_model_train = Kernel_Ridge_MM( Y_train=Response_train, Matrix_covariates_train=Predictors_train, method="RR-BLUP")

# Measure goodness of fit of non-optimized model on training set

f_hat_train_Linear_KRR = Predict_kernel_Ridge_MM( Linear_KRR_model_train, Matrix_covariates_target=Predictors_train )

dev.new()
plot(f_hat_train_Linear_KRR,Response_train, main="Linear RKHS regression with default rate of decay (not optimized)")
cor(f_hat_train_Linear_KRR,Response_train)

# Predict new entries for target set and measure prediction error

f_hat_target_Linear_KRR = Predict_kernel_Ridge_MM( Linear_KRR_model_train, Matrix_covariates_target=Predictors_target )

mean((f_hat_target_Linear_KRR-True_value_target)^2)


### Unoptimized Gaussian Kernel 

Gaussian_KRR_model_train = Kernel_Ridge_MM( Y_train=Response_train, Matrix_covariates_train=Predictors_train, method="RKHS")

# Measure goodness of fit of non-optimized model on training set

f_hat_train_Gaussian_KRR = Predict_kernel_Ridge_MM( Gaussian_KRR_model_train, Matrix_covariates_target=Predictors_train )

dev.new()
plot(f_hat_train_Gaussian_KRR,Response_train, main="Gaussian RKHS regression with default rate of decay (not optimized)")
cor(f_hat_train_Gaussian_KRR,Response_train)

# Predict new entries for target set and measure prediction error

f_hat_target_Gaussian_KRR = Predict_kernel_Ridge_MM( Gaussian_KRR_model_train, Matrix_covariates_target=Predictors_target )

mean((f_hat_target_Gaussian_KRR-True_value_target)^2)



\dontrun{


### Tuned Gaussian Kernel 

#Tune the model with the Tune_kernel_Ridge_MM function

Tuned_Gaussian_KRR_train = Tune_kernel_Ridge_MM( Y_train=Response_train, Matrix_covariates_train=Predictors_train,
                                                 method='RKHS', rate_decay_grid=seq(1,20,length.out=10), nb_folds=5, loss='mse' )

Tuned_Gaussian_KRR_model_train = Tuned_Gaussian_KRR_train$tuned_model
Tuned_Gaussian_KRR_train$optimal_h
Tuned_Gaussian_KRR_train$rate_decay_grid
Tuned_Gaussian_KRR_train$expected_loss_grid

dev.new()
plot(Tuned_Gaussian_KRR_train$rate_decay_grid, Tuned_Gaussian_KRR_train$expected_loss_grid, type="l", main="Tuning the rate of decay (for Gaussian kernel) with K-folds cross-validation")


# Measure goodness of fit of tuned model on training set 

f_hat_train_tuned_Gaussian_KRR = Predict_kernel_Ridge_MM( Tuned_Gaussian_KRR_model_train, Matrix_covariates_target=Predictors_train )

dev.new()
plot(f_hat_train_tuned_Gaussian_KRR,Response_train, main="Tuned Gaussian RKHS regression")
cor(f_hat_train_tuned_Gaussian_KRR,Response_train)


# Predict new entries for target set and measure prediction error

f_hat_target_tuned_Gaussian_KRR = Predict_kernel_Ridge_MM( Tuned_Gaussian_KRR_model_train, Matrix_covariates_target=Predictors_target )

mean((f_hat_target_tuned_Gaussian_KRR-True_value_target)^2)


}