R transformieren nnls in lm

Question

Ich war auf der Suche nach einer Möglichkeit, eine lineare Regression unter positiven Einschränkungen zu tun, kam daher über die nnls Ansatz. Wie auch immer ich mich gefragt habe, wie ich die gleichen Statistiken von den nnls bekommen konnte wie die von einem lm-Objekt. Genauer gesagt das R-Quadrat, das Akaike-Kriterium, die P-Werte.

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) 

Gibt es eine Möglichkeit in einem lm Rahmen reestimate unter Verwendung von Offset und zur Festsetzung der Koeffizienten nicht funktioniert ... Gibt es eine Möglichkeit, diese Statistiken zu erhalten? Oder eine andere Möglichkeit, ein LM-Objekt unter positiven Bedingungen für den Koeffizienten zu erstellen?

Danke Romain.

Answer 1

Was Sie vorhaben zu tun, ist eine massiv schlechte Idee, so sehr, dass ich widerwillig bin, Ihnen zu zeigen, wie es geht. Der Grund dafür ist, dass die Parameterschätzungen für OLS unter der Annahme, dass die Residuen normalerweise mit konstanter Varianz verteilt sind, einer multivariaten t-Verteilung folgen und Konfidenzgrenzen und p-Werte auf die übliche Weise berechnet werden können.

Wenn wir jedoch NNLS auf denselben Daten durchzuführen, die Residuen werden normalerweise nicht und die Standardtechniken für die Berechnung der p-Werte werden ditributed usw. produzieren garbage. Es gibt Methoden zur Schätzung von Konfidenzgrenzen für die Parameter einer NNLS-Anpassung (siehe z. B. this reference), aber sie sind Näherungswerte und beruhen normalerweise auf ziemlich restriktiven Annahmen über den Datensatz.

Auf der anderen Seite wäre es schön, wenn einige der Grundfunktionen für ein lm Objekt, wie predict(...), coeff(...), residuals(...) usw. auch für das Ergebnis eines NNLS fit gearbeitet. Ein Weg, um das zu erreichen, ist die Verwendung nls(...): Nur weil ein Modell in den Parametern linear ist, bedeutet das nicht, dass Sie nicht-lineare kleinste Quadrate verwenden können, um die Parameter zu finden. nls(...) bietet die Option, niedrigere (und obere) Grenzwerte für die Parameter festzulegen, wenn Sie den Algorithmus port verwenden.

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 

Wie Sie sehen können, ist das Ergebnis von nnls(...) verwenden und das Ergebnis nls(...) mit lower-c(0,0,0) verwenden, sind identisch. Aber nls(...) produziert ein nls Objekt, das (die meisten) die gleichen Methoden wie ein lm Objekt unterstützt. So können Sie schreiben precict(fit), coef(fit), residuals(fit), AIC(fit) etc. Sie auch summary(fit) schreiben und confint(fit)aber Vorsicht: Die Werte, die Sie erhalten, sind nicht sinnvoll !!!

Um den Punkt zu den Residuen zu verdeutlichen, vergleichen wir die Residuen für eine OLS, die mit diesen Daten übereinstimmt, mit den Residuen für die NNLS-Anpassung.

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)) 

In diesem Datensatz der stochastische Teil der Variabilität in y ist N (0,1) von Design, so dass die aus der Residuen OLS passen (QQ-Plot auf der linken Seite) sind normal . Aber die Reste aus dem gleichen Datensatz, die mit NNLS angepasst wurden, sind nicht im entferntesten normal. Dies liegt daran, dass die wahre Abhängigkeit von y auf x1-10 ist, aber die NNLS-Anpassung zwingt sie zu 0. Folglich ist der Anteil sehr großer Residuen (sowohl positiv als auch negativ) viel höher, als von der normalen Verteilung erwartet würde.