2012-02-04 21 views
5

Considere un modelo no lineal de mínimos cuadrados en R, por ejemplo de la siguiente forma):Splines dentro de mínimos cuadrados no lineales en R

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

(mi verdadero problema tiene varias variables y la función externa no es logístico sino una un poco más involucrado, este es más simple, pero creo que si puedo hacer esto, mi caso debería seguir casi de inmediato)

Me gustaría reemplazar el término "alpha + beta * x" con (digamos) una spline cúbica natural .

aquí algo de código para crear unos datos de ejemplo de una función no lineal dentro de la logística:

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) 

Sin la necesidad de una logística alrededor de ella, si yo estaba en la película, que podría sustituir un término lineal con una término spline fácilmente; por lo que un modelo lineal algo como esto:

lm(y ~ x) 

se convierte entonces en

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

generar valores ajustados es simple y consiguiendo valores con la ayuda de (por ejemplo ) el paquete rms parece bastante simple predijo.

De hecho, ajustar los datos originales con ese ajuste spline basado en lm no es tan malo, pero hay una razón por la que lo necesito dentro de la función logística (o más bien, el equivalente en mi problema).

El problema con nls es que necesito proporcionar nombres para todos los parámetros (estoy bastante contento de llamarlos say (b1, ..., b5) para un ajuste spline (y decir c1, ..., c6 para otra variable - Necesitaré poder hacer varios de ellos).

Hay una manera razonablemente clara de generar la fórmula correspondiente para nls para que pueda reemplazar el término lineal dentro de la función no lineal con un spline?

la única manera que se me ocurre que podría haber para hacerlo son un poco incómodo y torpe y no se generalizan muy bien sin necesidad de escribir un montón de código.

(edición para aclaración) Para este pequeño problema, puedo hacerlo a mano, por supuesto, escriba una expresión para el producto interno de cada variable en la matriz generada por ns, multiplicado por el vector de parámetros. Pero luego tengo que escribir todo el asunto término por término nuevamente para cada spline en cada otra variable, y de nuevo cada vez que cambio el df en cualquiera de las splines, y de nuevo si quiero usar cs en lugar de ns. Y luego, cuando quiero intentar hacer algo de predicción (/ interpolación), recibimos una gran cantidad de nuevos problemas para tratar. Tengo que seguir haciéndolo, una y otra vez, y potencialmente por una cantidad de nudos considerablemente mayor, y por varias variables, para analizarlas después del análisis, y me pregunté si habría una manera más ordenada y simple que escribir cada término individual, sin tener que escribir una gran cantidad de código. Puedo ver una manera bastante fácil de hacerlo que implicaría un poco de código para hacerlo bien, pero al ser R, sospecho que hay una manera mucho más ordenada (o más probablemente 3 o 4 formas más ordenadas) que es simplemente eludiéndome. De ahí la pregunta.

Pensé que había visto a alguien hacer algo como esto en el pasado de una manera bastante agradable, pero por mi vida no puedo encontrarlo ahora; He intentado muchas veces localizarlo.

[Más en particular, en general, me gustaría poder probar el ajuste de varias splines diferentes en cada variable, para intentar un par de posibilidades, para ver si puedo encontrar un modelo simple, pero aún así uno donde el ajuste es adecuado para este propósito (el ruido es realmente bastante bajo, algunos sesgos en el ajuste están bien para lograr un buen resultado suave, pero solo hasta cierto punto). Es más 'encontrar una función agradable, interpretable, pero adecuada' que cualquier inferencia que se aproxime y la minería de datos no es realmente un problema para este problema.]

Alternativamente, si esto fuera mucho más fácil en decir gnm o ASSIST o uno de los otros paquetes, eso sería conocimiento útil, pero luego algunos consejos sobre cómo proceder con el problema del juguete anterior con ellos ayudaría.

Respuesta

9

ns en realidad genera una matriz de predictores. Lo que puede hacer es dividir esa matriz en variables individuales y alimentarlas al 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: aquí tienes una vez para automatizar esto para diferentes valores de df. Esto construye la fórmula usando text munging, y luego usa do.call para llamar al nls. Advertencia: no probado.

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)) 
} 
+0

@Glen_b: Ok, he editado mi respuesta; ver si esto ayuda. –

2

Una realización que encontré al aclarar mi propia pregunta me hizo ver que hay una manera menos torpe que había visto antes.

Incluso con un poco de racionalización obvia que puede entrar, esto sigue siendo un poco poco elegante a mi vista, pero al menos lo suficientemente soportable como para usarlo repetidamente, por lo que lo considero una respuesta adecuada. Estoy todavía interesado en una manera más nítida que esta a continuación.

El truco de Hong Ooi de usar data.frame en la matriz generada por ns para auto-nombrar las columnas es algo lindo y lo he usado a continuación. Probablemente usaré pegar para compilarlos en general, porque tengo varias variables para jugar.

Suponiendo que el conjunto de datos que se ha dado en la pregunta -

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) 

Mi fórmula real tendrá varios términos como nspb. Mejoras sustantivas apreciadas; Preferiría no elegir mi propia respuesta, pero supongo que la escogeré si no hay nada más en un día o dos.

editar: la adición de Hong Ooi (que fue publicada mientras escribía la mía y usa ideas similares, pero agrega un par de extras) básicamente lo hace; es una respuesta aceptable, así que lo he comprobado.

Cuestiones relacionadas