R trasforma nnls in lm

Question

Stavo cercando un modo per eseguire una regressione lineare sotto vincoli positivi, quindi ho trovato l'approccio nnls. Comunque mi stavo chiedendo come ottenere le stesse statistiche da nnls come quella fornita da un oggetto lm. Più specificamente il R-quadrato, i criteri akaike, i valori-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) 

C'è un modo per nuova stima in un quadro lm, in offset e che fissa il coefficiente non ha funzionato ... C'è un modo per ottenere queste statistiche? O un altro modo per creare un oggetto lm sotto vincoli positivi sul coefficiente?

Grazie Romain.

Answer 1

Quello che stai proponendo di fare è una pessima idea, tanto che sono riluttante a mostrarti come farlo. Il motivo è che per OLS, supponendo che i residui siano normalmente distribuiti con varianza costante, le stime dei parametri seguono una distribuzione t multivariata e possiamo calcolare limiti di confidenza e valori p nel solito modo.

Tuttavia, se facciamo NNLS sugli stessi dati, i residui non saranno normalmente ditributed, e le tecniche standard di calcolo dei valori di p, ecc produrranno rifiuti. Esistono metodi per stimare i limiti di confidenza sui parametri di un adattamento NNLS (per esempio, vedere this reference), ma sono approssimativi e di solito si basano su ipotesi abbastanza restrittive sul set di dati.

D'altra parte, sarebbe bello se alcune delle funzioni più basilari di un oggetto lm, come ad esempio predict(...), coeff(...), residuals(...), ecc lavorato anche per il risultato di una misura NNLS. Quindi, un modo per ottenere ciò che si utilizza è nls(...): solo perché un modello è lineare nei parametri non significa che non è possibile utilizzare i minimi quadrati non lineari per trovare i parametri. nls(...) offre l'opzione di impostare limiti inferiori (e superiori) sui parametri se si utilizza l'algoritmo 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 

Come si può vedere, il risultato dell'utilizzo nnls(...) e il risultato dell'utilizzo di nls(...) con lower-c(0,0,0) sono identici. Ma nls(...) produce un oggetto nls, che supporta (la maggior parte degli) gli stessi metodi di un oggetto lm. Così si può scrivere precict(fit), coef(fit), residuals(fit), AIC(fit) ecc È inoltre possibile scrivere summary(fit) e confint(fit)ma attenzione: i valori che si ottengono non sono significativi !!!

Per illustrare il punto sui residui, confrontiamo i residui per un OLS adatto a questi dati, con i residui per l'adattamento 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)) 

In questo set di dati, la parte stocastica della variabilità in y è N (0,1) di progettazione, in modo che i residui dalle OLS fit (terreno QQ a sinistra) sono normali . Ma i residui della stessa serie di dati montati utilizzando NNLS non sono lontanamente normali. Questo perché la vera dipendenza di y su x1 è -10, ma l'adattamento NNLS lo sta forzando a 0. Di conseguenza, la proporzione di residui molto grandi (sia positivi che negativi) è molto più alta di quanto ci si aspetterebbe dalla distribuzione normale.