Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #--- Motecarlo ---#
- library(glmnet)
- library(mvtnorm)
- library(Matrix)
- # Functions
- correlation_matrix <- function(p) {
- a = runif(1,0.1,0.9)
- b = runif(1,0.1,0.9)
- corr_min = min(a,b)
- corr_max = max(a,b)
- Sigma <- diag(1, p)
- Sigma[lower.tri(Sigma)] <- runif((p^2-p)/2, corr_min, corr_max)
- Sigma[upper.tri(Sigma)] <- t(Sigma)[upper.tri(Sigma)]
- Sigma
- Sigma <- nearPD(Sigma, corr = TRUE, keepDiag = TRUE) # trova la matrice definita positiva più vicina alla matrice di input utilizzando l’algoritmo di Higham
- Sigma$mat
- return(list(matrix = Sigma$mat, mean_correlation = mean(Sigma$mat[upper.tri(Sigma$mat)]))) # restituisce matrice e correlazione media delle variabili
- }
- # Data generating process
- nsim <- 500 # numero simulazioni
- n <- 200 # numero osservazioni
- p <- 100 # numero varaibili indipendenti
- significative <- 20 # numero variabili indipendenti significative
- beta_vector <- runif(p,-5,5)
- # beta_value <- 8 # valore del coeff beta
- # beta_vector <- c(rep(beta_value, significative), rep(0,p-significative)) # vettore dei beta
- sigma_X <- (1/n)^.5 # deviazione standard della distribuzione normale utilizzata per generare le variabili indipendenti
- sigma_epsilon <- 1 # deviazione standard della distribuzione normale utilizzata per generare l'errore
- df <- data.frame(mean_sigma_inv = numeric(0), mean_correlation = numeric(0), determinant = numeric(0), mse_ols = numeric(0), mse_ridge = numeric(0), optimal_lambda = numeric(0))
- for (i in 1:nsim){
- set.seed(i+1)
- # Sigma <- clusterGeneration::rcorrmatrix(p, alphad=1) # matrice di varianza e covarianza che rappresenta appunto la correlazione tra le p variabili indipendenti
- res <- correlation_matrix(p)
- Sigma <- as.matrix(res[1]$matrix)
- mean_corr <- res[2]
- Sigma
- mean_corr
- # Multicollinearità
- determinant <- det(Sigma) # http://www.stat.unipg.it/~bart/metodi/lezione14.pdf Quando c’è una forte correlazione tra due o più variabili esplicative in un modello di regressione, il determinante della matrice di varianza e covarianza può diventare molto piccolo
- Sigma_inv <- solve(Sigma) # Matrice inversa
- mean_sigma_inv <- mean(Sigma_inv[upper.tri(Sigma_inv)])
- mean_sigma_inv
- epsilon <- rnorm(n,0,sigma_epsilon) # genero l'errore da una normale N(0,1)
- X <- mvtnorm::rmvnorm(n,sigma = Sigma) # genero la matrice di dati nxp, ogni riga rappresenta un'osservazione, ogni colonna una variabile indipendente
- Y <- X%*%beta_vector + epsilon # genero la variabile dipendente come combinazione lineare delle X e del vettore dei beta con l'aggiunta di un errore
- # OLS Estimator
- fit.ols <- lm(Y ~ X)
- summary(fit.ols)
- stime_ols <- fit.ols$coefficients
- stime_ols
- # Scelgo il lambda ottimale mediante CV
- cv.ridge <- cv.glmnet(x=X,y=Y,family="gaussian", alpha=0)
- # cv.ridge <- cv.glmnet(x[train,],y[train], nfolds = 10, type.measure = "mse",alpha=0)
- # plot(cv.ridge)
- optimal_lambda <- cv.ridge$lambda.min
- optimal_lambda
- # Ridge Estimator
- # Prova senza family
- fit.ridge <- glmnet(x=X,y=Y,family="gaussian", alpha=0, lambda = optimal_lambda)
- stime_ridge <- fit.ridge$beta
- # MSE
- mse_ols <- mean((stime_ols[-1] - beta_vector)^2) # il primo valore in stime_ols è l'intercetta (?)
- mse_ridge <- mean((stime_ridge - beta_vector)^2)
- new_row <- data.frame(mean_sigma_inv=mean_sigma_inv,mean_correlation=mean_corr, determinant = determinant, mse_ols = mse_ols, mse_ridge = mse_ridge, optimal_lambda=optimal_lambda)
- df <- rbind(df, new_row)
- }
- library(ggplot2)
- plot(df$mean_correlation, df$mse_ols)
- plot(df$mean_correlation, df$mse_ridge)
- ggplot(df, aes(x = df$determinant)) +
- geom_point(aes(y = df$determinant), color = "red", alpha=0.7) +
- scale_y_log10()
- ggplot(df, aes(x = df$mean_correlation)) +
- geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
- geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
- scale_y_log10()
- ggplot(df, aes(x = df$mean_sigma_inv)) +
- geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
- geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
- scale_y_log10()
- ggplot(df, aes(x = df$mean_sigma_inv)) +
- geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
- geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3)
- ggplot(df, aes(x = df$determinant)) +
- geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
- geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
- geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3)
- ggplot(df, aes(x = df$mean_correlation)) +
- geom_point(aes(y = df$optimal_lambda), color = "red", alpha=0.7) +
- geom_smooth(aes(y = df$optimal_lambda), color = "red",alpha=0.3) +
- scale_y_log10()
- ggplot(df, aes(x = df$mean_sigma_inv)) +
- geom_point(aes(y = df$optimal_lambda), color = "red", alpha=0.7) +
- geom_smooth(aes(y = df$optimal_lambda), color = "red",alpha=0.3) +
- scale_y_log10()
- # La multicollinearità va guardata sull'inversa
- # Posso calcolare il determinante della matrie inversa per capire il livello di correlazione
- # --- Scenario 2 --- #
- library(glmnet)
- library(mvtnorm)
- library(Matrix) # Per trovare la semidefinita positiva piu vicina
- # Functions
- correlation_matrix <- function(p) {
- a = runif(1,0.1,0.9)
- b = runif(1,0.1,0.9)
- corr_min = min(a,b)
- corr_max = max(a,b)
- Sigma <- diag(1, p)
- Sigma[lower.tri(Sigma)] <- runif((p^2-p)/2, corr_min, corr_max)
- Sigma[upper.tri(Sigma)] <- t(Sigma)[upper.tri(Sigma)]
- Sigma
- Sigma <- nearPD(Sigma, corr = TRUE, keepDiag = TRUE) # trova la matrice definita positiva più vicina alla matrice di input utilizzando l’algoritmo di Higham
- Sigma$mat
- return(Sigma$mat) # restituisce matrice e correlazione media delle variabili
- }
- # Data generating process
- nsim <- 500 # numero simulazioni
- n <- 100 # numero osservazioni
- p <- 200 # numero varaibili indipendenti
- significative <- 20 # numero variabili indipendenti significative
- beta_vector <- runif(p,-5,5)
- # beta_value <- 8 # valore del coeff beta
- # beta_vector <- c(rep(beta_value, significative), rep(0,p-significative)) # vettore dei beta
- sigma_X <- (1/n)^.5 # deviazione standard della distribuzione normale utilizzata per generare le variabili indipendenti
- sigma_epsilon <- 1 # deviazione standard della distribuzione normale utilizzata per generare l'errore
- df <- data.frame(mse_ridge=numeric(0),mse_lasso=numeric(0))
- for (i in 1:nsim){
- set.seed(i+1)
- # Sigma <- clusterGeneration::rcorrmatrix(p, alphad=1) # matrice di varianza e covarianza che rappresenta appunto la correlazione tra le p variabili indipendenti
- Sigma <- correlation_matrix(p)
- Sigma <- as.matrix(Sigma)
- Sigma
- epsilon <- rnorm(n,0,sigma_epsilon) # genero l'errore da una normale N(0,1)
- X <- mvtnorm::rmvnorm(n,sigma = Sigma) # genero la matrice di dati nxp, ogni riga rappresenta un'osservazione, ogni colonna una variabile indipendente
- Y <- X%*%beta_vector + epsilon # genero la variabile dipendente come combinazione lineare delle X e del vettore dei beta con l'aggiunta di un errore
- # Ridge Estimator
- fit.ridge <- glmnet(x=X,y=Y, alpha=0)
- stime_ridge <- fit.ridge$beta
- # Lasso Estimator
- fit.lasso <- glmnet(x=X,y=Y, alpha=1)
- stime_lasso <- fit.lasso$beta
- stime_lasso
- # MSE
- mse_ridge <- mean((stime_ridge - beta_vector)^2)
- mse_lasso <- mean((stime_lasso - beta_vector)^2)
- new_row <- data.frame(mse_ridge=mse_ridge,mse_lasso=mse_lasso)
- df <- rbind(df, new_row)
- }
- mean_mse_ridge <- mean(df$mse_ridge)
- mean_mse_lasso <- mean(df$mse_lasso)
- mean_mse_ridge
- mean_mse_lasso
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement