Week01
Shannon
April 2, 2016
Summary of Combined Land-Surface Air and Sea-Surface Water Temperature Anomalies (Land-Ocean Temperature Index, LOTI)
library(moments)
setwd('/Users/shannoncall/Desktop/ESCI597A/Week01')
load("/Users/shannoncall/Desktop/ESCI597A/Week01/climate_Time_Series_Extravaganza.Rdata")
summary(loti.ts)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.77000 -0.22000 -0.05000 0.02167 0.21750 1.35000Plot Dataset
plot(loti.ts)
####Linear Model
tsp(loti.ts) #Verify the structure of the data set using the time-series properties command. This data set starts in 1880, ends in 2016, and has a frequency of 12. This means there are 12 observations per year (Note: 2016 is not over, I may restrict the dataset to 1880-2015)## [1] 1880.000 2016.083 12.000kt <- kurtosis(loti.ts)
sk<- skewness(loti.ts)
my.xlim <- range(loti.ts)
h<- hist(loti.ts, breaks=10, col="magenta", xlab="Loti Index", main="",xlim=my.xlim)
xfit<-seq(min(loti.ts),max(loti.ts),length=100)
yfit<-dnorm(xfit,mean=mean(loti.ts),sd=sd(loti.ts))
yfit <- yfit*diff(h$mids[1:2])*length(loti.ts)
lines(xfit, yfit, col="gold2", lwd=2)
boxplot(loti.ts, horizontal=TRUE, outline=TRUE, axes=FALSE, ylim=my.xlim, col = "gold2", add=TRUE, lwd=3, boxwex=3)
text(x = 0.75, y=450, labels = paste("Kurtosis=", round(kt,2)),pos=4)
text(x = 0.75, y=400, labels = paste("Skewness=", round(sk,2)),pos=4)
Linear Model
LotiTime <- time(loti.ts) #Create an object called *LotiTime* and make the loti.ts a time series. Though, in hindsight it was already a time series. :)
Loti.lm <- lm(loti.ts ~ LotiTime) #Run a linear model of the LotiTime object
summary(Loti.lm)##
## Call:
## lm(formula = loti.ts ~ LotiTime)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57565 -0.14028 -0.00237 0.13641 0.85608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.350e+01 2.351e-01 -57.41 <2e-16 ***
## LotiTime 6.940e-03 1.207e-04 57.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1918 on 1632 degrees of freedom
## Multiple R-squared: 0.6696, Adjusted R-squared: 0.6694
## F-statistic: 3307 on 1 and 1632 DF, p-value: < 2.2e-16This linear model has an adjusted R-squared value of 0.67, indicating the Loti index is changing over time. However, the residuals on this model are not identical and independently distributed. The coefficients are accurate though. Looks like the land-surface air and sea-surface water temperature anomalies have decreased by 0.007 degrees C over time.
Add some filters,
ma20 <- filter(x=loti.ts, filter=rep(x=1/20,times=20), sides=2)
ma10 <- filter(x=loti.ts, filter=rep(x=1/10,times=10), sides=2)Plot those filters
plot(loti.ts,col="dodgerblue4")
lines(ma20,col="firebrick4",lwd=2)
lines(ma10,col="goldenrod3",lwd=1)
abline(Loti.lm, col="darkslateblue",lwd=2, lty="dashed") #Trend line from the linear model.
head(cbind(loti.ts,ma20,ma10), n=10)## loti.ts ma20 ma10
## [1,] -0.29 NA NA
## [2,] -0.20 NA NA
## [3,] -0.18 NA NA
## [4,] -0.27 NA NA
## [5,] -0.14 NA -0.196
## [6,] -0.28 NA -0.185
## [7,] -0.22 NA -0.185
## [8,] -0.07 NA -0.176
## [9,] -0.16 NA -0.162
## [10,] -0.15 -0.1475 -0.146Decomposition and Seasonality - I’m betting there is some seasonality in this index.
f <- frequency(loti.ts)
filt <- c(0.5, rep(1, times=f-1), 0.5)/f
m.hat <- filter(x=loti.ts, filter=filt) #This creates a filter where we create a centered moving average
s.hat <- loti.ts - m.hat
par(mfcol=c(2,1))
plot(s.hat,ylab=expression("Temperature C"))
boxplot(s.hat~cycle(s.hat), xlab="Year", ylab=expression("Temperature C"))
These plots show the seasonality over time so that we might be able to detect the trend in the data without seasonal variation confounding our data.
Decomposition - sans belabored code
loti.ts.d <- decompose(loti.ts)
plot(loti.ts.d) 
Second: Earthquakes Per Year
read.csv("eqperyear.csv",header=T)## Year Number.of.earthquakes.per.year.magnitude.7.0.or.greater..1900.1998
## 1 1900 13
## 2 1901 14
## 3 1902 8
## 4 1903 10
## 5 1904 16
## 6 1905 26
## 7 1906 32
## 8 1907 27
## 9 1908 18
## 10 1909 32
## 11 1910 36
## 12 1911 24
## 13 1912 22
## 14 1913 23
## 15 1914 22
## 16 1915 18
## 17 1916 25
## 18 1917 21
## 19 1918 21
## 20 1919 14
## 21 1920 8
## 22 1921 11
## 23 1922 14
## 24 1923 23
## 25 1924 18
## 26 1925 17
## 27 1926 19
## 28 1927 20
## 29 1928 22
## 30 1929 19
## 31 1930 13
## 32 1931 26
## 33 1932 13
## 34 1933 14
## 35 1934 22
## 36 1935 24
## 37 1936 21
## 38 1937 22
## 39 1938 26
## 40 1939 21
## 41 1940 23
## 42 1941 24
## 43 1942 27
## 44 1943 41
## 45 1944 31
## 46 1945 27
## 47 1946 35
## 48 1947 26
## 49 1948 28
## 50 1949 36
## 51 1950 39
## 52 1951 21
## 53 1952 17
## 54 1953 22
## 55 1954 17
## 56 1955 19
## 57 1956 15
## 58 1957 34
## 59 1958 10
## 60 1959 15
## 61 1960 22
## 62 1961 18
## 63 1962 15
## 64 1963 20
## 65 1964 15
## 66 1965 22
## 67 1966 19
## 68 1967 16
## 69 1968 30
## 70 1969 27
## 71 1970 29
## 72 1971 23
## 73 1972 20
## 74 1973 16
## 75 1974 21
## 76 1975 21
## 77 1976 25
## 78 1977 16
## 79 1978 18
## 80 1979 15
## 81 1980 18
## 82 1981 14
## 83 1982 10
## 84 1983 15
## 85 1984 8
## 86 1985 15
## 87 1986 6
## 88 1987 11
## 89 1988 8
## 90 1989 7
## 91 1990 13
## 92 1991 10
## 93 1992 23
## 94 1993 16
## 95 1994 15
## 96 1995 25
## 97 1996 22
## 98 1997 20
## 99 1998 16eq<-read.csv("eqperyear.csv",header=T)
numb<-eq$Number.of.earthquakes.per.year.magnitude.7.0.or.greater..1900.1998
earthquakes<-ts(numb,start=c(1900,2),end=c(1998,2),frequency=2)Plot Earthquakes
plot(earthquakes,ylab="Earthquakes per year")
Figure 7. Earthquakes per year, since 1900 thru 1998
my.xlim <- range(earthquakes)
h<-hist(earthquakes, breaks=5, col="deeppink4", xlab="Number of Earthquakes",
main="",xlim=my.xlim)
xfit<-seq(min(earthquakes),max(earthquakes),length=8)
yfit<-dnorm(xfit,mean=mean(earthquakes),sd=sd(earthquakes))
yfit <- yfit*diff(h$mids[1:2])*length(earthquakes)
lines(xfit, yfit, col="darkturquoise", lwd=2)
boxplot(earthquakes, horizontal=TRUE, outline=TRUE, axes=FALSE,
ylim=my.xlim, col = "peachpuff2", add = TRUE, boxwex=3)
Figure 8. Histogram of the number of earthquakes per year.
Linear Model
eqtime <- time(earthquakes)
earthquakes.lm <- lm(earthquakes ~ eqtime)
summary(earthquakes.lm)##
## Call:
## lm(formula = earthquakes ~ eqtime)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.189 -5.090 -0.053 3.945 21.504
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.27796 35.37817 1.930 0.0551 .
## eqtime -0.02474 0.01815 -1.364 0.1743
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.242 on 195 degrees of freedom
## Multiple R-squared: 0.009446, Adjusted R-squared: 0.004366
## F-statistic: 1.859 on 1 and 195 DF, p-value: 0.1743Looks like there is a downward trend of -0.02474 earthquakes over time. Weird, not sure what that means.
Add Moving Averages and Trend Lines
ma3 <- filter(x=earthquakes, filter=rep(x=1/4,times=4), sides=2)
ma12 <- filter(x=earthquakes, filter=rep(x=1/10,times=10), sides=2)plot(earthquakes,col="steelblue4")
lines(ma3,col="lightsalmon3",lwd=2)
lines(ma12,col="olivedrab3",lwd=2)
abline(earthquakes.lm, col="magenta2",lwd=2, lty="dashed")
Figure 9. Filter with a four and ten year centered moving average.
A trend isn’t very clear. Let’s decompose to see what happens…
Decomposition
earthquakes.d<-decompose(earthquakes)
plot(earthquakes.d)