Spline all'interno lineari minimi quadrati a R

Question

consideri un lineare almeno modello piazze R, per esempio la seguente forma):

y ~ theta/(1 + exp(-(alpha + beta * x))) 

(mio vero problema ha diverse variabili e la funzione esterna non è logistico ma un un po 'più impegnato, questo è più semplice ma penso che se posso farlo il mio caso dovrebbe seguire quasi immediatamente)

Vorrei sostituire il termine "alpha + beta * x" con (diciamo) una spline cubica naturale .

ecco qualche codice per creare alcuni dati di esempio con una funzione non lineare all'interno della logistica:

set.seed(438572L) 
x <- seq(1,10,by=.25) 
y <- 8.6/(1+exp(-(-3+x/4.4+sqrt(x*1.1)*(1.-sin(1.+x/2.9))))) + rnorm(x, s=0.2) 

Senza la necessità di una logistica intorno ad esso, se ero in lm, ho potuto sostituire un termine lineare con un termine spline facilmente; così un modello di qualcosa di lineare, come questo:

lm(y ~ x) 

diventa allora

library("splines") 
lm(y ~ ns(x, df = 5)) 

generare valori stimati è semplice e ottenere i valori con l'ausilio di (per esempio) il pacchetto rms sembra abbastanza semplice predetto.

In effetti, il montaggio dei dati originali con quella spline basata su lm non è male, ma c'è un motivo per cui ho bisogno all'interno della funzione logistica (o piuttosto, l'equivalente nel mio problema).

Il problema con nls è necessario fornire i nomi per tutti i parametri (sono abbastanza contento di chiamarli dire (b1, ..., b5) per una spline in forma (e dire c1, ..., c6 per un'altra variabile - Dovrò essere in grado di crearne parecchi.

Esiste un modo abbastanza ragionevole per generare la formula corrispondente per nls in modo da poter sostituire il termine lineare all'interno della funzione non lineare con un spline?

Gli unici modi riesco a capire che ci potrebbe essere per farlo sono un po 'imbarazzante e goffo e non ben generalizzare senza scrivere un sacco di codice.

(modifica per chiarimenti) Per questo piccolo problema, posso farlo a mano ovviamente - scrivere un'espressione per prodotto interno di ogni variabile nella matrice generata da ns, volte il vettore dei parametri. Ma poi devo scrivere di nuovo l'intero termine di termine per ogni spline in ogni altra variabile, e ancora ogni volta che cambio il df in una qualsiasi delle spline, e di nuovo se voglio usare cs invece di ns. E poi quando voglio provare a fare qualche predizione (/ interpolazione), abbiamo una nuova serie di problemi da affrontare. Ho bisogno di continuare a farlo, più e più volte, e potenzialmente per un numero sostanzialmente maggiore di nodi, e su più variabili, per analisi dopo analisi - e mi chiedevo se esistesse un modo più semplice e semplice di scrivere ogni singolo termine, senza dover scrivere una grande quantità di codice. Riesco a vedere un modo abbastanza bull-at-a-gate per farlo che implicherebbe un bel po 'di codice per avere ragione, ma essendo R, sospetto che ci sia un modo molto più ordinato (o più probabilmente 3 o 4 modi più ordinari) che è semplicemente eludendomi. Da qui la domanda.

Pensavo di aver visto qualcuno fare qualcosa del genere in passato in un modo abbastanza carino, ma per la vita di me non riesco a trovarlo ora; Ho provato un sacco di volte per individuarlo.

[In particolare, mi piacerebbe in genere essere in grado di provare l'adattamento di più spline diverse in ciascuna variabile, per provare un paio di possibilità, per vedere se potevo trovare un modello semplice, ma ancora uno in cui la vestibilità è adeguata allo scopo (il rumore è davvero piuttosto basso, un po 'di pregiudizio nella vestibilità va bene per ottenere un risultato gradevole, ma solo fino a un certo punto). È più "trovare una bella, interpretabile, ma adeguata funzione di adattamento" di qualsiasi cosa che si avvicina all'inferenza e il data mining non è davvero un problema per questo problema.]

In alternativa, se questo sarebbe molto più facile in say gnm o ASSIST o uno degli altri pacchetti, sarebbe una conoscenza utile, ma poi alcuni suggerimenti su come procedere sul problema del giocattolo sopra con loro aiuterebbero.

Answer 1

ns genera effettivamente una matrice di predittori. Quello che puoi fare è dividere quella matrice in singole variabili e darle da mangiare a nls.

m <- ns(x, df=5) 
df <- data.frame(y, m) # X-variables will be named X1, ... X5 
# starting values should be set as appropriate for your data 
nls(y ~ theta * plogis(alpha + b1*X1 + b2*X2 + b3*X3 + b4*X4 + b5*X5), data=df, 
     start=list(theta=1, alpha=0, b1=1, b2=1, b3=1, b4=1, b5=1)) 

ETA: ecco un tentativo di automatizzare questo per diversi valori di df. Ciò costruisce la formula utilizzando il testo munging e quindi utilizza do.call per chiamare nls. Avvertenza: non testata.

my.nls <- function(x, y, df) 
{ 
    m <- ns(x, df=df) 
    xn <- colnames(m) 
    b <- paste("b", seq_along(xn), sep="") 
    fm <- formula(paste("y ~ theta * plogis(1 + alpha + ", paste(b, xn, sep="*", 
      collapse=" + "), ")", sep="")) 
    start <- c(1, 1, rep(1, length=length(b))) 
    names(start) <- c("theta", "alpha", b) 
    do.call(nls, list(fm, data=data.frame(y, m), start=start)) 
} 

Answer 2

Una realizzazione a cui sono giunto mentre chiarendo la mia domanda mi ha fatto vedere che c'è un modo meno goffo di quello che avevo visto prima.

Anche con un po 'di snellimento ovvio che può entrare, questo è ancora un po' inelegante per il mio occhio, ma almeno sopportabile abbastanza da usare su una base ripetuta, quindi la considero una risposta adeguata. Sono ancora interessato a un modo più ordinato di questo qui sotto.

Il trucco di Hong Ooi di usare data.frame sulla matrice generata da ns per auto-nominare le colonne è abbastanza carino e l'ho usato sotto. Probabilmente userò paste per costruirle in generale, perché ho diverse variabili con cui giocare.

Supponendo che i dati di set-up data alla questione -

lin.expr <- function(p,xn) { 
    pn<-paste(p, 1:length(xn), sep = "") 
    paste(paste(pn,xn,sep=" * "),collapse=" + ") 
    } 


m <- ns(x, df=3) 
mydf <- data.frame(y, m) # X-variables will be named X1, X2, ... 
xn <- names(mydf)[2:dim(mydf)[2]] 

nspb <- lin.expr("b",xn) 

c.form <- paste("y ~ theta * plogis(a + ",nspb,")",sep="") 
stl <- list(theta=2, a=-5,b1=10, b2=10, b3=10) 
nls(c.form, data=mydf, start= stl) 

La mia formula attuale avrà diversi termini come nspb. Miglioramenti sostanziali apprezzati; Preferirei non scegliere la mia risposta, ma suppongo che la sceglierò se non ci sarà più nulla in un giorno o due.

modifica: L'aggiunta di Hong Ooi (che è stata pubblicata mentre inserivo il mio in e utilizza idee simili, ma aggiungo un paio di bei extra) praticamente lo fa; è una risposta accettabile, quindi l'ho verificato.