Hay una aplicación de loess en R con más de 3 paramétrico de predictores o un truco para un efecto similar?

StackOverflow https://stackoverflow.com/questions/6370361

Pregunta

Llamando a todos los expertos en local de regresión y/o R!

He corrido en una limitación de la norma loess función en R y espero que usted tiene algunos consejos.La implementación actual sólo admite 1-4 predictores.Permítanme nuestra aplicación de escenario para mostrar por qué esto puede ser un problema tan pronto como queremos emplear globalmente ajuste paramétrico covariables.

Esencialmente, tenemos una distorsión espacial s(x,y) superpuesto sobre una serie de medidas z:

z_i = s(x_i,y_i) + v_{g_i}

Estas mediciones z se pueden agrupar por el mismo subyacente, sin distorsión de los valores de medición v para cada grupo g.La pertenencia al grupo g_i para cada medición, pero la subyacente sin distorsión de los valores de medición v_g para los grupos que no son conocidos y debe ser determinado por (global, no local) de regresión.

Necesitamos calcular las dos dimensiones espaciales tendencia s(x,y), que luego nos quieren quitar.En nuestra aplicación, decir que hay 20 grupos de al menos 35 de las mediciones de cada uno, en el más simple escenario.Las mediciones están colocadas al azar.Tomando el primer grupo como referencia, por lo tanto hay 19 desconocido compensaciones.

El código de abajo para juguete de datos (con una tendencia espacial en una dimensión x) obras para dos o tres desplazamiento de grupos.

Por desgracia, el loess llamada falla de cuatro o más desplazamiento de los grupos con el mensaje de error

Error in simpleLoess(y, x, w, span, degree, parametric, drop.square,
normalize,  :
  only 1-4 predictors are allowed"

He intentado reemplazar la restricción y consiguió

k>d2MAX in ehg136.  Need to recompile with increased dimensions.

Lo fácil que sería hacer?No puedo encontrar una definición de d2MAX en cualquier lugar, y parece que esta podría ser codificado -- el error es al parecer provocada por la línea #1359 en loessf.f

if(k .gt. 15)   call ehg182(105)

Alternativamente, ¿alguien sabe de una implementación de locales de regresión con global (paramétrico) desplazamiento de los grupos que se podría aplicar aquí?

O hay una mejor manera de lidiar con esto?He intentado lme con relación a las estructuras sino que parece ser mucho, mucho más lento.

Cualquier comentario sería muy apreciada!

Muchas gracias,
David

###
#
# loess with parametric offsets - toy data demo
#

x<-seq(0,9,.1);
x.N<-length(x);

o<-c(0.4,-0.8,1.2#,-0.2  # works for three but not four
     );  # these are the (unknown) offsets
o.N<-length(o);
f<-sapply(seq(o.N),
          function(n){
            ifelse((seq(x.N)<= n   *x.N/(o.N+1) &
                    seq(x.N)> (n-1)*x.N/(o.N+1)),
                    1,0);
          });
f<-f[sample(NROW(f)),];

y<-sin(x)+rnorm(length(x),0,.1)+f%*%o;
s.fs<-sapply(seq(NCOL(f)),function(i){paste('f',i,sep='')});
s<-paste(c('y~x',s.fs),collapse='+');
d<-data.frame(x,y,f)
names(d)<-c('x','y',s.fs);

l<-loess(formula(s),parametric=s.fs,drop.square=s.fs,normalize=F,data=d,
         span=0.4);
yp<-predict(l,newdata=d);
plot(x,y,pch='+',ylim=c(-3,3),col='red');  # input data
points(x,yp,pch='o',col='blue');           # fit of that

d0<-d; d0$f1<-d0$f2<-d0$f3<-0;
yp0<-predict(l,newdata=d0);
points(x,y-f%*%o);     # spatial distortion
lines(x,yp0,pch='+');  # estimate of that

op<-sapply(seq(NCOL(f)),function(i){(yp-yp0)[!!f[,i]][1]});

cat("Demo offsets:",o,"\n");
cat("Estimated offsets:",format(op,digits=1),"\n");
¿Fue útil?

Solución

¿Por qué no usa un modelo aditivo para esto? Paquete MGCV Manejará este tipo de modelo, si entiendo su pregunta, bien. Podría tener esto mal, pero el código que muestra se relaciona x ~ y, pero su pregunta menciona z ~ s (x, y) + g. Lo que muestro a continuación para gam() es para la respuesta z modelado por un suave espacial en x y y con g siendo estimado paramétricamente, con g almacenado como factor en el marco de datos:

require(mgcv)
m <- gam(z ~ s(x,y) + g, data = foo)

¿O he entendido mal lo que querías? Si desea publicar un pequeño fragmento de datos, puedo dar un ejemplo adecuado usando MGCV...?

Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top