Вычисление p-значений в ограниченных (неотрицательных) наименьших квадратах

Я использовал Matlab для выполнения неограниченных наименьших квадратов (обычных наименьших квадратов), и он автоматически выводит коэффициенты, тестовую статистику и p-значения.

Мой вопрос заключается в том, что при выполнении ограниченных наименьших квадратов (строго неотрицательных коэффициентов) он выводит только коэффициенты, БЕЗ тестовой статистики, p-значения.

Можно ли рассчитать эти значения для обеспечения значимости? И почему он не доступен напрямую в программном обеспечении (или любом другом программном обеспечении в этом отношении?)

Решение неотрицательных наименьших квадратов (NNLS) основано на алгоритме, который отличает его от обычных наименьших квадратов.

Алгебраическое выражение для стандартной ошибки (не работает)

С помощью регулярных наименьших квадратов вы можете выразить p-значения, используя t-критерий в сочетании с оценками дисперсии коэффициентов.

$\hat\theta$

$\theta$

Обратная матрица информации Фишера (не применимо)

Дисперсия / распределение оценки коэффициентов также асимптотически приближается к наблюдаемой информационной матрице Фишера :

Но я не уверен, хорошо ли это здесь применимо. Оценка NNLS не является объективной оценкой.

Метод Монте-Карло

Когда выражения становятся слишком сложными, вы можете использовать вычислительный метод для оценки ошибки. Используя метод Монте-Карло, вы моделируете распределение случайности эксперимента, моделируя повторы эксперимента (пересчитывая / моделируя новые данные), и на основе этого вы оцениваете дисперсию коэффициентов.

$\hat\theta$$\hat\sigma$ и на основе этого расчета новых данных (несколько тысяч повторений, но это зависит от того , насколько точной вы хотите) , из которого вы можете наблюдать Распределение (и вариация и производная от этого оценка погрешности) для коэффициентов. (и есть более сложные схемы для выполнения этого моделирования)

Если бы вы были в порядке, используя RI, думаю, вы могли бы также использовать bbmle «s mle2функцию для оптимизации наименьших квадратов функцию правдоподобия и вычислить 95% доверительные интервалы на неотрицательными коэффициентами NNLS. Кроме того, вы можете принять во внимание, что ваши коэффициенты не могут стать отрицательными, оптимизируя журнал ваших коэффициентов, так что в обратном преобразовании они никогда не станут отрицательными.

Вот числовой пример, иллюстрирующий этот подход, здесь в контексте деконволюции суперпозиции гауссовидных хроматографических пиков с гауссовским шумом на них: (любые комментарии приветствуются)

Сначала давайте смоделируем некоторые данные:

require(Matrix)
n = 200
x = 1:n
npeaks = 20
set.seed(123)
u = sample(x, npeaks, replace=FALSE) # peak locations which later need to be estimated
peakhrange = c(10,1E3) # peak height range
h = 10^runif(npeaks, min=log10(min(peakhrange)), max=log10(max(peakhrange))) # simulated peak heights, to be estimated
a = rep(0, n) # locations of spikes of simulated spike train, need to be estimated
a[u] = h
gauspeak = function(x, u, w, h=1) h*exp(((x-u)^2)/(-2*(w^2))) # shape of single peak, assumed to be known
bM = do.call(cbind, lapply(1:n, function (u) gauspeak(x, u=u, w=5, h=1) )) # banded matrix with theoretical peak shape function used
y_nonoise = as.vector(bM %*% a) # noiseless simulated signal = linear convolution of spike train with peak shape function
y = y_nonoise + rnorm(n, mean=0, sd=100) # simulated signal with gaussian noise on it
y = pmax(y,0)
par(mfrow=c(1,1))
plot(y, type="l", ylab="Signal", xlab="x", main="Simulated spike train (red) to be estimated given known blur kernel & with Gaussian noise")
lines(a, type="h", col="red")

Давайте теперь деконволюцию измеренного зашумленного сигнала yс помощью полосчатой ​​матрицы, содержащей сдвинутую копию известного ядра гауссовой формы размытияbM (это наша матрица ковариации / дизайна).

Во-первых, давайте деконволюцию сигнала с неотрицательными наименьшими квадратами:

