# Purpose: to demonstrate polynomial regression # # # Created by Hui Lan on 19 April 2018 ### Functions ### f <- function(x) { # the true function # y <- 0.5 + 1.5*x + 0.5*x^2 # also try quadratic function y <- 0.5 + 1.5*x # y is implicitly returned. don't need a return statement } ### main ### sd.err <- 2 N <- 12 x1 <- rnorm(N, mean=5, sd=3) # my data for the the predictor variable X1 e = rnorm(N, mean=0, sd=sd.err) # error term y <- f(x1) + e # my data for the response variable Y model <- lm(y~x1) # call linear model fit. use ?lm to get details of this function model summary(model) D <- data.frame(Y=y, X1=x1, X2=x1^2, X3=x1^3, X4=x1^4, X5=x1^5, X6=x1^6,X7=x1^7) model.polynomial <- lm(Y ~ X1+X2+X3+X4+X5+X6+X7, data=D) par(mfrow=c(1,1)) plot(x1, y, pch='+') abline(model, lwd=1, col='green') curve(f, from=min(x1), to=max(x1), col='red', add=T) num.interval <- 100 new.x <- seq(min(x1),max(x1),length=num.interval) test.x <- data.frame(X1=new.x, X2=new.x^2, X3=new.x^3, X4=new.x^4, X5=new.x^5, X6=new.x^6, X7=new.x^7) lines(new.x,predict(model.polynomial, newdata=test.x), col="blue")