2

我有一个关于人、他们的位置以及他们是否彼此认识的小型数据集。它是包含 1000 人的数据集的子集。鉴于每个人都可以认识任何其他人，潜在链接的数量会增长到 n^2 以下。我想用小子集拟合模型，以获得作为距离函数的链接概率，以便我可以使用更广泛的数据集执行模拟。

我有两个问题：

  1. 我不确定如何从拟合的 GAM 对象创建预测区间。
  2. 使用后验模拟或使用R-sig-mixed 中的这种技术生成预测区间在计算上是令人望而却步的。

以下是我的问题的一个示例，使用 R-sig-mixed 的技术生成间隔。请注意，最后一步将抛出一个错误，即无法分配一个巨大的向量，除非你在一台非常令人印象深刻的机器上。

#Some fake location data
set.seed(13)
x = runif(50)*2
y = runif(50)*2
d = cbind(ID = 1:50,as.matrix(dist(data.frame(x,y))))

我想将链接建模为距离的函数。更多虚假数据：

library(reshape)
mdata <- melt(as.data.frame(d), id=c("ID"),measure.vars = colnames(d)[2:ncol(d)],variable.name="distance") 
mdata$popularity = rnorm(25,sd=.3)
colnames(mdata)[colnames(mdata)=="variable"] = "knows"
colnames(mdata)[colnames(mdata)=="value"] = "distance"
mdata = subset(mdata,ID!=knows)
a = exp(1/(mdata$distance/runif(nrow(mdata))^mdata$distance)+mdata$popularity+rnorm(nrow(mdata),sd=.001))
mdata$prlink = a/(1+a)
with(mdata,plot(distance,prlink))
mdata$link = runif(nrow(mdata))<mdata$prlink
mdata$ID = as.factor(mdata$ID)
mdata$knows = as.factor(mdata$knows)
mdata$dum=1 #this facilitates predicting from the population of the model, later

现在，我对数据进行建模：

library(mgcv)
mod = gam(link~s(distance)+s(ID,bs="re",by=dum)+s(knows,bs="re",by=dum),data=mdata,family=binomial(link="logit"))
plot(mod,pages=1)
summary(mod)

现在，我想将拟合模型应用于我的主数据集：

x = runif(1000)*2
y = runif(1000)*2
d = cbind(ID = 1:1000,as.matrix(dist(data.frame(x,y))))
mdata <- melt(as.data.frame(d),id.vars = "ID") 
colnames(mdata)[colnames(mdata)=="variable"] = "knows"
colnames(mdata)[colnames(mdata)=="value"] = "distance"
mdata = subset(mdata,ID!=knows)
mdata$dum=0; mdata$ID=1; mdata$knows=2 #These are needed for prediction, even though I am predicting from the population of the model, not one of the levels.

一些计时工具...

tic <- function(gcFirst = TRUE, type=c("elapsed", "user.self", "sys.self"))
{
   type <- match.arg(type)
   assign(".type", type, envir=baseenv())
   if(gcFirst) gc(FALSE)
   tic <- proc.time()[type]         
   assign(".tic", tic, envir=baseenv())
   invisible(tic)
}

toc <- function()
{
   type <- get(".type", envir=baseenv())
   toc <- proc.time()[type]
   tic <- get(".tic", envir=baseenv())
   print(toc - tic)
   invisible(toc)
}
tic()
p = predict(mod,newdata=mdata,type="response")
toc()

在我的机器上仅预测点估计需要 31 秒。现在尝试获取预测区间，首先获取设计矩阵......

tic()
Designmat = predict(mod,newdata=mdata,type="lpmatrix")
toc()

这花了我 47 秒，并在我的电脑工作时冻结了它。

现在这里是获取我在 R-sig-mixed 上找到的预测区间的技术......

注意：以下代码将尝试分配大量内存并可能导致您的机器崩溃。

tic()
predvar <- diag(Designmat %*% vcov(mod) %*% t(Designmat))
SE <- sqrt(predvar) 
SE2 <- sqrt(predvar+mod$sig2) 
tfrac <- qt(0.975, mod$df.residual)
interval = tfrac*SE2
toc()

>Error: cannot allocate vector of size 7435.7 Gb

有没有别的办法？？？

4

1 回答 1

1

您需要避免计算Designmat %*% vcov(mod) %*% t(Designmat). 毕竟你只需要对角线。试试这个：

tmp <- Designmat %*% vcov(mod)

library(compiler)
diagMult <- cmpfun(function(m1, m2) sapply(seq_len(nrow(m1)), 
                                            function(i) m1[i,] %*% m2[,i]))
predvar <-  diagMult(tmp, t(Designmat))

（未彻底测试。如果某些包中还没有编译版本，则该功能应使用 Rcpp 实现以提高速度。）

于 2013-09-20T08:13:20.957 回答