library(glasso) library(mvtnorm) library(glmnet) ## EXERCICE 1 ---------------------------------------------------------------------------------------------- # a) generate_gamma = function(p){ gamma = matrix(rep(0,p^2),ncol=p) for (i in 1:p){ for (j in 1:p){ gamma[i,j] = min(i,j) } } return(gamma) } p = 10 gamma = generate_gamma(p) # c) K = solve(gamma) # matrice tridiagonale : # 1 si i=j=p # 2 si i=j!=p # -1 si abs(j-i)=1 # 0 sinon K%*%gamma # identité # d) # 1--2--3--4--5 etc jusqu'a p # e) n = 5 p = 10 gamma = generate_gamma(p) A = rmvnorm(n,rep(0,p),gamma) A ## 2 # f) g1 = glasso(var(A),rho=0.1) g2 = glasso(var(A),rho=1) g3 = glasso(var(A),rho=5) g4 = glasso(var(A),rho=10) # les Zi sont independants 2 a 2 g5 = glasso(var(A),rho=100) # les Zi sont independants 2 a 2 # Pour rho grand, on ne retrouve pas le modèle graphique de Z a cause de l'independance entre Zi et Zj (i!=j) # Pour rho petit, on ne retrouve pas le modèle graphique de Z car il y a beaucoup de dependances (par ex, rho=0.1, toutes les variables sont dependantes) # g) gpath = glassopath(var(A)) gpath # h) n = 10 p = 100 gamma = generate_gamma(p) A = rmvnorm(n,rep(0,p),gamma) A glasso(var(A),rho=5) # glassopath(var(A)) # tres long # On peut aussi passer n à 50 pour avoir une meilleure convergence: n = 50 p = 100 gamma = generate_gamma(p) A = rmvnorm(n,rep(0,p),gamma) A glasso(var(A),rho=10) ## 3 # i) n = 50 p = 10 gamma = generate_gamma(p) A = rmvnorm(n,rep(0,p),gamma) A i = 4 lasso4 = glmnet(A[,-i],A[,i],alpha=1) coef(lasso4,FALSE) # ici on a des coeffictients non nuls pour les variables 2,3,7 et 8 # Les interpretations changent en fonction de la compilation de A lambdamin4 = cv.glmnet(A[,-i],A[,i])$lambda.min lambdamin4 for (i in 1:p){ lambdamin = cv.glmnet(A[,-i],A[,i])$lambda.min print(lambdamin) } # WARNINGS car le jeu de données impose de tailles de folds trop petites # lambdamin dépend de i plot(lasso4) abline(v=c(lambdamin4),lty=4,col="orange") # On ne retrouve pas le modele graphique de Z, on le voit dans les coefs du lasso on devrait avoir des coefficients # nuls partout sauf pour 3,4 : ce n'est pas le cas # On peut aussi passer n tres grand pour avoir une meilleure convergence: n = 1000 p = 10 gamma = generate_gamma(p) A = rmvnorm(n,rep(0,p),gamma) i = 4 lasso4 = glmnet(A[,-i],A[,i],alpha=1) coef(lasso4,FALSE) # ici on a des coeffictients tres proches de 0 pour les variables autres que 3 et 4 lambdamin4 = cv.glmnet(A[,-i],A[,i])$lambda.min lambdamin4 glasso(var(A),rho=10) plot(lasso4) abline(v=c(lambdamin4),lty=4,col="orange") # j) glasso(var(A),rho=5,approx=TRUE) # Cette methode renvoie une matrice de variance/covariance nulle ! # k) # rho est un argument obligatoire pour compiler glasso # on peut se passer de rho avec la fonction glassopath mais elle est tres longue # Si on ne donne pas de valeur lambda a la fonction lasso, elle le fixe elle meme en fonction du jeu de donnees # On peut trouver le meilleur lambda a lui donner en parametre avec une validation croisee # Le parametre a donner est differents entre lasso et glasso mais joue un grand role dans l'algorithme ## EXERCICE 2 ---------------------------------------------------------------------------------------------- # a) library(boot) data("frets") head(frets) g1 = glasso(var(frets),rho=0.1) g2 = glasso(var(frets),rho=1) g3 = glasso(var(frets),rho=5) g4 = glasso(var(frets),rho=10) g5 = glasso(var(frets),rho=100) # rho in {0.1, 1, 5, 10} : Toutes les variables sont correlees # rho=100 : Aucune correlation entre les variables g = glassopath(var(frets)) cov2cor(g1$wi) # matrice des correlations cov2cor(g2$wi) cov2cor(g3$wi) cov2cor(g4$wi) cov2cor(g5$wi) L = list(g4$wi,g5$wi,g$wi[,,10]) L[[1]] L[[2]] L[[3]] for (i in 1:3){ a = L[[i]] a[abs(a)<0.01] = 0 L[[i]] = a } L[[1]] # matrice diagonale L[[2]] # matrice nulle L[[3]] # matrice nulle # c) # Pour rho in {0.1, 1, 5, 10} : tous les sommets sont relies par des aretes # Pour rho = 100 : le graphe ne contient aucune arete # d) correlations = cov2cor(glasso(var(frets),rho=0)$wi) correlations # e) g1 = glasso(var(frets),rho=0.1,approx=TRUE) g2 = glasso(var(frets),rho=1,approx=TRUE) g3 = glasso(var(frets),rho=5,approx=TRUE) g4 = glasso(var(frets),rho=10,approx=TRUE) g5 = glasso(var(frets),rho=100,approx=TRUE) # PRENDRE UN RHO INTERMEDIAIRE : # PAS TROP GRAND POUR NE PAS PERDRE TOUTE L'INFORMATION # PAS TROP PETIT POUR NE PAS CONSERVER TOUTE LES COVARIANCES ET ELIMINER DES ARETES # Rho est fondamental dans l'algorithme, le resultat depend enormement de son choix # Fonction qui fait varier rho pour obtenir un bon resultat # X tableau de donnees # del_rate la proportion de correlations a eliminer mein_bull = function(X,del_rate=0.5){ v = var(X) n = nrow(v) N = n*(n-1)/2 R = 1 turn_ok = TRUE while (turn_ok){ g = glasso(v,rho=R,approx=TRUE)$wi u = c() for (i in 1:(n-1)){ for (j in (i+1):n){ u = c(u,g[i,j]) } } s = sum(u==0) if (s>N*del_rate){ print(R) turn_ok = FALSE } else { R = R+1 } } return(g) } mein_bull(frets) mein_bull(frets,0.25) mein_bull(frets,0.75) # On retient principalement la correlation entre V1 et V3 soit : Head length of first son et Head length of second son # On retient de plus la correlation entre V1-V2 et V3-V4 soit : Head length et Head breadth pour chaque enfant # la reponse attendue etait : # L1 ___________ L2 soit V1 ________ V2 # | | | | # B1 ___________ B2 V3 ________ V4 # On ne retrouve donc pas le lien entre V2 et V3 avec notre interpretation