--- title: Testing a linear regression specification against non-parametric alternatives date: 16 February 2016 --- # First case: linear specification is wrong Make up sample data; it'll simplify things later if they're sorted for plotting ```{r} x <- sort(runif(300, 0, 3)) ``` Impose true regression function \$\log{(x+1)}\$, with Gaussian noise: ```{r} yg <- log(x+1) + rnorm(length(x), 0, 0.15) ``` Bind into a data frame: ```{r} gframe <- data.frame(x=x, y=yg) ``` Plot it, plus the true regression curve ```{r} plot(y~x, data=gframe, xlab="x", ylab="y", pch=16, cex=0.5) curve(log(1+x), col="grey", add=TRUE, lwd=4) ``` Fit the linear model, add to the plot: ```{r} glinfit <- lm(y~x, data=gframe) print(summary(glinfit), signif.stars=FALSE, digits=2) plot(x, yg, xlab="x", ylab="y") curve(log(1+x), col="grey", add=TRUE, lwd=4) abline(glinfit, lwd=4) ``` MSE of linear model, in-sample: `r signif(mean(residuals(glinfit)^2),3)`. We'll need to do that a lot, so make it a function: ```{r} mse.residuals <- function(model) { mean(residuals(model)^2) } ``` - EXERCISE: Write comments giving the inputs and outputs of the function - EXERCISE: Check that the function works Fit the non-parametric alternative: ```{r, message=FALSE} library(mgcv) ``` We'll use spline smoothing as provided by the `gam` function: ```{r} gnpr <- gam(y~s(x),data=gframe) ``` Add the fitted values from the spline to the plot: ```{r} plot(x,yg,xlab="x",ylab="y") curve(log(1+x),col="grey",add=TRUE,lwd=4) abline(glinfit, lwd=4) lines(x,fitted(gnpr),col="blue",lwd=4) ``` Calculate the difference in MSEs: ```{r} t.hat <- mse.residuals(glinfit) - mse.residuals(gnpr) ``` - "t" for test, NOT the t statistic for a Gaussian mean - EXERCISE: Why do we believe this is the right number? Simulate from the parametric model, assuming Gaussian noise ```{r} sim.lm <- function(linfit, test.x) { n <- length(test.x) sim.frame <- data.frame(x=test.x) sigma <- sqrt(mse.residuals(linfit)) # There are other ways to get sigma y.sim <- predict(linfit,newdata=sim.frame) y.sim <- y.sim + rnorm(n,0,sigma) sim.frame <- data.frame(sim.frame,y=y.sim) return(sim.frame) } ``` - EXERCISE: Write comments describing input and output of this function - EXERCISE: How can we check that this is working? - EXERCISE: Non-Gaussian noise? Calculate difference in MSEs between parametric and nonparametric models on a data frame: ```{r} calc.T <- function(df) { MSE.p <- mse.residuals(lm(y~x,data=df)) MSE.np <- mse.residuals(gam(y~s(x),data=df)) return(MSE.p - MSE.np) } ``` - EXERCISE: Comments once again (with feeling this time) - EXERCISE: Check this function Calculate the MSE difference on one simulation run: ```{r} calc.T(sim.lm(glinfit,x)) ``` Calculate the MSE difference on 200 simulation runs, so we get a sample from the null hypothesis: ```{r} null.samples.T <- replicate(200,calc.T(sim.lm(glinfit,x))) ``` How often does the simulation produce gaps bigger than what we really saw? ```{r} sum(null.samples.T > t.hat) ``` Plot histogram of the sampling distribution, and the observed value: ```{r} hist(null.samples.T,n=31,xlim=c(min(null.samples.T),1.1*t.hat),probability=TRUE) abline(v=t.hat, col="red") ``` Another approach: smooth the residuals of the parametric model: ```{r} plot(gam(residuals(glinfit) ~ s(x)),residuals=TRUE, lwd=4,ylab="Residuals of the linear model") abline(h=0,col="grey") ``` - We'd need to calculate how far that curve is from 0, and then simulate to get a null distribution for that test statistic # Second case: linear model is properly specified Deliberately left uncommented so that you can explore ```{r} y2 <- 0.2+0.5*x + rnorm(length(x),0,0.15) # Why don't I have to make up x values again? y2.frame <- data.frame(x=x,y=y2) plot(x,y2,xlab="x",ylab="y",pch=16) abline(0.2,0.5,col="grey",lwd=4) ``` ```{r} y2.fit <- lm(y~x,data=y2.frame) null.samples.T.y2 <- replicate(200,calc.T(sim.lm(y2.fit,x))) t.hat2 <- calc.T(y2.frame) hist(null.samples.T.y2,n=31,probability=TRUE) abline(v=t.hat2, col="red") ```