R przekształcić nnls w lm

Question

Szukałem sposobu na regresję liniową pod pozytywnymi ograniczeniami, dlatego natknąłem się na podejście NNL. Jednak zastanawiałem się, jak mogę uzyskać te same statystyki od nnls, jak te dostarczone przez obiekt LM. Dokładniej: R-kwadrat, kryteria akaike, wartości p.

library(arm) 
library(nnls) 


data = runif(100*4, min = -1, max = 1) 
data = matrix(data, ncol = 4) 
colnames(data) = c("y", "x1", "x2", "x3") 
data = as.data.frame(data) 
data$x1 = -data$y 

A = as.matrix(data[,c("x1", "x2", "x3")]) 
b = data$y 

test = nnls(A,b) 
print(test) 

Czy istnieje sposób reestimate w ramach lm, offsetowego oraz ustalające współczynnik nie działa ... Czy istnieje sposób, aby uzyskać te statystyki? Lub inny sposób tworzenia obiektu Lm pod pozytywnymi ograniczeniami na współczynnik?

Dzięki Romain.

Answer 1

Co masz zamiar zrobić, to bardzo zły pomysł, tak bardzo, że nie chcę ci pokazać, jak to zrobić. Powodem jest to, że dla OLS, zakładając, że reszty są normalnie dystrybuowane ze stałą wariancją, wówczas oszacowania parametrów następują po wielu odmianach t-rozkładu i możemy obliczyć limity ufności i wartości p w zwykły sposób.

Jeśli jednak wykonać NNLS na tych samych danych, reszty nie będą normalnie ditributed i standardowe techniki do obliczenia wartości p itp spowoduje śmieci. Istnieją metody szacowania limitów ufności dla parametrów dopasowania NNLS (patrz na przykład this reference), ale są one przybliżone i zazwyczaj opierają się na dość restrykcyjnych założeniach dotyczących zbioru danych.

Z drugiej strony, byłoby miło, gdyby niektóre z bardziej podstawowych funkcji dla obiektu lm, takich jak predict(...), coeff(...), residuals(...) itp pracował również dla wyniku pasowania NNLS. Jednym ze sposobów osiągnięcia tego jest użycie nls(...): tylko dlatego, że model jest liniowy w parametrach, nie oznacza to, że nie można użyć nieliniowych najmniejszych kwadratów do znalezienia parametrów. nls(...) oferuje opcję ustawienia niższych (i wyższych) limitów parametrów, jeśli używasz algorytmu port.

set.seed(1) # for reproducible example 
data <- as.data.frame(matrix(runif(1e4, min = -1, max = 1),nc=4)) 
colnames(data) <-c("y", "x1", "x2", "x3") 
data$y <- with(data,-10*x1+x2 + rnorm(2500)) 

A <- as.matrix(data[,c("x1", "x2", "x3")]) 
b <- data$y 
test <- nnls(A,b) 
test 
# Nonnegative least squares model 
# x estimates: 0 1.142601 0 
# residual sum-of-squares: 88391 
# reason terminated: The solution has been computed sucessfully. 

fit <- nls(y~b.1*x1+b.2*x2+b.3*x3,data,algorithm="port",lower=c(0,0,0)) 
fit 
# Nonlinear regression model 
# model: y ~ b.1 * x1 + b.2 * x2 + b.3 * x3 
# data: data 
# b.1 b.2 b.3 
# 0.000 1.143 0.000 
# residual sum-of-squares: 88391 

Jak widać, wynik za pomocą nnls(...) a wynik za pomocą nls(...) z lower-c(0,0,0) są identyczne. Ale nls(...) tworzy obiekt nls, który obsługuje (większość) te same metody, co obiekt lm. Więc można napisać precict(fit), coef(fit), residuals(fit), AIC(fit) itd. Można również napisać summary(fit) i confint(fit)ale uważaj: wartości można dostać nie są miarodajne !!!

Aby zilustrować punkt dotyczący reszt, porównujemy reszty dla dopasowania OLS do tych danych, z resztami dla dopasowania NNLS.

par(mfrow=c(1,2),mar=c(3,4,1,1)) 
qqnorm(residuals(lm(y~.,data)),main="OLS"); qqline(residuals(lm(y~.,data))) 
qqnorm(residuals(fit),main="NNLS"); qqline(residuals(fit)) 

w tych bazach danych, część stochastyczny zmienności y oznacza N (0,1) za pomocą konstrukcji, a więc z resztek OLS Fit (QQ wykres po lewej stronie) jest normalne . Ale reszty z tego samego zestawu danych dopasowanego przy użyciu NNLS nie są zdalnie normalne. Dzieje się tak dlatego, że prawdziwa zależność od na x1 wynosi -10, ale dopasowanie NNLS wymusza jej 0. W rezultacie udział bardzo dużych reszt (zarówno dodatnich, jak i ujemnych) jest znacznie wyższy niż można by oczekiwać od rozkładu normalnego.