R help - Jun 2011 - Is there an implementation loess with more than 4 parametric predictors or a trick to similar effect?

Dear R experts,

I have a problem that is a related to the question raised in this earlier
post
    https://stat.ethz.ch/pipermail/r-help/2007-January/124064.html

My situation is different in that I have only 2 predictors (coordinates x,y)
for local regression but a number of global ("parametric") offsets
that I
need to consider.

Essentially, I have a spatial distortion overlaid over a number of
measurements. These measurements can be grouped by the same underlying value
for each group. The groups are known, but the values are not. We need to
estimate the spatial trend, which we then want to remove. In our
application, the spatial trend is two-dimensional (x,y), and there are about
20 groups of about 50 measurements each. The measurements are randomly
placed. Taking the first group as reference, there are thus 19 unknown
offsets.

The below code for toy data (spatial trend in one dimension x) works for two
or three offset groups, although I have not yet found out how to extract the
fitted values for the globally fit parameters (the offsets).

Unfortunately, the loess call fails for four or more offset groups with the
error message
"Error in simpleLoess(y, x, w, span, degree, parametric, drop.square,
normalize,  :
  only 1-4 predictors are allowed"

Does anyone know of an implementation of local regression with global
(parametric) offset groups that could be applied here, or is there a better
way of dealing with this?

Any comments would be greatly appreciated!

Best regards,
David Kreil.



###
#
# 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-f%*%o,ylim=c(-3,3));     # spatial distortion
lines(x,yp-f%*%o,pch='+');        # estimate of that
points(x,y,pch='+',col='red');    # input data
points(x,yp,pch='o',col='blue');  # fit of that

	[[alternative HTML version deleted]]

I have revised my earlier question to and would be grateful for any
comments!

----------

Subject: Is there an implementation of loess with more than 3 parametric
predictors or a trick to a similar effect?

Dear R experts,

I have a problem that is a related to the question raised in this earlier
post
    https://stat.ethz.ch/pipermail/r-help/2007-January/124064.html

My situation is different in that I have only 2 predictors (coordinates x,y)
for local regression but a number of global ("parametric") offsets
that I
need to consider.

Essentially, I have a spatial distortion overlaid over a number of
measurements. These measurements can be grouped by the same underlying
undistorted measurement value for each group. The groups are known, but the
values are not. We need to estimate the spatial trend, which we then want to
remove. In our application, the spatial trend is two-dimensional (x,y), and
there are about 20 groups of about 50 measurements each, in the most simple
scenario. The measurements are randomly placed. Taking the first group as
reference, there are thus 19 unknown offsets.

The below code for toy data (spatial trend in one dimension x) works for two
or three offset groups.

Unfortunately, the loess call fails for four or more offset groups with the
error message
"Error in simpleLoess(y, x, w, span, degree, parametric, drop.square,
normalize,  :
  only 1-4 predictors are allowed"

I tried overriding the restriction and got
"k>d2MAX in ehg136.  Need to recompile with increased dimensions."

How easy would that be to do? I cannot find a definition of d2MAX anywhere,
and it seems this might be hardcoded -- the error is apparently triggered by
line #1359 in loessf.f
      if(k .gt. 15)   call ehg182(105)

Alternatively, does anyone know of an implementation of local regression
with global (parametric) offset groups that could be applied here?

Or is there a better way of dealing with this? I tried lme with correlation
structures but that seems to be much, much slower.

Any comments would be greatly appreciated!

Best regards,
David Kreil.



###
#
# 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");

	[[alternative HTML version deleted]]