-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathXAMPLE2.R
117 lines (117 loc) · 3.28 KB
/
XAMPLE2.R
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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
setwd("c:/grant/backup/books/lcurve/rprogram/examples")
####################
# load and plot data
####################
co2data = read.table("mloco2.dat")
t = co2data[,1]
x = co2data[,2]
plot(t,x,pch=".",cex=3,
xlab="Time (years)",
ylab="CO2 concentration (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
#########################################
# simplest model: linear function of time
#########################################
xfit = lm(x ~ t)
summary(xfit)
##########################
# superimpose linear model
##########################
lines(t,xfit$fitted.values,col="red",lwd=3)
################
# plot residuals
################
plot(t,xfit$residuals,pch=".",cex=3,
xlab="Time (years)",
ylab="Residuals (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
##############################
# quadratic model of residuals
##############################
tsquared = t^2
xfit2 = lm(xfit$res ~ t+tsquared)
summary(xfit2)
###########################
# superimpose model on data
###########################
lines(t,xfit2$fit,col="red",lwd=3)
#############################
# quadratic model of raw data
#############################
xfit3 = lm(x ~ t+tsquared)
summary(xfit3)
#################################
# plot data and superimpose model
#################################
plot(t,x,pch=".",cex=3,
xlab="Time (years)",
ylab="CO2 concentration (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
lines(t,xfit3$fit,col="red",lwd=3)
################
# plot residuals
################
plot(t,xfit3$resid,pch=".",cex=3,
xlab="Time (years)",
ylab="Residuals (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
######################
# residuals since 2000
######################
plot(t,xfit3$resid,type="o",pch=".",cex=3,
xlim=c(2000,2010),
xlab="Time (years)",
ylab="Residuals (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
############################
# quadratic + sinusoid model
############################
c1 = cos(2*pi*t)
s1 = sin(2*pi*t)
xfit4 = lm(x ~ t+tsquared+c1+s1)
plot(t,x,pch=".",cex=3,
xlab="Time (years)",
ylab="CO2 concentration (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
lines(t,xfit4$fit,col="red",lwd=2)
####################
# closeup since 2000
####################
plot(t,x,pch=".",cex=3,
xlim=c(2000,2010),
ylim=c(365,390),
xlab="Time (years)",
ylab="CO2 concentration (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
lines(t,xfit4$fit,col="red",lwd=2)
##################################
# quadratic + sinusoid w/harmonics
##################################
c2 = cos(4*pi*t)
s2 = sin(4*pi*t)
c3 = cos(6*pi*t)
s3 = sin(6*pi*t)
xfit5 = lm(x ~ t+tsquared+c1+s1+c2+s2+c3+s3)
plot(t,x,pch=".",cex=3,
xlim=c(2000,2010),
ylim=c(365,390),
xlab="Time (years)",
ylab="CO2 concentration (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)
lines(t,xfit5$fit,col="red",lwd=2)
###########
# residuals
###########
plot(t,xfit5$resid,type="o",pch=".",cex=3,
xlab="Time (years)",
ylab="Residuals (ppmv)",
main="Mauna Loa",
cex.axis=1.5,cex.lab=1.5,cex.main=2)