Ridge GLM pénalisé en utilisant l'augmentation de ligne?

J'ai lu que la régression de crête pourrait être obtenue en ajoutant simplement des lignes de données à la matrice de données d'origine, où chaque ligne est construite en utilisant 0 pour les variables dépendantes et la racine carrée de $k$ ou zéro pour les variables indépendantes. Une ligne supplémentaire est ensuite ajoutée pour chaque variable indépendante.

Je me demandais s'il était possible d'en tirer une preuve pour tous les cas, y compris pour la régression logistique ou d'autres GLM.

La régression de crête minimise .$\sum_{i=1}^n (y_i-x_i^T\beta)^2+\lambda\sum_{j=1}^p\beta_j^2$

(Souvent, une constante est requise, mais pas rétrécie. Dans ce cas, elle est incluse dans le et les prédicteurs - mais si vous ne voulez pas la réduire, vous n'avez pas de ligne correspondante pour la pseudo-observation. Ou si vous ne voulez réduire, vous faire une ligne pour elle. Je vais l' écrire comme si elle ne compte pas dans le p , et non rétréci, comme il est le cas le plus compliqué. l'autre cas est un changement trivial de cette situation . )$\beta$$p$

Nous pouvons écrire le deuxième terme sous la forme de pseudo-observations si nous pouvons écrire chaque "y" et chacun des vecteurs ( p + 1 ) correspondants "x" de telle sorte que$p$$(p+1)$

$(y_{n+j}-x_{n+j}^T\beta)^2=\lambda\beta_j^2\,,\quad j=1,\ldots,p$

Mais par inspection, il suffit de laisser , soit x n + j , j = $y_{n+j}=0$ et laissez tous les autresx n + j , k =0(y comprisx n + j , 0 =0typiquement).$x_{n+j,j}=\sqrt{\lambda}$$x_{n+j,k}=0$$x_{n+j,0}=0$

alors

$(y_{n+j}-[x_{n+j,0}\beta_0+x_{n+j,1}\beta_1+x_{n+j,2}\beta_2+...+x_{n+j,p}\beta_p])^2=\lambda\beta_j^2$

Cela fonctionne pour la régression linéaire. Cela ne fonctionne pas pour la régression logistique, car la régression logistique ordinaire ne minimise pas la somme des résidus au carré.

[La régression de crête n'est pas la seule chose qui peut être faite via de telles astuces de pseudo-observation - elles surviennent dans un certain nombre d'autres contextes]

Il n'est en effet pas difficile de généraliser cette recette aux GLM, car les GLM sont généralement adaptés à l'aide des moindres carrés itérativement repondérés . Par conséquent, à l'intérieur de chaque itération, on peut substituer l'étape des moindres carrés pondérés régulière avec une étape des moindres carrés pondérés pénalisée par une crête pour obtenir un GLM pénalisé par une crête. En fait, en combinaison avec des pénalités de crête adaptatives, cette recette est utilisée pour ajuster les GLM pénalisés L0 (c'est-à-dire le meilleur sous-ensemble, c'est-à-dire les GLM où le nombre total de coefficients non nuls est pénalisé). Ceci a été implémenté par exemple dans le paquet l0ara , voir cet article et celui-ci pour plus de détails.

Il convient également de noter que le moyen le plus rapide de résoudre une régression de crête régulière utilise

lmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  solve(crossprod(X) + diag(lambdas), crossprod(X, y))[, 1]
}

dans le cas où n>=p, ou en utilisant

lmridge_solve_largep = function (X, Y, lambda) (t(X) %*% solve(tcrossprod(X)+lambda*diag(nrow(X)), Y))[,1]

quand p>net pour un modèle sans interception.

C'est plus rapide que d'utiliser la recette d'augmentation de ligne , c'est-à-dire de faire

lmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  qr.solve(rbind(X, diag(sqrt(lambdas))), c(y, rep(0, ncol(X))))
}

S'il vous arrivait d'avoir besoin de contraintes de non négativité sur vos coefficients ajustés, alors vous pouvez simplement faire

library(nnls)

nnlmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x
}

ce qui donne ensuite un résultat légèrement plus précis que

nnlmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=rbind(X,diag(sqrt(lambdas))), b=c(Y,rep(0,ncol(X))))$x 
}

(et à proprement parler seule la solution nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x est alors la bonne).

Je n'ai pas encore compris comment le cas contraint à la non-négativité pourrait être encore optimisé pour le p > ncas - faites-moi savoir si quelqu'un pourrait savoir comment faire cela ... [ lmridge_nnls_largep = function (X, Y, lambda) t(X) %*% nnls(A=tcrossprod(X)+lambda*diag(nrow(X)), b=Y)$xne fonctionne pas]