title: Coding a prediction interval from a generalized additive model with a very large dataset
body: <p>I have a small dataset of people, thier locations, and whether or not they know each other.  It is a subset of a dataset with 1000 people.  Given that each person could know any other person, the number of potential links grows just under n^2.  I want to fit a model with the small subset, to get probability of linkage as a function of distance, so that I can perform simulations with the wider dataset.</p>
  <p>I've got two problems:</p>
  <ol>
  <li>I'm not sure how to create a prediction interval from a fitted GAM object.</li>
  <li>Generating a prediction interval using either <a href="http://r.789695.n4.nabble.com/Prediction-interval-with-GAM-td3460175.html" rel="nofollow noreferrer">posterior simulation</a> or using <a href="https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003585.html" rel="nofollow noreferrer">this technique</a> from R-sig-mixed is computationally prohibitive.</li>
  </ol>
  <p>The following is an example of my problem, generating the interval using the technique from R-sig-mixed.  Be warned that the last step will throw an error about not being able to allocate a huge vector, unless you're on a really impressive machine.</p>
  <pre><code>#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))))
  </code></pre>
  <p>I want to model links as a function of distance.  More fake data:</p>
  <pre><code>library(reshape)
  mdata &lt;- melt(as.data.frame(d), id=c(&quot;ID&quot;),measure.vars = colnames(d)[2:ncol(d)],variable.name=&quot;distance&quot;) 
  mdata$popularity = rnorm(25,sd=.3)
  colnames(mdata)[colnames(mdata)==&quot;variable&quot;] = &quot;knows&quot;
  colnames(mdata)[colnames(mdata)==&quot;value&quot;] = &quot;distance&quot;
  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))&lt;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
  </code></pre>
  <p>Now, I model the data:</p>
  <pre><code>library(mgcv)
  mod = gam(link~s(distance)+s(ID,bs=&quot;re&quot;,by=dum)+s(knows,bs=&quot;re&quot;,by=dum),data=mdata,family=binomial(link=&quot;logit&quot;))
  plot(mod,pages=1)
  summary(mod)
  </code></pre>
  <p>Now, I want to apply the fitted model to my main dataset:</p>
  <pre><code>x = runif(1000)*2
  y = runif(1000)*2
  d = cbind(ID = 1:1000,as.matrix(dist(data.frame(x,y))))
  mdata &lt;- melt(as.data.frame(d),id.vars = &quot;ID&quot;) 
  colnames(mdata)[colnames(mdata)==&quot;variable&quot;] = &quot;knows&quot;
  colnames(mdata)[colnames(mdata)==&quot;value&quot;] = &quot;distance&quot;
  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.
  </code></pre>
  <p>Some timekeeping tools...</p>
  <pre><code>tic &lt;- function(gcFirst = TRUE, type=c(&quot;elapsed&quot;, &quot;user.self&quot;, &quot;sys.self&quot;))
  {
     type &lt;- match.arg(type)
     assign(&quot;.type&quot;, type, envir=baseenv())
     if(gcFirst) gc(FALSE)
     tic &lt;- proc.time()[type]         
     assign(&quot;.tic&quot;, tic, envir=baseenv())
     invisible(tic)
  }
  
  toc &lt;- function()
  {
     type &lt;- get(&quot;.type&quot;, envir=baseenv())
     toc &lt;- proc.time()[type]
     tic &lt;- get(&quot;.tic&quot;, envir=baseenv())
     print(toc - tic)
     invisible(toc)
  }
  tic()
  p = predict(mod,newdata=mdata,type=&quot;response&quot;)
  toc()
  </code></pre>
  <p>Just predicting the point estimates takes 31 seconds on my machine.  Now to try to get the prediction intervals, first get the design matrix...</p>
  <pre><code>tic()
  Designmat = predict(mod,newdata=mdata,type=&quot;lpmatrix&quot;)
  toc()
  </code></pre>
  <p>That took me 47 seconds and froze my computer while it was working.</p>
  <p>Now here is the technique for getting a prediction interval that I found on R-sig-mixed...</p>
  <h1>CAUTION:  THE FOLLOWING CODE WILL ATTEMPT TO ALLOCATE A MASSIVE AMOUNT OF MEMORY AND MIGHT CRASH YOUR MACHINE.</h1>
  <pre><code>tic()
  predvar &lt;- diag(Designmat %*% vcov(mod) %*% t(Designmat))
  SE &lt;- sqrt(predvar) 
  SE2 &lt;- sqrt(predvar+mod$sig2) 
  tfrac &lt;- qt(0.975, mod$df.residual)
  interval = tfrac*SE2
  toc()
  
  &gt;Error: cannot allocate vector of size 7435.7 Gb
  </code></pre>
  <p>Is there another way???</p>