I am aware that the "prediction interval", as defined in most textbooks on linear models, is focused on the uncertainty in the model being fit and is used to estimate an output prediction range for exact inputs. However, as is normally the case, what about situations where the input isn't exact? How do you calculate a "prediction interval" that accounts for both the uncertainty in the model and the uncertainty in the input data?
For example,
Assume the model is: y = a0 + (a1 * x1) + (a2 * x2)
Where y
, x1
, and x2
are time series vectors.
I only have one observation of the y
time series with its associated x1
and x2
time series. That data is used to fit the model. However, I also have 1000 additional observations for x1
and x2
. I can easily calculate individual "prediction intervals", using the model and each x1
x2
pair, however I want to estimate a "prediction interval" that allows for all x1
and x2
observations. To be more specific, I generated the example below. The questionable code is below the #======
comment line.
Basically, I used lm(..)
to fit the y
vector to its associated x1
and x2
vectors. Next, I used predict.lm(...interval="prediction", level=alpha)
to generate its typical fit
, lwr
, and upr
vectors for each additional x1
x2
pair (see Graph 1, fit
is the solid line, upr
and lwr
are the dashed lines for two observations). I then collected these vectors for all x1
x2
pairs and used alpha
to extract the upper and lower tails of all of the lwr
and upr
vectors (see Graph 2).
Is this scheme valid? Does the alpha
that was used in predict.lm(... level=alpha)
apply directly to counting-up/sorting the results to generate a "prediction interval" that allows for all x1
x2
pairs? Can a 5% range from a competing model (for example an ARIMA model) be compared to this 5% "prediction interval"?
I'm fairly sure that the following scheme isn't right, but so far, I haven't figured out what I need to fix.
set.seed(1)
numpoi <- 10 #Number of data points in a time series vector
#First independent variable "x1", first observation
x1mea <- 0.03
x1sta <- 0.05
x1 <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))
#Second independent variable "x2", first observation
x2mea <- -0.01
x2sta <- 0.1
x2 <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))
#Dependent variable "y", first observation
a0 <- 3
a1 <- 2
a2 <- 1
noimea <- 0.0
noista <- 0.1
y <- a0 + (a1 * x1) + (a2 * x2) + rnorm(numpoi, mean=noimea, sd=noista)
#Build a data frame of the "first observation" data
datfra <- data.frame(y=y, x1=x1, x2=x2)
#Fit the model for the "first observation"
mod <- lm(y ~ x1 + x2, data=datfra)
summary(mod)
#Set up desired alpha value
alpha <- 0.95
onetai <- (1 - alpha)/2 #Convert the two tail "alpha" to a one tail value for use later
#Generate some new data "a" for a second observation of "x1" and "x2".
x1a <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))
x2a <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))
datfraa <- data.frame(y=rep(NA, numpoi), x1=x1a, x2=x2a)
modprea <- predict(mod, newdata=datfraa, interval="prediction", level=alpha)
#Generate some new data "b" for a third observation of "x1" and "x2".
x1b <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))
x2b <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))
datfrab <- data.frame(y=rep(NA, numpoi), x1=x1b, x2=x2b)
modpreb <- predict(mod, newdata=datfrab, interval="prediction", level=alpha)
#Plot the results for both new data "a" and "b"
plot(modprea[, 1], type="l", ylim=c(min(modprea, modpreb), max(modprea, modpreb)), main="Graph 1 - Second and Third Observations for x1 and x2", lwd=2, col="red")
lines(modprea[, 2], lwd=2, lty=2, col="red")
lines(modprea[, 3], lwd=2, lty=2, col="red")
lines(modpreb[, 1], lwd=2, col="green")
lines(modpreb[, 2], lwd=2, lty=2, col="green")
lines(modpreb[, 3], lwd=2, lty=2, col="green")
#===========================================================================
#The code below is where my question lies. Is this the appropriate method
#to account for all observations?
#Run the above calculation scheme on "all" observations.
numtri <- 1000 #All observations
modprecfit <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "fit" vectors
modpreclwr <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "lwr" vectors
modprecupr <- matrix(0, nrow=numpoi, ncol = numtri) #Matrix to hold all "upr" vectors
#Fill up the modprecfit, modpreclwr, and modprecupr matricies.
for (i in 1:numtri) {
#Generate the new data and put it in a dataframe
x1c <- cumsum(rnorm(numpoi, mean=x1mea, sd=x1sta))
x2c <- cumsum(rnorm(numpoi, mean=x2mea, sd=x2sta))
datfrac <- data.frame(y=rep(NA, numpoi), x1=x1c, x2=x2c)
#Predict the new "y" for new input data
modprec <- predict(mod, newdata=datfrac, interval="prediction", level=alpha)
#Store the "ft", "lwr", and "upr" vectors so they can be processed later
modprecfit[, i] <- modprec[, 1]
modpreclwr[, i] <- modprec[, 2]
modprecupr[, i] <- modprec[, 3]
}
#Extract the average "fit", the lower quantile "lwr", and the upper quantile "upr"
modprecfitfin <- apply(modprecfit, 1, quantile, 0.5)
modpreclwrfin <- apply(modpreclwr, 1, quantile, onetai)
modprecuprfin <- apply(modprecupr, 1, quantile, (1 - onetai))
plot(modprecfitfin, type="l", ylim=c(min(modpreclwrfin), max(modprecuprfin)), main="Graph 2 - All Observations for x1 and x2", lwd=2, col="red")
lines(modpreclwrfin, lwd=2, lty=2, col="red")
lines(modprecuprfin, lwd=2, lty=2, col="red")