Come eseguire i minimi quadrati non lineari con parametri condivisi in R?

Aug 22 2020

Vorrei eseguire la regressione dei minimi quadrati non lineari in R dove minimizzo simultaneamente i residui al quadrato di tre modelli (vedi sotto). Ora, i tre modelli condividono alcuni dei parametri, nel mio esempio, parametri be d.

C'è un modo per farlo con uno nls(), o, i pacchetti minpack.lmo nlsr?

Così, idealmente, vorrei generare la funzione obiettivo (la somma dei minimi quadrati di tutti i modelli insieme) e regredire tutti i parametri in una volta: a1, a2, a3, b, c1, c2, c3e d.

(Sto cercando di evitare di eseguire tre regressioni indipendenti e quindi di eseguire una media su be d.)

my_model <- function(x, a, b, c, d) {
  a * b ^ (x - c) + d
}

# x values
x <- seq(0, 10, 0.2)

# Shared parameters
b <- 2
d <- 10

a1 <- 1
c1 <- 1
y1 <- my_model(x,
               a = a1,
               b = b,
               c = c1,
               d = d) + rnorm(length(x))

a2 <- 2
c2 <- 5
y2 <- my_model(x,
               a = a2,
               b = b,
               c = c2,
               d = d) + rnorm(length(x))

a3 <- -2
c3 <- 3
y3 <- my_model(x,
               a = a3,
               b = b,
               c = c3,
               d = d) + rnorm(length(x))

plot(
  y1 ~ x,
  xlim = range(x),
  ylim = d + c(-50, 50),
  type = 'b',
  col = 'red',
  ylab = 'y'
)
lines(y2 ~ x, type = 'b', col = 'green')
lines(y3 ~ x, type = 'b', col = 'blue')

Risposte

2 G.Grothendieck Aug 22 2020 at 19:14

Di seguito eseguiamo nls(utilizzando un modello leggermente modificato) e nlxb(da nlsr) ma si nlxbferma prima della convergenza. Desidero che questi problemi entrambi diano comunque risultati che si adattano visivamente bene ai dati. Questi problemi suggeriscono che ci sono problemi con il modello stesso, quindi nella sezione Altro , guidati nlxbdall'output, mostriamo come correggere il modello fornendo un sottomodello del modello originale che si adatta facilmente ai dati con entrambi nlse nlxbe dà anche un buon adattamento . Alla fine nella sezione Note forniamo i dati in forma riproducibile.

nls

Assumendo la configurazione mostrata riproducibile nella Nota alla fine, riformulare il problema per l'algoritmo plineare nls definendo una matrice di destra le cui colonne moltiplicano ciascuno dei parametri lineari, rispettivamente a1, a2, a3 e d. plinear non richiede valori iniziali per chi semplifica il setup. Li riporterà rispettivamente come .lin1, .lin2, .lin3 e .lin4.

Per ottenere i valori di partenza abbiamo utilizzato un modello più semplice senza raggruppamento e una griglia di ricerca su b da 1 a 10 ec anche da 1 a 10 utilizzando nls2nel pacchetto con lo stesso nome. Abbiamo anche scoperto che nlsproduceva ancora errori ma, utilizzando absnella formula, come mostrato, è stato completato.

I problemi con il modello suggeriscono che c'è un problema fondamentale con esso e nella sezione Altro si discute come risolverlo.

xx <- c(x, x, x)
yy <- c(y1, y2, y3)

# startingi values using nls2
library(nls2)
fo0 <- yy ~ cbind(b ^ abs(xx - c), 1)
st0 <- data.frame(b = c(1, 10), c = c(1, 10))
fm0 <- nls2(fo0, start = st0, alg = "plinear-brute")

# run nls using starting values from above
g <- rep(1:3, each = length(x))   
fo <- yy ~ cbind((g==1) * b ^ abs(xx - c[g]), 
                 (g==2) * b ^ abs(xx - c[g]),  
                 (g==3) * b ^ abs(xx - c[g]), 
                 1) 
st <- with(as.list(coef(fm0)), list(b = b, c = c(c, c, c)))
fm <- nls(fo, start = st, alg = "plinear")

plot(yy ~ xx, col = g)
for(i in unique(g)) lines(predict(fm) ~ xx, col = i, subset = g == i)

