Generating and graphing points to which we will fit the model. Signal is shown through the red curve.
x <- seq(0.1,15,0.1)
y <- -0.72*(x-7.5)^3 + 20
noise <- rnorm(length(x), 7.5, 55)
y.noise <- y + noise
plot(x,y.noise, ylab = "y")
lines(x,y, lwd = 3, col = "red")
Now we will fit a polynomial model to the data using lm and then approximate our bounds using confint to perform a confidence interval. Residual plot included along with final plot showing signal, bounds, and polynomial regression curve.
model <- lm(y.noise ~ poly(x,3))
confint(model, level=0.95)
## 2.5 % 97.5 %
## (Intercept) 15.52304 33.12152
## poly(x, 3)1 -1311.93139 -1096.39491
## poly(x, 3)2 -99.01748 116.51900
## poly(x, 3)3 -743.97899 -528.44250
plot(fitted(model),residuals(model), xlab = "x", ylab = "residual")
Below is the code showing the plot with the signal and predicted curve included. Note that we performed polynomial regression of degree three to prevent overfitting, as one should be cautious of performing polynomial regression of degree four or higher.
predicted.intervals <- predict(model,data.frame(x),interval='confidence',
level=0.95)
plot(x,y.noise, ylab = "y")
lines(x,predicted.intervals[,2], col = "black")
lines(x,predicted.intervals[,3], col = "black")
lines(x,predicted.intervals[,1],lwd = 2, col = "green")
lines(x,y, lwd = 2, col = "red")
legend("bottomleft",c("Signal","Predicted"),
col=c("red","green"), lwd=3)