În timpul liber, m-am gândit puțin la metodele de control sintetice ca abordare a evaluării și lucrez la o postare pe blog. Dar o problemă secundară care a apărut a fost suficient de interesantă pentru a fi tratată separat.
Din punct de vedere tehnic, este o chestiune de ajustare a unei regresii în care coeficienții sunt constrânși să fie nenegativi și să însumeze unul; cu alte cuvinte, coeficienții sunt a simplex. Motivul pentru care apare acest lucru este că, în metodele de control sintetic, o astfel de regresie este utilizată pentru a determina ponderile de utilizat în construirea unei medii ponderate a mai multor unități (țări, firme, oameni) care este cât mai comparabilă cu unitatea care a primit. intervenția în curs de evaluare.
Argumentul pentru adunarea ponderilor la unu este evitarea extrapolării din intervalul datelor reale. Nu mi se pare un argument convingător, dar îl voi lăsa pentru o postare ulterioară, dacă reușesc. Deocamdată, sunt interesat de întrebarea tehnică a modului de a obține un astfel de set de greutăți diferite de zero însumând până la una; sau abstragând-o departe de aplicație, cum să se potrivească o regresie în care coeficienții sunt constrânși să fie un simplex.
Desigur, există mai multe moduri în care ați putea face acest lucru. A existat o întrebare și un răspuns bun la acest lucru pe validarea încrucișată încă din 2012. Două dintre cele patru metode prin care voi trece mai jos provin aproape direct din răspunsurile sau comentariile de acolo. Celelalte două sunt invenția mea proprie (nu pretind că sunt primul, doar că nu i-am văzut pe alții făcând asta).
Mai întâi, să începem prin a simula unele date. O să fac o matrice numită X cu cinci coloane, toate din distribuții diferite. Apoi creez un vector de coeficienți adevărați care vor lega y Sunt pe cale să fac asta X. Pentru a face exercițiul interesant și realist, coeficienții adevărați vor implica un număr negativ și nu se vor adăuga exact la 1. Unul dintre ei va fi zero.
library(glmnet)
library(MASS) # for eqscplot
library(pracma) # for lsqlincon
library(rstan)
options(mc.cores = min(8, parallel::detectCores())
par(bty = "l", family = "Roboto")
n <- 1000
X <- cbind(
x1 = rnorm(n , 10, 2),
x2 = rnorm(n, 7, 1),
x3 = rnorm(n, 15, 3),
x4 = runif(n, 5, 20),
x5 = rnorm(n, 8, 1.5)
)
# note true data generating process coefficients add up
# to more than 1 and have a negative, to make it interesting
true_coefs <- c(0.2, 0.6, -0.1, 0, 0.4)
y <- X %*% true_coefs + rnorm(n, 0, 2)
Acum, voi potrivi diverse regresii (fără interceptare) la aceste date și voi extrage coeficienții și îi voi compara cu adevăratul proces de generare a datelor. Mai întâi, să începem cu cele mai mici pătrate obișnuite, fără constrângeri
#-------------method 1 - unconstrained regression------
# This is the best way to recover the true DGP
mod1 <- lm(y ~ X - 1)
coefs1 <- coef(mod1)
# plotting function we will use for each model
f <- function(newx, main = ""){
eqscplot(true_coefs, newx, bty = "l",
xlab = "Coefficients from original data generating process",
ylab = "Coefficients from model",
main = main)
abline(0, 1, col = "grey")
}
f(coefs1, "OLS")
Această abordare este foarte bună pentru a se apropia de adevăratul proces de generare a datelor, așa cum ne-am aștepta. Dar, desigur, ca și procesul adevărat, are un coeficient negativ pentru x3 și însumează mai mult de 1:
> round(coefs1, 2) Xx1 Xx2 Xx3 Xx4 Xx5 0.20 0.61 -0.11 -0.02 0.45 > sum(coefs1) (1) 1.122484
Acesta este graficul pe care îl obținem, comparând coeficienții modelului cu cei din procesul real de generare a datelor:
OK, totul foarte bine, dar nu face niciun efort (în afară de doar speranța că funcționează) pentru a constrânge coeficienții să se adună la unul. Așadar, iată prima mea încercare reală, care se bazează pe o sugestie pe care whuber a făcut-o în comentariile la postarea cu validare încrucișată. El subliniază cu o algebră de bază că, dacă scădeți una dintre coloanele X din celelalte și din y și apoi potriviți o regresie cu MCO la acele date transformate, obțineți coeficienții corecti pe toate coloanele de X pe care le-ați lăsat în , și îl poate calcula pe cel rămas scăzându-le pe toate din 1.
#-----------method 2 - transform-----------
# This way will get you coefficients that add up to 1,
# but might not be in (0, 1). This is suggested by whuber in
# his comment on Elvis' solution at:
# https://stats.stackexchange.com/questions/21565/how-do-i-fit-a-constrained-regression-in-r-so-that-coefficients-total-1
X2 <- X(, 1:4) - X(, 5)
y2 <- y - X(, 5)
mod2 <- lm(y2 ~ X2 - 1)
coefs2 <- c(coef(mod2), 1 - sum(coef(mod2)))
f(coefs2, "Transformation") # draw plot
De data aceasta, coeficienții se adună exact la unu după cum este necesar, dar, din păcate, există încă un negativ acolo
> round(coefs2, 2) X2x1 X2x2 X2x3 X2x4 0.24 0.40 -0.04 0.00 0.40 > sum(coefs2) (1) 1
Aceasta nu este o știre – a fost identificată ca o problemă în discuția despre Cross-Validated, dar cu datele simulate folosite acolo nu s-a manifestat. De fapt, acesta este unul dintre motivele pentru care datele mele simulate sunt puțin mai murdare, cu coeficienți originali care nu se adaugă la unul și care au o valoare negativă autentică printre ei.
Cam atât pentru modelul 2.
Modelul 3 este unul pe care l-am inventat eu și folosește faptul că glmnet nu vă permite doar să vă regularizați coeficienții (strângeți-i spre zero, pentru a ajuta la gestionarea multicoliniarității, supraajustării și a problemelor conexe), ci vă permite să specificați limite superioare și inferioare pentru coeficienți. Așadar, am putea obține primul nostru set complet de coeficienți de regresie (nenegativi, adăugând la unul) utilizând abordarea transformării datelor în model2, dar potrivind modelul cu glmnet în schimb. Ca aceasta:
#----------method 3 - transform, regularise and constrain-----------
# This method is guaranteed to get you coefficients that add up to 1,
# and they will also be in the (0, 1) range, plus they will be regularised
# alpha = 0 is ridge regression, but could use 1 for lasso (or anything in between)
mod3 <- cv.glmnet(X2, y2, intercept = FALSE, alpha = 0,
lower.limits = 0, upper.limits = 1)
# Coefficients from glmnet include the intercept even though we forced it to be zero
# which is why we need the (-1) in the below
coefs3 <- c(as.vector(coef(mod3)(-1)), 1 - sum(coef(mod3)))
f(coefs3, "Transform and ridge") # draw plot
Rezultatele sunt toate nenegative și se adaugă corect la unul:
> round(coefs3, 2) (1) 0.11 0.24 0.00 0.00 0.64 > sum(coefs3) (1) 1
Ele diferă semnificativ de coeficienții adevărați:
Vom reveni la asta.
În continuare este ceea ce am putea numi soluția „adecvată”, care este să tratăm problema ca pe cazul clasic al programării pătratice. Programarea pătratică a fost dezvoltată exact pentru a face față acestei situații – optimizați o funcție pătratică multivariată (cum ar fi cele mai mici pătrate într-o regresie) supuse unor constrângeri liniare (cum ar fi, coeficienții trebuie să fie nenegativi și să adună până la unul). Pachetul R al lui Hans Borchers („Funcții matematice numerice practice”) are o interfață ușor de înțeles pentru a face acest lucru fără a fi nevoie să inversați propriile matrice, așa că singurul lucru pe care trebuie să-l faceți este să exprimați constrângerile în mod corect. Iată cum facem asta cu problema noastră:
#---------------method 4 - quadratic programming-----------------
# This is based on dariober's solution, which itself is Elvis' solution
# but using the easier API from pracma::lsqlincon
# lsqlincon "solves linearly constrained linear least-squares problems"
# Equality constraint: We want the sum of the coefficients to be 1.
# I.e. Aeq x == beq
Aeq <- matrix(rep(1, ncol(X)), nrow= 1)
beq <- c(1)
# Lower and upper bounds of the parameters, i.e (0, 1)
lb <- rep(0, ncol(X))
ub <- rep(1, ncol(X))
# And solve:
coefs4 <- lsqlincon(X, y, Aeq= Aeq, beq= beq, lb= lb, ub= ub)
f(coefs4, "Quadratic programming")
La fel ca metoda anterioară, rezultatul nostru este nenegativ și se însumează la zero:
> round(coefs4, 2) (1) 0.17 0.44 0.00 0.00 0.38 > sum(coefs4) (1) 1
De data aceasta coeficienții arată puțin mai aproape de cei inițiali:
În cele din urmă, cealaltă metodă pe care am inventat-o a fost să merg all-in Bayesian și să specific pur și simplu în Stan că coeficienții trebuie să fie simplex și că au o distribuție Dirichlet. Aceasta mi se pare de fapt cea mai intuitivă soluție – ce modalitate mai bună de a face ceva un simplex decât să consideri că are o distribuție Dirichlet?
Așa că am scris acest program în Stan:
// Stan program for a regression where the coefficients are a simplex
// Probable use case is when y is some kind of weighted average of the X
// and you want to guarantee the weights add up to 1 (note - note sure I think
// this is a good idea but it is out there)
//
// Peter Ellis November 2024
data {
int n; // number of data items must match number rows of X
int k; // number of predictors must match number columns of X
matrix(n, k) X; // predictor matrix. Might include intercept (column of 1s)
vector(n) y; // outcome vector
}
parameters {
simplex(k) b; // coefficients for predictors
real sigma; // error scale
real alpha; // parameter for dirichlet distribution
}
model {
alpha ~ beta(2, 2);
b ~ dirichlet(rep_vector(alpha, k));
y ~ normal(X * b, sigma); // likelihood
}
… și l-a numit de la R cu asta:
#----------------method 5 - Bayesian model with a Dirichlet distribution------
stan_data <- list(
n = nrow(X),
k = ncol(X),
X = X,
y = as.vector(y)
)
mod5 <- stan("0281b-dirichlet-coefs.stan", data = stan_data)
coefs5 <- apply(extract(mod5, "b")$b,2, mean)
Aceasta este acum a treia metodă pentru care rezultatul nostru este corect tot nenegativ și se însumează la zero:
> round(coefs5, 2) (1) 0.17 0.45 0.00 0.00 0.38 > sum(coefs5) (1) 1
… și așa arată, practic foarte asemănător cu soluția de programare pătratică:
OK, care metodă este cea mai bună? Avem nevoie de un fel de benchmark pentru a le testa, mai bine decât coeficienții inițiali din procesul de generare a datelor în conformitate cu restricțiile fundamentale. Având în vedere că scopul final este de a construi o medie ponderată, dorim ceva care să aibă proporții similare cu coeficienții inițiali. Așadar, dorind doar să fiu pragmatic în acest sens, am creat un set „normalizat” de coeficienți eliminându-i pe cel negativ și scalând restul pentru a adăuga la unu.
Acest lucru ne permite să definim un criteriu de succes – metoda câștigătoare este cea care îndeplinește constrângerile (excluzând modelele 1 și 2) și se potrivește cel mai bine cu acele proporții. Iată un complot frumos de perechi care compară totul dintr-o privire:
Ceea ce vedem este că modelele 4 și 5 sunt câștigătoare clare, efectiv egale. Nu numai că sunt aproape identici, ci se potrivesc și cu coeficienții reali normalizați (corelații de 0,98 sau mai sus). Din păcate, metoda glmnet nu funcționează foarte bine la recuperarea proporțiilor inițiale. Poate că s-ar fi descurcat mai bine dacă aș lăsa coloanele lui X să fie corelate între ele, mai degrabă decât independente (care se ocupă cu multi-colinearitatea fiind terenul de acasă al regresiei nete elastice). Poate că, de asemenea, s-ar apropia de soluția de programare pătratică dacă aș forța cantitatea de regularizare la aproape zero. Dar, având deja două soluții bune, nu sunt înclinat să experimentez în continuare în acest sens.
Abordarea de programare pătratică a fost mai ușor de scris și de depanat și mai rapid de rulat decât modelarea explicită a distribuției Dirichlet în Stan, așa că pragmatic devine câștigătoare; singurul avertisment fiind acel semn de întrebare din ultimul paragraf despre ce se întâmplă dacă coloanele lui X sunt multi-coliniare.
Oh, să terminăm cu codul pentru a desena acel complot de perechi. Dintr-un motiv sau altul, am folosit o bază veche bună R pairs() pentru a face asta, mai degrabă GGally pachetul pe care îl folosesc în mod mai normal pentru matrice de dispersie. Îmi place controlul ușor pairs() vă oferă peste triunghiurile superioare și inferioare ale matricei scatterplot, prin definirea propriei funcții din mers. Deci, pentru cei mai vechi, iată acel cod R de bază:
#---------------comparisons------------
norm_coefs <- pmax(0, true_coefs) # remove negative values
norm_coefs <- norm_coefs / sum(norm_coefs) # force to ad dup to one
pairs(cbind('Actual datangenerating process' = true_coefs,
'Normalised tonadd to one and > 0' = norm_coefs,
"Ordinary least squares,nwon't add to 1, notnconstrained to (0,1)" = coefs1,
'Transform and OLS, add to 1nbut can be negative' = coefs2,
'Transform and ridge regression,nconstrain to within (0,1)' = coefs3,
'Quadratic programming' = coefs4,
'Dirichlet & Stan' = coefs5),
lower.panel = function(x, y){
abline(0, 1, col = "grey")
points(x, y, pch = 19)
},
upper.panel = function(x, y){
text(mean(range(x)), mean(range(y)), round(cor(x, y), 2),
col = "steelblue", font = 2)
},
main = "Comparison of different methods for creating a weighted average with weights adding up to one")
