打电话给当地回归和/或R的所有专家!
我在R中遇到了标准黄土函数的局限性,希望你能给我一些建议。当前的实现只支持1-4个预测器。让我列出我们的应用场景,以说明为什么只要我们想要使用全局匹配的参数协变量,这就很容易成为一个问题。
本质上,我们有一个空间畸变s(x,y)覆盖在一些测量z上。
z_i = s(x_i,y_i) + v_{g_i}
这些测量值z可以按每个组g的相同的不失真测量值v进行分组。对于每个度量,组成员数g_i是已知的,但是组的基本不失真测量值v_g是未知的,应该通过(全局而不是局部)回归来确定。
我们需要估计二维空间趋势s(x,y),然后我们要删除它。在我们的应用程序中,假设在最简单的场景中,每组至少有35次度量。测量是随机放置的。因此,以第一组为参考,有19个未知的偏移量。
下面的玩具数据代码(具有一维x中的空间趋势)适用于两个或三个偏移组。
不幸的是,对于包含错误消息的四个或多个偏移组,黄土调用失败。
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的定义,这似乎是硬编码的--这个错误显然是由loessf.f中的第1359行触发的。
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");
发布于 2011-06-16 06:08:11
你为什么不用加性模型来做这个呢?Package mgcv将处理这类模型,如果我理解您的问题,很好。我可能错了,但您显示的代码与x~ y有关,但您的问题提到了z~ s(x,y) + g。下面我为gam()
显示的是用x
和y
中的空间平滑建模的响应z
,并对g
进行参数估计,并将g
作为数据框架中的一个因素存储:
require(mgcv)
m <- gam(z ~ s(x,y) + g, data = foo)
还是我误解了你想要的?如果您想发布一个小的数据片段,我可以使用mgcv.给出一个适当的示例。
https://stackoverflow.com/questions/6370361
复制相似问题