1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
|
# 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")
|
