Приглашаем всех экспертов по локальной регрессии и/или R!
Я столкнулся с ограничением стандартной функции loess в R и надеюсь на ваш совет. Текущая реализация поддерживает только 1–4 предиктора. Позвольте мне изложить наш сценарий приложения, чтобы показать, почему это может легко стать проблемой, как только мы захотим использовать глобально подходящие параметрические копеременные.
По сути, у нас есть пространственное искажение s(x,y), наложенное на ряд измерений z:
z_i = s(x_i,y_i) + v_{g_i}
Эти измерения z могут быть сгруппированы по одному и тому же базовому неискаженному значению измерения v для каждой группы g. Членство в группе g_i известно для каждого измерения, но базовые неискаженные значения измерения v_g для групп неизвестны и должны определяться с помощью (глобальной, а не локальной) регрессии.
Нам нужно оценить двумерный пространственный тренд s(x,y), который мы затем хотим удалить. В нашем приложении, скажем, есть 20 групп не менее 35 измерений в каждой, в самом простом сценарии. Измерения расположены случайным образом. Принимая за основу первую группу, получаем 19 неизвестных смещений.
Приведенный ниже код для данных игрушек (с пространственным трендом в одном измерении x) работает для двух или трех групп смещения.
К сожалению, вызов loess завершается ошибкой для четырех или более групп смещений с сообщением об ошибке
Error in simpleLoess(y, x, w, span, degree, parametric, drop.square,
normalize, :
only 1-4 predictors are allowed"
Я попытался обойти ограничение и получил
k>d2MAX in ehg136. Need to recompile with increased dimensions.
Насколько легко это будет сделать? Я нигде не могу найти определение d2MAX, и, похоже, это может быть жестко закодировано — ошибка, по-видимому, вызвана строкой № 1359 в loessf.f.
if(k .gt. 15) call ehg182(105)
В качестве альтернативы, кто-нибудь знает о реализации локальной регрессии с глобальными (параметрическими) группами смещения, которые можно было бы применить здесь?
Или есть лучший способ справиться с этим? Я попробовал lme с корреляционными структурами, но это оказалось намного медленнее.
Будем очень признательны за любые комментарии!
Большое спасибо,
Дэвид
###
#
# 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");