library(nnls)
library(microbenchmark)
microbenchmark(a_nnls <- nnls(A=bM,b=y)$x) # 5.5 ms
plot(x, y, type="l", main="Ground truth (red), nnls estimate (blue)", ylab="Signal (black) & peaks (red & blue)", xlab="Time", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(-a_nnls, type="h", col="blue", lwd=2)
yhat = as.vector(bM %*% a_nnls) # predicted values
residuals = (y-yhat)
nonzero = (a_nnls!=0) # nonzero coefficients
n = nrow(bM)
p = sum(nonzero)+1 # nr of estimated parameters = nr of nonzero coefficients+estimated variance
variance = sum(residuals^2)/(n-p) # estimated variance = 8114.505

Теперь давайте оптимизируем отрицательное логарифмическое правдоподобие нашей цели потерь по Гауссу и оптимизируем логарифм ваших коэффициентов так, чтобы в обратном преобразовании они никогда не были отрицательными:

library(bbmle)
XM=as.matrix(bM)[,nonzero,drop=FALSE] # design matrix, keeping only covariates with nonnegative nnls coefs
colnames(XM)=paste0("v",as.character(1:n))[nonzero]
yv=as.vector(y) # response
# negative log likelihood function for gaussian loss
NEGLL_gaus_logbetas <- function(logbetas, X=XM, y=yv, sd=sqrt(variance)) {
  -sum(stats::dnorm(x = y, mean = X %*% exp(logbetas), sd = sd, log = TRUE))
}  
parnames(NEGLL_gaus_logbetas) <- colnames(XM)
system.time(fit <- mle2(
  minuslogl = NEGLL_gaus_logbetas, 
  start = setNames(log(a_nnls[nonzero]+1E-10), colnames(XM)), # we initialise with nnls estimates
  vecpar = TRUE,
  optimizer = "nlminb"
)) # takes 0.86s
AIC(fit) # 2394.857
summary(fit) # now gives log(coefficients) (note that p values are 2 sided)
# Coefficients:
#       Estimate Std. Error z value     Pr(z)    
# v10    4.57339    2.28665  2.0000 0.0454962 *  
# v11    5.30521    1.10127  4.8173 1.455e-06 ***
# v27    3.36162    1.37185  2.4504 0.0142689 *  
# v38    3.08328   23.98324  0.1286 0.8977059    
# v39    3.88101   12.01675  0.3230 0.7467206    
# v48    5.63771    3.33932  1.6883 0.0913571 .  
# v49    4.07475   16.21209  0.2513 0.8015511    
# v58    3.77749   19.78448  0.1909 0.8485789    
# v59    6.28745    1.53541  4.0950 4.222e-05 ***
# v70    1.23613  222.34992  0.0056 0.9955643    
# v71    2.67320   54.28789  0.0492 0.9607271    
# v80    5.54908    1.12656  4.9257 8.407e-07 ***
# v86    5.96813    9.31872  0.6404 0.5218830    
# v87    4.27829   84.86010  0.0504 0.9597911    
# v88    4.83853   21.42043  0.2259 0.8212918    
# v107   6.11318    0.64794  9.4348 < 2.2e-16 ***
# v108   4.13673    4.85345  0.8523 0.3940316    
# v117   3.27223    1.86578  1.7538 0.0794627 .  
# v129   4.48811    2.82435  1.5891 0.1120434    
# v130   4.79551    2.04481  2.3452 0.0190165 *  
# v145   3.97314    0.60547  6.5620 5.308e-11 ***
# v157   5.49003    0.13670 40.1608 < 2.2e-16 ***
# v172   5.88622    1.65908  3.5479 0.0003884 ***
# v173   6.49017    1.08156  6.0008 1.964e-09 ***
# v181   6.79913    1.81802  3.7399 0.0001841 ***
# v182   5.43450    7.66955  0.7086 0.4785848    
# v188   1.51878  233.81977  0.0065 0.9948174    
# v189   5.06634    5.20058  0.9742 0.3299632    
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# -2 log L: 2338.857 
exp(confint(fit, method="quad"))  # backtransformed confidence intervals calculated via quadratic approximation (=Wald confidence intervals)
#              2.5 %        97.5 %
# v10   1.095964e+00  8.562480e+03
# v11   2.326040e+01  1.743531e+03
# v27   1.959787e+00  4.242829e+02
# v38   8.403942e-20  5.670507e+21
# v39   2.863032e-09  8.206810e+11
# v48   4.036402e-01  1.953696e+05
# v49   9.330044e-13  3.710221e+15
# v58   6.309090e-16  3.027742e+18
# v59   2.652533e+01  1.090313e+04
# v70  1.871739e-189 6.330566e+189
# v71   8.933534e-46  2.349031e+47
# v80   2.824905e+01  2.338118e+03
# v86   4.568985e-06  3.342200e+10
# v87   4.216892e-71  1.233336e+74
# v88   7.383119e-17  2.159994e+20
# v107  1.268806e+02  1.608602e+03
# v108  4.626990e-03  8.468795e+05
# v117  6.806996e-01  1.021572e+03
# v129  3.508065e-01  2.255556e+04
# v130  2.198449e+00  6.655952e+03
# v145  1.622306e+01  1.741383e+02
# v157  1.853224e+02  3.167003e+02
# v172  1.393601e+01  9.301732e+03
# v173  7.907170e+01  5.486191e+03
# v181  2.542890e+01  3.164652e+04
# v182  6.789470e-05  7.735850e+08
# v188 4.284006e-199 4.867958e+199
# v189  5.936664e-03  4.236704e+06
library(broom)
signlevels = tidy(fit)$p.value/2 # 1-sided p values for peak to be sign higher than 1
adjsignlevels = p.adjust(signlevels, method="fdr") # FDR corrected p values
a_nnlsbbmle = exp(coef(fit)) # exp to backtransform
max(a_nnls[nonzero]-a_nnlsbbmle) # -9.981704e-11, coefficients as expected almost the same
plot(x, y, type="l", main="Ground truth (red), nnls bbmle logcoeff estimate (blue & green, green=FDR p value<0.05)", ylab="Signal (black) & peaks (red & blue)", xlab="Time", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(x[nonzero], -a_nnlsbbmle, type="h", col="blue", lwd=2)
lines(x[nonzero][(adjsignlevels<0.05)&(a_nnlsbbmle>1)], -a_nnlsbbmle[(adjsignlevels<0.05)&(a_nnlsbbmle>1)], 
      type="h", col="green", lwd=2)
sum((signlevels<0.05)&(a_nnlsbbmle>1)) # 14 peaks significantly higher than 1 before FDR correction
sum((adjsignlevels<0.05)&(a_nnlsbbmle>1)) # 11 peaks significant after FDR correction

Я не пытался сравнивать производительность этого метода относительно непараметрической или параметрической начальной загрузки, но это, безусловно, намного быстрее.

Я также был склонен думать, что я должен быть в состоянии рассчитать доверительные интервалы Вальда для неотрицательных nnlsкоэффициентов на основе наблюдаемой информационной матрицы Фишера, рассчитанной по логарифмически преобразованной шкале коэффициентов для обеспечения ограничений неотрицательности и оцененной наnnls оценкам.

Я думаю, что это так, и на самом деле должно быть формально идентично тому, что я использовал mle2выше:

XM=as.matrix(bM)[,nonzero,drop=FALSE] # design matrix
posbetas = a_nnls[nonzero] # nonzero nnls coefficients
dispersion=sum(residuals^2)/(n-p) # estimated dispersion (variance in case of gaussian noise) (1 if noise were poisson or binomial)
information_matrix = t(XM) %*% XM # observed Fisher information matrix for nonzero coefs, ie negative of the 2nd derivative (Hessian) of the log likelihood at param estimates
scaled_information_matrix = (t(XM) %*% XM)*(1/dispersion) # information matrix scaled by 1/dispersion
# let's now calculate this scaled information matrix on a log transformed Y scale (cf. stat.psu.edu/~sesa/stat504/Lecture/lec2part2.pdf, slide 20 eqn 8 & Table 1) to take into account the nonnegativity constraints on the parameters
scaled_information_matrix_logscale = scaled_information_matrix/((1/posbetas)^2) # scaled information_matrix on transformed log scale=scaled information matrix/(PHI'(betas)^2) if PHI(beta)=log(beta)
vcov_logscale = solve(scaled_information_matrix_logscale) # scaled variance-covariance matrix of coefs on log scale ie of log(posbetas) # PS maybe figure out how to do this in better way using chol2inv & QR decomposition - in R unscaled covariance matrix is calculated as chol2inv(qr(XW_glm)$qr)
SEs_logscale = sqrt(diag(vcov_logscale)) # SEs of coefs on log scale ie of log(posbetas)
posbetas_LOWER95CL = exp(log(posbetas) - 1.96*SEs_logscale)
posbetas_UPPER95CL = exp(log(posbetas) + 1.96*SEs_logscale)
data.frame("2.5 %"=posbetas_LOWER95CL,"97.5 %"=posbetas_UPPER95CL,check.names=F)
#            2.5 %        97.5 %
# 1   1.095874e+00  8.563185e+03
# 2   2.325947e+01  1.743600e+03
# 3   1.959691e+00  4.243037e+02
# 4   8.397159e-20  5.675087e+21
# 5   2.861885e-09  8.210098e+11
# 6   4.036017e-01  1.953882e+05
# 7   9.325838e-13  3.711894e+15
# 8   6.306894e-16  3.028796e+18
# 9   2.652467e+01  1.090340e+04
# 10 1.870702e-189 6.334074e+189
# 11  8.932335e-46  2.349347e+47
# 12  2.824872e+01  2.338145e+03
# 13  4.568282e-06  3.342714e+10
# 14  4.210592e-71  1.235182e+74
# 15  7.380152e-17  2.160863e+20
# 16  1.268778e+02  1.608639e+03
# 17  4.626207e-03  8.470228e+05
# 18  6.806543e-01  1.021640e+03
# 19  3.507709e-01  2.255786e+04
# 20  2.198287e+00  6.656441e+03
# 21  1.622270e+01  1.741421e+02
# 22  1.853214e+02  3.167018e+02
# 23  1.393520e+01  9.302273e+03
# 24  7.906871e+01  5.486398e+03
# 25  2.542730e+01  3.164851e+04
# 26  6.787667e-05  7.737904e+08
# 27 4.249153e-199 4.907886e+199
# 28  5.935583e-03  4.237476e+06
z_logscale = log(posbetas)/SEs_logscale # z values for log(coefs) being greater than 0, ie coefs being > 1 (since log(1) = 0) 
pvals = pnorm(z_logscale, lower.tail=FALSE) # one-sided p values for log(coefs) being greater than 0, ie coefs being > 1 (since log(1) = 0)
pvals.adj = p.adjust(pvals, method="fdr") # FDR corrected p values

plot(x, y, type="l", main="Ground truth (red), nnls estimates (blue & green, green=FDR Wald p value<0.05)", ylab="Signal (black) & peaks (red & blue)", xlab="Time", ylim=c(-max(y),max(y)))
lines(x,-y)
lines(a, type="h", col="red", lwd=2)
lines(-a_nnls, type="h", col="blue", lwd=2)
lines(x[nonzero][pvals.adj<0.05], -a_nnls[nonzero][pvals.adj<0.05], 
      type="h", col="green", lwd=2)
sum((pvals<0.05)&(posbetas>1)) # 14 peaks significantly higher than 1 before FDR correction
sum((pvals.adj<0.05)&(posbetas>1)) # 11 peaks significantly higher than 1 after FDR correction

Результаты этих вычислений и те, которые были возвращены mle2, почти идентичны (но намного быстрее), поэтому я думаю, что это правильно, и будет соответствовать тому, что мы неявно делали с mle2...

Простая установка ковариат с положительными коэффициентами в подбор с nnlsиспользованием регулярного линейного подбора модели, кстати, не работает, так как такое линейное моделирование не будет учитывать ограничения неотрицательности и, следовательно, приведет к бессмысленным доверительным интервалам, которые могут стать отрицательными. Это бумажный Джейсона Ли и Джонатана Тейлора «Точный вывод выбора после модели для предельного скрининга» также представляет метод, позволяющий сделать вывод выбора модели после неотрицательных коэффициентов nnls (или LASSO) и использует для этого усеченные гауссовские распределения. Я не видел какой-либо открыто доступной реализации этого метода для подгонки nnls - для подгонки LASSO есть селективный выводпакет, который делает что-то подобное. Если у кого-нибудь будет реализация, пожалуйста, дайте мне знать!

В вышеописанном методе можно также разделить данные в наборе обучения и проверки (например, нечетные и четные наблюдения) и вывести ковариаты с положительными коэффициентами из набора обучения, а затем вычислить доверительные интервалы и значения p из набора проверки. Это было бы немного более устойчивым к переоснащению, хотя это также привело бы к потере мощности, поскольку можно было бы использовать только половину данных. Я не делал этого здесь, потому что ограничение неотрицательности само по себе уже достаточно эффективно для защиты от переоснащения.

Чтобы быть более конкретным в отношении метода Монте-Карло, о котором упоминается @Martijn, вы можете использовать Bootstrap, метод повторной выборки, который включает выборку из исходных данных (с заменой) из множества наборов данных для оценки распределения оценочных коэффициентов и, следовательно, любой связанной статистики, включая доверительные интервалы и р-значения.

Широко используемый метод подробно описан здесь: Эфрон, Брэдли. «Методы начальной загрузки: еще один взгляд на складной нож». Прорывы в статистике. Springer, New York, NY, 1992. 569-593.

Это реализовано в Matlab, см. Https://www.mathworks.com/help/stats/bootstrp.html, в частности, раздел «Начальная загрузка модели регрессии».