제한된 (음수가 아닌) 최소 제곱에서 p- 값 계산 p- 값을 자동으로 출력합니다. 내 질문은

나는 Matlab을 사용하여 제한되지 않은 최소 제곱 (일반 최소 제곱)을 수행했으며 계수, 테스트 통계 및 p- 값을 자동으로 출력합니다.

내 질문은 제한된 최소 제곱 (음이 아닌 음의 계수)을 수행하면 테스트 통계없이 p- 값없이 계수 만 출력합니다.

유의성을 확보하기 위해 이러한 값을 계산할 수 있습니까? 소프트웨어에서 직접 사용할 수없는 이유는 무엇입니까?

답변

음이 아닌 최소 제곱 (NNLS)을 해결하는 것은 정규 최소 제곱과 다른 알고리즘 을 기반으로 합니다 .

표준 오차에 대한 대수 표현 (작동하지 않음)

규칙적으로 최소 제곱을 사용하면 계수 분산에 대한 추정값과 함께 t- 검정을 사용하여 p- 값을 표현할 수 있습니다.

계수의 추정의 샘플의 분산에 대한이 식 인 오차의 분산 일반적 알 것이다 있지만 잔차를 이용하여 추정 될 수있다 . 이 식은 측정면에서의 계수에 대해 식부터 수학적으로 유도 될 수있다 $\hat{θ}$

\hat{θ}

$\hat\theta$

V ㅏ 아르 자형 (\hat{θ}) = σ^{2} ({엑스}^{티} 엑스)^{- 1}

$Var(\hat\theta) = \sigma^2(X^TX)^{-1}$ $σ$

σ

$\sigma$ $y$

와이

$y$

\hat{θ} = ({엑스}^{티} 엑스)^{- 1} {엑스}^{티} 와이

$\hat\theta = (X^TX)^{-1} X^T y$

$θ$

θ

$\theta$

피셔 정보 매트릭스의 역수 (해당 사항 없음)

계수 추정치의 분산 / 분포는 또한 관측 된 Fisher 정보 매트릭스에 무조건 접근 합니다 .

(\hat{θ} - θ) \overset{디}{\to} 엔 (0, 나는 (\hat{θ}))

$(\hat\theta-\theta) \xrightarrow{d} N(0,\mathcal{I}(\hat\theta))$

그러나 이것이 여기에 잘 적용되는지 확실하지 않습니다. NNLS 추정치는 바이어스되지 않은 추정치가 아닙니다.

몬테 카를로 방법

식이 너무 복잡해지면 계산 방법을 사용하여 오류를 추정 할 수 있습니다. 으로 몬테카를로 방법 는 (새로운 데이터를 모델링 / 재 계산)과 이것에 근거하여 계수의 변동을 추정 실험을 반복하여 시뮬레이션 실험의 랜덤 분포를 시뮬레이션.

$\hat{θ}$

\hat{θ}

$\hat\theta$ $\hat{σ}$

\hat{σ}

$\hat\sigma$

답변

RI를 사용해도 괜찮다면 bbmle의 mle2함수를 사용 하여 최소 제곱 가능성 함수를 최적화하고 음이 아닌 nnls 계수에서 95 % 신뢰 구간을 계산할 수 있다고 생각합니다 . 또한 계수 로그를 최적화하여 계수가 마이너스가 될 수 없으므로 역변환 된 스케일에서는 절대 마이너스가 될 수 없습니다.

다음은 가우시안 노이즈가있는 가우스 모양의 크로마토 그래피 피크 중첩을 해제하는 맥락에서이 접근 방식을 설명하는 숫자 예제입니다. (모든 의견 환영)

먼저 일부 데이터를 시뮬레이트합시다.

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(이것은 공변량 / 설계 행렬)의 시프트 된 사본을 포함하는 밴딩 된 행렬로 측정 된 노이즈 신호 를 디컨 볼 루트하겠습니다 .

먼저 음이 아닌 최소 제곱으로 신호를 deconvolute하자

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관측 된 Fisher 정보 매트릭스를 기반으로하는 음이 아닌 계수에 대한 Wald 신뢰 구간을 계산할 수 있어야하며, 음 이 아닌 구속 조건을 시행하기 위해 로그 변환 된 계수 척도로 계산하고 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정규 선형 모형 적합 btw를 사용하여 적합에 양의 계수를 갖는 공변량을 재기동하는 것은 효과가 없습니다. 그러한 선형 모형 적합은 비 음성 구속 조건을 고려하지 않으므로 음이 될 수있는 무감각 신뢰 구간이 발생하기 때문입니다. Jason Lee & Jonathan Taylor 의이 논문 “마진 스크리닝에 대한 정확한 포스트 모델 선택 추론” 은 음이 아닌 nnl (또는 LASSO) 계수에 대한 포스트 모델 선택 추론을 수행하는 방법을 제시하고이를 위해 잘린 가우스 분포를 사용합니다. nnls에 맞는이 방법의 공개적으로 사용 가능한 구현을 보지 못했습니다 .- LASSO에 대해서는 selectiveInference가 있습니다.그런 일을하는 패키지. 누군가 구현을 원한다면 알려주세요!

위의 방법에서 훈련 & 검증 세트 (예 : 홀수 & 짝수 관측)의 데이터를 분할하고 훈련 세트에서 양의 계수로 공변량을 유추 한 다음 검증 세트에서 신뢰 구간 및 p 값을 계산할 수 있습니다. 데이터의 절반 만 사용하므로 전력 손실이 발생할 수 있지만 과적 합에 대해 조금 더 저항력이 있습니다. 비 음성 제약 자체가 이미 과적 합을 막는 데 상당히 효과적이기 때문에 여기서는 그렇게하지 않았습니다.

답변

@Martijn이 언급 한 Monte Carlo 방법에 대해 좀 더 구체적으로 말하면, 추정 계수의 분포와 관련 통계량을 추정하기 위해 원래 데이터 (대체로) 여러 데이터 세트에서 샘플링하는 리샘플링 방법 인 Bootstrap을 사용할 수 있습니다. 신뢰 구간 및 p- 값을 포함합니다.

널리 사용되는 방법은 Efron, Bradley입니다. “부트 스트랩 방법 : 또 다른 잭나이프 살펴보기” 통계의 혁신. 스프링거, 뉴욕, 뉴욕, 1992. 569-593.

Matlab이 구현했습니다. https://www.mathworks.com/help/stats/bootstrp.html 특히 회귀 모델 부트 스트랩 섹션을 참조 하십시오 .

How IT

언제든지 물어보세요.

제한된 (음수가 아닌) 최소 제곱에서 p- 값 계산 p- 값을 자동으로 출력합니다. 내 질문은

답변

표준 오차에 대한 대수 표현 (작동하지 않음)

피셔 정보 매트릭스의 역수 (해당 사항 없음)

몬테 카를로 방법

답변

답변

답변