fm

dando:

Nonlinear regression model
  model: yy ~ cbind((g == 1) * b^abs(xx - c[g]), (g == 2) * b^abs(xx -     c[g]), (g == 3) * b^abs(xx - c[g]), 1)
   data: parent.frame()
     b     c1     c2     c3  .lin1  .lin2  .lin3  .lin4 
 1.997  0.424  1.622  1.074  0.680  0.196 -0.532  9.922 
 residual sum-of-squares: 133

Number of iterations to convergence: 5 
Achieved convergence tolerance: 5.47e-06

(continua dopo la trama)

nlsr

Con nlsr sarebbe stato fatto così. Non era necessaria alcuna ricerca sulla griglia per i valori iniziali e absnemmeno l' aggiunta . I valori b e d sembrano simili alla soluzione nls ma gli altri coefficienti differiscono. Visivamente entrambe le soluzioni sembrano adattarsi ai dati.

D'altra parte dalla colonna JSingval vediamo che il giacobiano è carente di rango che lo ha fatto fermare e non produrre valori SE e la convergenza è dubbia (sebbene possa essere sufficiente dato che visivamente il grafico, non mostrato, buona vestibilità). Discuteremo come risolvere questo problema nella sezione Altro.

g1 <- g == 1; g2 <- g == 2; g3 <- g == 3
fo2 <- yy ~ g1 * (a1 * b ^ (xx - c1) + d) + 
            g2 * (a2 * b ^ (xx - c2) + d) + 
            g3 * (a3 * b ^ (xx - c3) + d)
st2 <- list(a1 = 1, a2 = 1, a3 = 1, b = 1, c1 = 1, c2 = 1, c3 = 1, d = 1)
fm2 <- nlxb(fo2, start = st2)
fm2

dando:

vn: [1] "yy" "g1" "a1" "b"  "xx" "c1" "d"  "g2" "a2" "c2" "g3" "a3" "c3"
no weights
nlsr object: x 
residual sumsquares =  133.45  on  153 observations
    after  16    Jacobian and  22 function evaluations
  name            coeff          SE       tstat      pval      gradient    JSingval   
a1               3.19575            NA         NA         NA    9.68e-10        4097  
a2               0.64157            NA         NA         NA   8.914e-11       662.5  
a3              -1.03096            NA         NA         NA  -1.002e-09       234.9  
b                1.99713            NA         NA         NA   -2.28e-08       72.57  
c1               2.66146            NA         NA         NA   -2.14e-09       10.25  
c2               3.33564            NA         NA         NA  -3.955e-11   1.585e-13  
c3                2.0297            NA         NA         NA  -7.144e-10   1.292e-13  
d                9.92363            NA         NA         NA  -2.603e-12   3.271e-14

Possiamo calcolare gli SE usando nls2 come secondo stadio, ma questo ancora non risolve il problema con tutto ciò che suggeriscono i valori singolari.

summary(nls2(fo2, start = coef(fm2), algorithm = "brute-force"))

dando:

Formula: yy ~ g1 * (a1 * b^(xx - c1) + d) + g2 * (a2 * b^(xx - c2) + d) + 
    g3 * (a3 * b^(xx - c3) + d)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
a1  3.20e+00   5.38e+05     0.0        1    
a2  6.42e-01   3.55e+05     0.0        1    
a3 -1.03e+00   3.16e+05     0.0        1    
b   2.00e+00   2.49e-03   803.4   <2e-16 ***
c1  2.66e+00   9.42e-02    28.2   <2e-16 ***
c2  3.34e+00   2.43e+05     0.0        1    
c3  2.03e+00   8.00e+05     0.0        1    
d   9.92e+00   4.42e+05     0.0        1    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.959 on 145 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: NA

Altro

Quando nlsha problemi ad adattare un modello, spesso suggerisce che c'è qualcosa che non va nel modello stesso. Giocandoci un po ', guidati dalla colonna JSingval nell'output nlsr sopra che suggerisce quei cparametri o dpotrebbe essere il problema, scopriamo che se fissiamo tutti ci valori dei parametri a 0, il modello è facile da adattare dati valori di partenza sufficientemente buoni e dà ancora una somma residua bassa dei quadrati.