Esiste un'implementazione di Loess in R con più di 3 predittori parametrici o un trucco a un effetto simile?

Question

Chiamando tutti gli esperti regressione locale e/o R!

Ho incontrato una limitazione dello standard loess funzione in R E spero che tu abbia qualche consiglio. L'attuale implementazione Supporta solo 1-4 predittori. Lasciami impostare il nostro scenario di applicazione per mostrare perché questo può facilmente diventare un problema non appena vogliamo impiegare covariabili parametriche a livello globale.

In sostanza, abbiamo una distorsione spaziale s (x, y) sovrapposto su una serie di misurazioni z:

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

Queste misurazioni z può essere raggruppato dallo stesso valore di misurazione non distorto sottostante v per ogni gruppo g. L'appartenenza al gruppo G_i è noto per ciascuna misurazione, ma i valori di misurazione non distorti sottomistenti v_g per i gruppi non sono noti e dovrebbero essere determinati da (globale, non locale) regressione.

Dobbiamo stimare la tendenza spaziale bidimensionale s (x, y), che quindi vogliamo rimuovere. Nella nostra applicazione, diciamo che ci sono 20 gruppi di almeno 35 misurazioni ciascuno, nello scenario più semplice. Le misurazioni sono posizionate in modo casuale. Prendendo il primo gruppo come riferimento, ci sono quindi 19 offset sconosciuti.

Il codice seguente per i dati dei giocattoli (con una tendenza spaziale in una dimensione X) funziona per due o tre gruppi di offset.

Sfortunatamente, il loess La chiamata non riesce per quattro o più gruppi di offset con il messaggio di errore

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

Ho provato a prevalere sulla restrizione e ho ottenuto

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

Quanto sarebbe facile fare? Non riesco a trovare una definizione di d2max Ovunque, e sembra che questo possa essere codificato - L'errore è apparentemente attivato dalla linea #1359 in loessf.f

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

In alternativa, qualcuno sa di un'implementazione della regressione locale con gruppi di offset globali (parametrici) che potrebbero essere applicati qui?

O esiste un modo migliore per affrontare questo? Provai lme Con strutture di correlazione ma sembra essere molto, molto più lento.

Qualsiasi commento sarebbe molto apprezzato!

Grazie molto,
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");

Answer 1

Perché non usi un modello additivo per questo? Pacchetto mgcv Gestirerà questo tipo di modello, se capisco la tua domanda, benissimo. Potrei averlo sbagliato, ma il codice che mostri è relativo a x ~ y, ma la tua domanda menziona z ~ s (x, y) + g. Quello che mostro di seguito gam() è per la risposta z modellato da un liscio spaziale in x e y insieme a g essere stimato parametricamente, con g memorizzato come fattore nella cornice dei dati:

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

O ho frainteso quello che volevi? Se vuoi pubblicare un piccolo frammento di dati, posso fare un esempio adeguato usando mgcv...?