An Introduction to R for Geoscientists

A short documented course to get you started in R.

[r-module-geo-blog]

Introduction: Why Choose R?

As a young geoscientist, you will sooner or later face the difficult question: “Which software do you use for your data processing and analysis?”. Conventionally scientists and geoscientists have been forced to use the software provided by their department or research group [1]. This may be an in-house software package specially developed for the scientific tasks faced on a daily basis, or additionally a complicated and complex proprietary software offering.

Often getting started in these programs can be a virtual nightmare - help files are full of jargon, and often virtually no tutorials exist. Perhaps at best as a 'newbie' you may be provided with some scripts by a supervisor. Regardless starting in an unknown software package is always daunting and can be simply frustrating!!

This online module offers the possibility of using an open source, community developed software that is built by scientists for scientists. R is a fully developed, and tested suite of tools, functions and packages to deal with virtually any scientific task you can expect to face.

Developed in 1993, at the University of Auckland, New Zealand, The R Project For Statistical Computing, or simply R as it is commonly know has grown into one of the most commonly used and popular statistical environments available. The active community of users have themselves contributed over 5000 additional packages to the base installation [2].

Across the diverse fields of mathematics, quantitative finance, biology, astrophysics, and statistics - one can access and exploit specifically designed tools to work with many types of data.

Starting Advice

This online course, or condensed manual is meant to serve as a quick 'point of referral' for your own continued exploration of the R environment. It does not aim to overwhelm the new user with technicalities and unnecessary jargon but instead provide the 'nuts and bolts" to navigate R, and complete common tasks.

It is structured to progressively teach basic skills first, and then expand developing the techniques needed to successfully complete more lengthy tasks. Thus it is recommended to read and follow through the examples sequentially, as each ensuing chapter builds on the fundamentals presented in earlier sections.

[1]: For specific tasks such as processing of geophysical data, or in cases where packages were developed by industry (i.e. Oil and Gas) this may be unavoidable. For any other tasks R is well-suited.

[2]: The full range of packages can be found here:

[r-module-geo-blog]

Chapter 3 "Beginnings": The Interface

The R interface, regardless of which platform it is run on consists of two main parts. The first is the Console - where commands are input and run, and error messages displayed. This serves as the main working space for an R session. Secondly there is the main Menu Bar - here one can access standard File, Edit and Help commands.

Under the File Menu you can save and open a new script, or alternatively save your current R history - all of the commands and objects you hvae currently worked with, or load a saved file. Importantly you can also change your Working Directory via the File Menu, to navigate to where your data or scripts are stored. One can only directly load data from within the current working directory.

The Console - First Uses as A Calculator

Now lets tap into some of R's capabilities. Lets crunch some numbers…. :)

Anything we type into R is known as an expression. Each expression below is commented to provide additional detail or information about that line of code. It may explain the function, arguments or provide special notes to consider.

When working in any scientific endevaour it is good practice to consistently comment in your workflow. This not only helps remind youself what you were doing, but also acts is an valuable aid to any one whom may assess or run your code later.

25 * 2

## [1] 50

# create a sequence of numbers, assign to the variable 'x' x <- 1:10 # print the variable 'x' x

## [1] 1 2 3 4 5 6 7 8 9 10

# create a sequeunce of numbers x <- seq(1, 10, by = 2) # print the variable 'x' x

## [1] 1 3 5 7 9

# concatenate a vector i.e. create a sequence X <- c(1, 2, 3, 4, 5) # notice X and x are two different variables - naming is case sensitve! X

## [1] 1 2 3 4 5

R acts as both a simple calculator as shown above, and additionally a more advanced manipulator of so-called objects. Lets learn more about the different types of object you will encounter.

Objects: vectors, matrics, lists and data.frames

Believe it or not but in our first example we have already encountered our first objects!! Objects are something (and anything!) that can be assigned to a variable (in ordianry language we can think of this as “giving the object a name”). For example by using the <- we assigned a simple sequence of numbers (the object) to the variable called 'x'.

This programming technique allows us to recall and display the object at will - by simply typing in the variable name. In our first examples we used the letters 'x' and 'X' as our variable names.

This is the fundamental nature or R - it is an object-orientated language. This means that our data manipulation is based on the fact that data we work with is represented in the form of objects which can be assigned and saved as to a variable. Later we can access and manipulate these data in our workspace via their respective variable names: performing calculations, running functions and much more.

There are several main object types. Lets introduce each in turn:

Vectors

Vectors underpin the majority of calculations in R. They are the simplest type of object - vectors can be numeric, character, and logical [TRUE, FALSE].

# a numeric vector c(100, 200, 300, 400, 500)

## [1] 100 200 300 400 500

# a character vector c("Red", "Green", "Blue")

## [1] "Red" "Green" "Blue"

# a logical vector c(T, F, F, T)

## [1] TRUE FALSE FALSE TRUE

Constructing vectors

Here are a few more ways to construct vectors using the concatenate function. As shown below we can combine several arguments within the parenthesis to produce more eleborate sequences:

# combining the standard concatenate function with seq c(1:3, seq(15, 20, length = 10))

## [1] 1.00 2.00 3.00 15.00 15.56 16.11 16.67 17.22 17.78 18.33 18.89 ## [12] 19.44 20.00

# replicate the vector a given no. of times rep(1:3, 3)

## [1] 1 2 3 1 2 3 1 2 3

# replicate combined with length rep(1:5, length = 18)

## [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3

Manipulating Vectors

Vectors can be intuitively manipulated with common mathematical expressions:

# multiplied or divided c(1, 10, 2, 20) * 5

## [1] 5 50 10 100

# addition of the elements of two vectors c(1, 2, 3) + c(4, 2, 4)

## [1] 5 4 7

# test a expression i.e >, <, >=, <=, ==. yields a LOGICAL vector c(25, 100, 50, 200) > 75

## [1] FALSE TRUE FALSE TRUE

We usually want to work with multiple vectors at the same time. So as mentioned earlier its essential to assign a suitable variable name to each vector when we create it…

# assign a vector to the variable one.to.five one.to.five <- c(1, 2, 3, 4, 5) # assign a sequence of numbers to the variable my.vector my.vector <- c(1:20)

We can also apply some pre-designated functions to these vectors as shown:

# calculate the sum sum(one.to.five)

## [1] 15

# calculate the mean mean(one.to.five)

## [1] 3

# calculate the median median(one.to.five)

## [1] 3

# calculate the std dev sd(one.to.five)

## [1] 1.581

Object management

Numeric vectors can contain integers, decimals as well as not-a-number (NaN), which is often used to denote missing data.

# create an 'incomplete' vector some.data <- c(1, 2, NaN, 4, NaN) # print contents of variable some.data

## [1] 1 2 NaN 4 NaN

We can keep track of which objects are currently in our workspace by using the ls() command:

# show the 4 objects, each is of the type vector ls()

## [1] "my.vector" "one.to.five" "some.data" "x" "X"

We can check the length of our objects

# show length of vector length(some.data)

## [1] 5

Using the str() we can find out the properties of any object

# a numeric vector str(one.to.five)

## num [1:5] 1 2 3 4 5

With some practice creating and manipulating vectors will become second-nature.

Matrices and Arrays

Matrics can be seen as vectors with additional dimensions.

Lets create our first matrix:

# create a vector x <- c(1:6) # designate the dimensions, creating a matrix dim(x) <- c(3, 2)

Lets inspect our matrix.

# notice the matrix is filled column-wise! x

## [,1] [,2] ## [1,] 1 4 ## [2,] 2 5 ## [3,] 3 6

# show no of rows in x nrow(x)

## [1] 3

# show no of columns in x ncol(x)

## [1] 2

We could now create a 3-D array using the same assignment technique as above;

Note: As geoscientists we often deal with 3D dimensional data and so it can be useful to simultaneously think in terms of x, y and z (or t) space when constructing arrays.

# array with 1 row, 2 columns and 3 'levels' dim(x) <- c(1, 2, 3) x

## , , 1 ## ## [,1] [,2] ## [1,] 1 2 ## ## , , 2 ## ## [,1] [,2] ## [1,] 3 4 ## ## , , 3 ## ## [,1] [,2] ## [1,] 5 6

[r-module-geo-blog]

Sources, Credit and Appendix of Commands body, td { font-family: sans-serif; background-color: white; font-size: 12px; margin: 8px; } tt, code, pre { font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace; } h1 { font-size:2.2em; } h2 { font-size:1.8em; } h3 { font-size:1.4em; } h4 { font-size:1.0em; } h5 { font-size:0.9em; } h6 { font-size:0.8em; } a:visited { color: rgb(50%, 0%, 50%); } pre { margin-top: 0; max-width: 95%; border: 1px solid #ccc; white-space: pre-wrap; } pre code { display: block; padding: 0.5em; } code.r, code.cpp { background-color: #F8F8F8; } table, td, th { border: none; } blockquote { color:#666666; margin:0; padding-left: 1em; border-left: 0.5em #EEE solid; } hr { height: 0px; border-bottom: none; border-top-width: thin; border-top-style: dotted; border-top-color: #999999; } @media print { * { background: transparent !important; color: black !important; filter:none !important; -ms-filter: none !important; } body { font-size:12pt; max-width:100%; } a, a:visited { text-decoration: underline; } hr { visibility: hidden; page-break-before: always; } pre, blockquote { padding-right: 1em; page-break-inside: avoid; } tr, img { page-break-inside: avoid; } img { max-width: 100% !important; } @page :left { margin: 15mm 20mm 15mm 10mm; } @page :right { margin: 15mm 10mm 15mm 20mm; } p, h2, h3 { orphans: 3; widows: 3; } h2, h3 { page-break-after: avoid; } }

Sources, Credit and Appendix of Commands

Many Thanks…

Without the reference and insight of several key contributors and sources this short guide would not be possible.

My sincere thanks go out to all those who have contributed by sharing their advice, insight and experiences!!

Important sources

Quick R

Revolutions

Using R

Rtips Revival 2012!

R Tutorials: Carolina Coastal University

A special thanks to Kelly O'Day at Climate Charts and Graphs - from which specific inspiration and large parts of the implementation of Chapter 7 were originally derived.

Appendix of Commands

Some commands used within this module are included here, in alphabetical order.

Each is listed with a brief description and the arguments required.

abline(h=15,col=“red”,lty=“3”) - add a horizontal line

apropos(“plot”) - list all functions with 'plot' in title

c(i,ii,iii) - concatenate a vector - i can be numeric, character or logical

cwd() - show the current working directory

dim(my.array) <-c(x,y) - assign my.array dimensions of x and y

getwd() - get current working directory

head(my.data) - show first few rows of object my.data

hist(new.data$temp) - plot a histogram of new.data, column 'temp'

history() - display the last 25 commands

length(my.vector) - show length of my.vector

load(file=“some.data.Rdata”) - load workspace objects

load.history(file=“myfile”) - load command history

ls() - list all objects in the workspace

mean(my.vector) - calculate mean of my.vector

nrow(x) - show no of rows in array x

ncol(X) - show no of cols in array x

numeric(x) - create numeric vector of length x

par(new=T) - set overplotting to on, use after plot() command

par(new=F) - set overplotting to off

plot(one,two) - a simple x-y plot of object one and two

points(mydata_x,mydata_y,col=“red”) - add points to a x-y plot

read.table(“mydata.csv”,sep=“,”,header=T) - load csv file, with header and sep “,”

rep(i:ii,iii) - replicate vector seq i:ii, iii times

save.history(file=“myfile”) - save command history

save.image() - save all workspace objects to file .Rdata

scan(mytextfile.ascii) - scan in .ascii file one item at a time

scan(mytextfile.txt) - scan in .txt file one item at a time

sd(my.vector) - calculate std dev of my.vector

seq(i,ii,by=iii) - create a vector from i to ii, by iii

str(my.object) - list attributes of object my.object

sum(my.vector) - calculate sum of my.vector

tail(my.data) - show last rows of object my.data

text(x,y,“my.message”) - add text to plot at pos x and y

[r-module-geo-blog]

Chapter 2: Obtaining R, and Installation

The R Project for Statistical Computing can obtained free, at the homepage. Under the Download/Packages tab one can find a list of mirror sites where the respective download files are stored.

Any one of these mirrors will have R available for Linux, Mac OS X and Windows as source code and pre-compiled binaries. For the average user the pre-compiled binaries are the easy and smart choice. For reference the R environment used in this guide is Version 2.13 (Windows).

Regardless of which operating system you choose, all features are common across the platforms with minor aesthetic differences. It shouldn't be a problem following these examples on any platform or version.

A Note For Windows Users

Windows users have the option of 32-or/and 64-bit versions depending on their operating system's capabilities. It is possible to have both versions installed at one time; generally the advantages of the 64-bit version is access to more memory (>4Gb) and better/faster handling of larger objects. If is recommended “using the 'native' build, that is the 32-bit version on 32-bit Windows and the 64-bit version on 64-bit Windows”.

A Note About Languages

Important Note: R is available in a range of languages, but in Windows it will by default try to run in the language you run Windows in. For example if your Windows runs in German, R will run by default in this language despite the fact you may have downloaded and installed the English, Italian or Hebrew version!!

To change this behaviour: right click on the desktop shortcut, created during installation and select 'Properties'. Then in the tab named Target: scrolll to the end of the given path (“C:\Program Files\R-2.13.2\bin\x64\Rgui.exe” and add ´LANGUAGE = en´ after this path.

Updating and Add-On Packages

One of the great advantages of R over similar statistical software environments is the ease at which the user-community can access the tools and functionality they desire, and require for their work.

As mentioned earlier, R has a a dedicated group of thousands of users who invent, produce and share new functions, routines, and complete add-on packages with the rest of the community. These packages provide essential functionality, for example new methods of statistical analysis, production of graphics, or simpler tools to convert units and measurements. More advanced packages provide a whole new GUI interface to the R environment, or the ability to access certain data types etc.

In this course, we will actively introduce, and use some of the more common packages. A brief introduction to installing any package is presented below:

The easiest way to install packages is within the program itself. Clicking on the menu Packages and then Install Package/s. This will open a window where the nearest server to your location can be chosen, from here an alphabetical list of the available packages on CRAN TaskView is displayed.

Scroll through to choose your desired package, and click OK. The package (usually <1000kb) will be downloaded and 'installed'.

Important Note: To access the package during a working session one needs to load it. Type require(packagename) in the command line or alternatively use the Packages => Load Package menu.

[r-module-geo-blog]

Chapter 4: Optimizing Your Productivity

Utilizing Your Workspace

Now that we are familiar with handling many of the object types in the R interface it is time to use and save scripts, and export our workspace.

Scripts

R is a programming language as we have seen, and despite the ease of typing out commands it is always better to save our commands as we type them… then the can be run later, a tremendous convenience and time-saver!

Scripts are simply a type of document that store any commands we type.. additionally we can comment within our scripts using the hash-tag character.

Using the File Menu, select File-> New -> R Script. Not surprisingly a new command line interface opens up. It is currently named Untitled. Lets give it a name (Test.Script for example). You may notice the extension for R scripts is the filename.R. Therefore, we could write our scripts in any Notepad/Wordpad type editor and save the file with the extension .R to create a working R script which can be executed from the R console.

For now lets create some dummy data to illustrate this..

# create a simple vector vector.1 <- (1:15) vector.1

## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

# create another simple vector time <- (1983:1997) # x vs y plot of data plot(time, vector.1)

You may notice that after inputting your commands nothing happens!! ;)

A script does not execute your commands after pressing Return or Enter. Instead we must choose to run our script via the Console. The shortcut key is Ctrl+R (to run each line successively!). Try it to generate our script now and show its content.

Lets save the script now… using Ctrl+S or File -> Save.

We can load it into any R session now using File -> Open Script.. or Ctrl+O.

The REAL use of scripts in being able to formulate a working analysis in a logical and progressive manner. Line by line we can build up a 'typical workflow' for any one type of task.

Once we have built the entire script successively we can REUSE; and REAPPLY it when needed. A well commented and well programmed script also allows others to learn, repeat, and reproduce your work. An essential part of any scientific investigation.

Note: It is important to note that with any script commands are executed regardless of who uses them! :) This may sound funny but a script is just a set of commands, like any computer program - it can only do what it has been told to!

For example if a script requires a data file (a .txt file of numbers for instance) and cannot locate that file, or it 'pointed' to the wrong directory you can expect error messages!! Therefore it is necessary to distribute your scripts bundled with the data that they may require. The one exception to this rule is data accessed over the Internet.

Saving Time and Energy

Navigating Through Folders

To find out which is the current directory and navigate to another we can use the following sets of commands:

# get current working directory getwd()

## [1] "C:/Users/Daniel/Desktop/Guide Markdown/Chp 4"

# use forward slash / even on Windows!! setwd("C:/Users/Daniel/Desktop/Guide Markdown/Chp 3")

Saving Command History

To keep track of your workspace history (without using scripts) you can alternatively save your command history to a file…

# display the last 25 commands history()

## Error in history : function not supported in RStudio

# defualt is '.Rhistory' savehistory(file = "myfile")

## Error in savehistory : function not supported in RStudio

loadhistory(file = "myfile")

## Error in loadhistory : function not supported in RStudio

Saving Workspace Data and Objects

It may be even more useful to save your workspace data or individual objects, for further later analysis…

# save the workspace to the file .Rdata save.image() # saves selected object/s to file save(vector.1, file = "some_data.Rdata") # load file, cwd assumed unless path is specified! load(file = "some_data.Rdata")

Getting Help

Anytime you are in R and get stuck there are a few instantaneous options to get help. They should be your first stop before trying online forums like www.stackoverflow.com.

# general advice help.start()

## If nothing happens, you should open ## 'http://127.0.0.1:20908/doc/html/index.html' yourself

# help about function plot help(plot) # list all functions containing string 'plot' apropos("plot")

## [1] ".__C__recordedplot" "assocplot" "barplot" ## [4] "barplot.default" "biplot" "boxplot" ## [7] "boxplot.default" "boxplot.matrix" "boxplot.stats" ## [10] "cdplot" "coplot" "fourfoldplot" ## [13] "hook_plot_custom" "hook_plot_html" "hook_plot_md" ## [16] "hook_plot_rst" "hook_plot_tex" "interaction.plot" ## [19] "lag.plot" "matplot" "monthplot" ## [22] "mosaicplot" "plot" "plot.default" ## [25] "plot.density" "plot.design" "plot.ecdf" ## [28] "plot.function" "plot.lm" "plot.mlm" ## [31] "plot.new" "plot.spec" "plot.spec.coherency" ## [34] "plot.spec.phase" "plot.stepfun" "plot.ts" ## [37] "plot.TukeyHSD" "plot.window" "plot.xy" ## [40] "preplot" "qqplot" "recordPlot" ## [43] "replayPlot" "savePlot" "screeplot" ## [46] "spineplot" "sunflowerplot" "termplot" ## [49] "ts.plot"

# show an example of function plot example(plot.ts)

## ## plt.ts> require(graphics) ## ## plt.ts> ## Multivariate ## plt.ts> z <- ts(matrix(rt(200 * 8, df = 3), 200, 8), ## plt.ts+ start = c(1961, 1), frequency = 12) ## ## plt.ts> plot(z, yax.flip = TRUE)

## ## plt.ts> plot(z, axes = FALSE, ann = FALSE, frame.plot = TRUE, ## plt.ts+ mar.multi = c(0,0,0,0), oma.multi = c(1,1,5,1))

## ## plt.ts> title("plot(ts(..), axes=FALSE, ann=FALSE, frame.plot=TRUE, mar..., oma...)") ## ## plt.ts> z <- window(z[,1:3], end = c(1969,12)) ## ## plt.ts> plot(z, type = "b") # multiple

## ## plt.ts> plot(z, plot.type="single", lty=1:3, col=4:2)

## ## plt.ts> ## A phase plot: ## plt.ts> plot(nhtemp, c(nhtemp[-1], NA), cex = .8, col="blue", ## plt.ts+ main = "Lag plot of New Haven temperatures")

## ## plt.ts> ## a clearer way to do this would be ## plt.ts> ## Not run: ## plt.ts> ##D plot(nhtemp, lag(nhtemp, 1), cex = .8, col="blue", ## plt.ts> ##D main = "Lag plot of New Haven temperatures") ## plt.ts> ## End(Not run) ## plt.ts> ## plt.ts> ## xy.lines and xy.labels are FALSE for large series: ## plt.ts> plot(lag(sunspots, 1), sunspots, pch = ".")

## ## plt.ts> SMI <- EuStockMarkets[, "SMI"] ## ## plt.ts> plot(lag(SMI, 1), SMI, pch = ".")

## ## plt.ts> plot(lag(SMI, 20), SMI, pch = ".", log = "xy", ## plt.ts+ main = "4 weeks lagged SMI stocks -- log scale", xy.lines= TRUE)

[r-module-geo-blog]

Chapter 5: Data Exploration

FINALLY! We can begin to work with real geoscientific data, and dive into the possibilities of data exploration and visualisation!

In this chapter we will begin our journey by loading in data from outside R - then we can use some standard functions and plots to explore and display our data.

Start by navigating to the required directory, as explained in Chapter 3.

Reading in .csv files

A paragraph here. The dataset we will we using is stored in the so-called comma-separated values (.csv) format. This format is commonly used to store and arrange categorical data. One can use the read.table function to load the file:

# load csv file using read.table function radio_data <- read.table("geob1105-14c-s.csv", sep = ",", header = TRUE)

Lets review our progress: we have loaded the dataset and assigned it to the object radio_data. Note the character string (the name of the file!) is enclosed in double quotes “ ”. The sep argument tells us the data is separated by a comma, and the header argument indicates our first line in the file contains the headers/titles of the in-file variables..

Lets find out what type of data we are working with:

# used to inspect the structure of an object in R str(radio_data)

## 'data.frame': 17 obs. of 10 variables: ## $ core : Factor w/ 2 levels "GeoB 1105-3",..: 1 1 1 2 2 2 2 2 2 2 ... ## $ Sample.ID : Factor w/ 17 levels "GrA-1215","KIA4095",..: 2 3 4 1 5 6 17 7 8 9 ... ## $ Bremen.ID : Factor w/ 8 levels "HB???","HB011",..: 1 1 1 2 1 1 8 1 1 1 ... ## $ species : Factor w/ 3 levels "G. crassaf.",..: 2 2 2 2 2 2 2 2 2 2 ... ## $ depth : int 1 10 19 23 63 83 103 123 163 183 ... ## $ composite.depth : num 0.01 0.1 0.19 0.38 0.78 0.98 1.18 1.38 1.78 1.98 ... ## $ conventional.radiocarbon.age: int 2460 4300 5980 9770 14670 17670 19320 21790 29310 31730 ... ## $ error.... : int 30 30 40 80 80 110 110 180 450 330 ... ## $ error.....1 : int 30 30 40 80 80 110 110 170 420 310 ... ## $ d13C : num 3.96 2.61 3.08 2.84 2.54 1.5 -0.67 2.27 2.35 -0.22 ...

Our variable radio_data contains all the data of the original .csv file. It is stored in a data.frame - containing ten variables (columns) and 17 observations (rows).

Additionally we can see we have some categorical data. The variables core, Sample.ID and Bremen.ID and species are each a Factor. This just means the data have been further sub-divided (classified) into various groups.

For example the variable species is a Factor with 3 levels - the levels corresponding to the 3 species of foraminfera classified.

We see that each level is identified by its corresponding name, and is additionally assigned a numeric id (i.e. 1, 2, 3 etc).

Lets inspect our data to get a better feel for the actual structure.

# show the first few rows of data head(radio_data)

## core Sample.ID Bremen.ID species depth composite.depth ## 1 GeoB 1105-3 KIA4095 HB??? G. sacculifer 1 0.01 ## 2 GeoB 1105-3 KIA4096 HB??? G. sacculifer 10 0.10 ## 3 GeoB 1105-3 KIA4097 HB??? G. sacculifer 19 0.19 ## 4 GeoB 1105-4 GrA-1215 HB011 G. sacculifer 23 0.38 ## 5 GeoB 1105-4 KIA4099 HB??? G. sacculifer 63 0.78 ## 6 GeoB 1105-4 KIA4100 HB??? G. sacculifer 83 0.98 ## conventional.radiocarbon.age error.... error.....1 d13C ## 1 2460 30 30 3.96 ## 2 4300 30 30 2.61 ## 3 5980 40 40 3.08 ## 4 9770 80 80 2.84 ## 5 14670 80 80 2.54 ## 6 17670 110 110 1.50

…and now….

# show the last few rows of data tail(radio_data)

## core Sample.ID Bremen.ID species depth composite.depth ## 12 GeoB 1105-4 KIA4107 HB??? G. sacculifer 243 2.58 ## 13 GeoB 1105-4 KIA6594 HB567 G. crassaf. 43 0.58 ## 14 GeoB 1105-4 KIA6595 HB568 N. dutertrei 43 0.58 ## 15 GeoB 1105-4 KIA7250 HB672 N. dutertrei 103 1.18 ## 16 GeoB 1105-4 KIA6592 HB565 G. crassaf. 103 1.18 ## 17 GeoB 1105-4 KIA6593 HB566 N. dutertrei 103 1.18 ## conventional.radiocarbon.age error.... error.....1 d13C ## 12 43360 2070 1650 2.79 ## 13 12970 70 70 -0.46 ## 14 13460 80 80 -0.29 ## 15 19460 110 110 -0.04 ## 16 19880 140 140 -1.28 ## 17 19570 140 130 -1.80

Histograms, Boxplots and Scatterplots

Lets start our data analysis using the approach of Summary Statistics. Lets investigate the distribution of our d13C data.

# select desired column using the $ notation hist(radio_data$d13C)

One line of code and we have some output! The only argument in the hist plotting function is designating which data to plot.

Admittedly this is not the most impressive plot ever - lets modify it for clarity…

# plot histogram, add custom title and label hist(radio_data$d13C, breaks = 10, main = "Isotopic carbon ratios of 3 foraminfera species", xlab = "d13C per mil", col = "red")

With a few lines we can dramatically change the appearance of any plot: the breaks argument divides our data into a number of “bins”. The main and xlab arguments assign labels and the col argument changes our plotting colour.

Lets go a bit further…

# designate hist breaks, scale axes and change color hist(radio_data$d13C, breaks = 10, font = 3, font.main = 4, main = "Isotopic carbon ratios of 3 foraminfera species", xlab = "d13C per mil", cex.main = 0.9, cex.lab = 0.85, cex.axis = 0.8, col = "259")

[r-module-geo-blog]

Chapter 6: Project Two - SOI and ENSO

In this chapter we attempt to replicate the scientific results generated by several leading climate institutes such as NCAR, NASA and others.

We will use the long term record of sea-level pressure (SLP) measurements to calculate the Southern Oscillation Index (SOI).

The SOI represents the differences in pressure difference between Tahiti and Darwin. The index is representative of the relative intensity and position of the trade winds in the central Pacific, associated atmospheric-oceanic coupling (upwelling strength and thermocline depth) and thus a measure of switches between the two ENSO states of El Nino (warm) and La Nina (cool).

Having a SOI at your disposal is useful as it can be quickly and easily correlated with other changes in the climatic system (i.e. rainfall patterns in Eastern Africa) and be used as a diagnostic tool for severe weather disaster prevention and preparation.

Handling Raw Data

This chapter focuses on some techniques to handle so-called “raw data”. Typically this is data which is not stored in a straightforward manner for direct analysis.

Firstly there is a preliminary phase of re-arranging and subsetting the data into a more usable form, before anaylsis begins.

Lets load in our required data, this time using a slightly different technique:

# read in txt file item by item Tahiti_raw <- scan("tahiti.ascii") Darwin_raw <- scan("darwin.ascii")

Here we have used the scan function - which can directly read in a .dat, .txt or in our case .ascii file. The data now exists as one long numeric vector, of length 1911 items.

Inspect our data (either in a simple text editor like Notepad++) or in R using the head()/tail() command.

We notice the data is arranged so that each line begins with the year of the the data (1866 etc.) and then each successive column is a monthly SLP value.

We require the ordered monthly values (to form one long vector) but the 'counters' of the years are a nuisance.

Subsetting using [] and expressions >, <, = etc.

We could manipulate the data into a matrix and extract only the monthly data but this requires a bit of effort and thought! Instead lets rather solve our problem by subsetting our data using a simple trick.

Notice the SLP values are only double-digits (12, 17 etc), and the years are recorded as 1989…i.e. four-digits.

Therefore lets choose to ONLY select the low values of our “raw data”…

# subset raw data into new object only selecting values of less than 100 # in this case Tahiti_slp <- (Tahiti_raw[Tahiti_raw < 100]) Darwin_slp <- (Darwin_raw[Darwin_raw < 100])

The [] or subset notation of square brackets is used to designate or subset a particular part of a object, whether it be a vector, array or dataframe.

In this case we use the < (less than) symbol to select the part/s of Tahiti_raw where the *values are less than 100.

Alternatively one could also subset data by using a vector selection i.e.[1:25] selecting the first 25 *positions of the object.

Lets check if it worked:

# show the length of an object length(Tahiti_slp)

## [1] 1764

# confirm we have a vector, now shorter str(Tahiti_slp)

## num [1:1764] 10.6 10.9 11.5 11.4 12.5 13.6 13.5 15.5 13.9 13.8 ...

Subsetting using vector expressions

You may have noted that the last few months of our data are missing or assigned NA. Lets just exclude these NA's in this exercise - we will only use the full record up until the end of 2011.

# select subset of data from index positions 1 through 1752 1866-2011, # excludes partially filled data from 2012 Tahiti_slp <- Tahiti_slp[1:1752] Darwin_slp <- Darwin_slp[1:1752]

Now we can get our first term in our Southern Oscillation Index (SOI) equation by a straight-forward subtraction of the SLP's from each other:

# vectors are the same length each element is subtracted from the # corresponding one Pdiff <- Tahiti_slp - Darwin_slp

Time Series Package

Now we need to calculate the average pressure differences between Tahiti and Darwin for each of the months of the year. In essence we need to calculate a climatology - the average differences over an extended period of time.

# freq = 12 for monthly data, start date of ts Pdiff.ts <- ts(Pdiff, frequency = 12)

Using the time.series package included in base R we have taken our vector and transformed it into a time.series object - use str() to check!

Using tapply and several other functions such as cycle and mean we can average our data for each successive month Jan, Feb, Mar, Apr….

# apply mean calculation iteratively tapply(Pdiff.ts, cycle(Pdiff.ts), mean)

## 1 2 3 4 5 6 7 8 9 10 11 12 ## 4.660 4.729 4.005 2.341 1.552 1.222 1.047 1.645 2.465 2.928 3.027 3.829

# aggregate average data into a data.frame Pdiff_avg.df <- aggregate(c(Pdiff.ts), list(month = cycle(Pdiff.ts)), mean)

# extract the data which lies in the un-named column 'x' Pdiff_avg <- Pdiff_avg.df$x

An Alternative Method: Using Matrices

An ALTERNATIVE METHOD in achieving the same outcome would have been to use matrix manipulation techniques:

# calculate length of data coverage in years no_years <- length(Tahiti_slp)/12 # create an empty array with 12 rows and 146 columns to test empty <- array(, dim = c(12, no_years)) # fill array of same dim with Pdiff data filled <- array(Pdiff, dim = c(12, no_years))

The matrix has now been filled with the Pdiff data (columnwise…) i.e the first column is filled with the first 12 months of data, and the second column the second 12 months and so on…

We can use the rowMean function to calculate the average for each required month; Jan, Feb, Mar (i.e each row) across the columns (the 146 year record)…

# rowMeans applied to filled matrix Pdiff_avg.mat <- rowMeans(matrix(filled, nrow = 12, ncol = 146)) # confirm results are the same! :D Pdiff_avg.mat

## [1] 4.660 4.729 4.005 2.341 1.552 1.222 1.047 1.645 2.465 2.928 3.027 ## [12] 3.829

Using tapply

Lets calculate the standard deviation of the Pdiff_avg.

# apply std dev calculation iteratively tapply(Pdiff.ts, cycle(Pdiff.ts), sd)

## 1 2 3 4 5 6 7 8 9 10 11 12 ## 2.223 2.134 2.056 1.470 1.296 1.514 1.622 1.617 1.791 1.676 1.579 1.925

# aggregate std dev data into a data.frame Pdiff_avgsd.df <- aggregate(c(Pdiff.ts), list(month = cycle(Pdiff.ts)), sd)

# extract the data which lies in the un-named column 'x' Pdiff_avgsd <- Pdiff_avgsd.df$x

Above the tapply family of functions were applied again to calculate the standard deviation, in an identical manner as that used used for the calcualtion of the mean pressures.

The final answers were then stored in an appropriately named vector for later use.

Loops: The Advantages of Vectors

The last part of our exercise would traditionally require a loop (and more thinking!). However we can take advantage of the fact that out data is stored in vectors, and these vectors are multiples of each other.

NB: Generally it is advisable to avoid loops when ever possible, as they can be complex to design and slow in execution!! In comparison vector calculations in R are orders of magnitude faster!!

Lets create an empty vector to store our calculation results.

# length of data to be stored - total months len <- no_years * 12 # vector to hold data: of type numeric and length, len Index <- numeric(len)

Now we can proceed with the actual calculation. This is where the power and usefulness of vectors becomes aparrent!! :D

R will evaluate ALL vectors in a continuous manner - we have already seen this before when we subtracted two vectors from each other! In that case each element in one vector was subtracted from the corresponding element in the other vector.

Furthermore shorter vectors will be continuously and successively evaluated ( i.e. 'looped') over and over within a longer vector.

Here is a simple example to illustrate this:

# a 4 element vector long_vec <- c(10, 20, 30, 40) # a 2 element vector short_vec <- c(2, 4) # divide the 'long' vector by the 'short' vector quick_calc <- long_vec/short_vec # show the results quick_calc

## [1] 5 5 15 10

We can apply this principle in our evaluation of the SOI. We need to calculate the index for each data-point in the time series - 1752 data points in total.

We remember our average and standard deviation data are stored in vectors of length 12 - each element representing one month of the year (J,F,M etc.). Within the ensuing calculation these vectors of length 12 will be successively evaluated over and over for each Pdiff value.

# SOI Index Equation Index[] <- 10 * (Pdiff - Pdiff_avg)/Pdiff_avgsd Index

## [1] -5.66711 -0.60328 -3.91821 -4.36159 1.91316 -4.76928 ## [7] -3.98943 5.29032 -2.03857 1.62220 10.59206 -2.22723 ## [13] 1.08165 0.33373 0.45979 5.16308 5.77118 -2.12672 ## [19] 3.40685 4.67192 0.75347 -2.55383 -7.13945 -8.98008 ## [25] -2.06777 -2.94581 -11.70132 1.08108 -13.51894 -14.01823 ## [31] -4.60579 -11.40659 -12.08993 -10.30932 -14.73866 4.00617 ## [37] -17.36497 -2.47730 -3.43177 17.40910 16.57365 11.74672 ## [43] 12.65220 9.61915 0.75347 9.97426 6.79245 2.96727 ## [49] 8.28033 2.67626 -4.89110 2.44175 -1.94487 -6.75119 ## [55] 0.94142 10.23755 -0.36335 -4.34356 -7.77271 -12.61623 ## [61] -10.61621 -1.07179 -15.10643 -3.00093 -3.48808 2.49776 ## [67] 2.17413 5.29032 5.22074 0.42904 1.09304 8.16177 ## [73] 23.12762 11.57787 14.56668 -6.40260 14.25883 18.35311 ## [79] 21.28119 10.23755 29.79071 18.92289 30.22337 24.78416 ## [85] 21.32794 10.64086 14.56668 14.00743 5.77118 -21.94591 ## [91] -12.00207 -3.36734 -20.46605 -14.48535 -2.70657 11.27846 ## [97] 7.83042 -2.94581 -2.94532 6.52375 13.48723 -12.03631 ## [103] -2.14036 9.00075 19.18095 21.90577 5.52591 10.75902 ## [109] 5.58083 -0.60328 -6.83688 9.24509 8.08600 18.35311 ## [115] -5.22214 -0.27532 7.45437 12.36056 -13.47213 -5.86338 ## [121] 9.63009 10.17235 -0.02665 7.88442 7.31439 13.72864 ## [127] -5.83850 12.09277 9.12960 -7.92302 -3.33984 -4.82448 ## [133] -8.36662 -13.72144 -7.80977 -26.81262 -13.51894 -28.55230 ## [139] -10.15300 -7.69617 -16.55719 -15.67850 -12.83886 -14.69403 ## [145] -9.26645 -21.68604 -14.13354 -7.08293 2.68476 -2.78736 ## [151] 14.50127 12.71117 15.83050 10.57083 13.75840 16.47296 ## [157] 10.97984 13.45189 11.64801 10.60576 2.68476 13.06800 ## [163] 20.04848 21.98723 16.94731 14.74686 8.69225 -7.42173 ## [169] 9.18017 6.89281 12.62090 4.48275 12.71562 7.12224 ## [175] 0.94142 13.94798 6.89596 4.60507 6.15918 -3.78558 ## [181] -7.91670 -6.22535 1.43268 0.40074 -3.48808 -4.10864 ## [187] -5.83850 -10.78819 -13.20674 -23.43398 6.15918 5.56452 ## [193] -7.46678 -1.54029 3.86490 1.08108 8.08600 -10.71503 ## [199] -20.01470 -25.01148 -14.32356 1.62220 1.72631 8.16177 ## [205] 5.13091 -8.56788 -23.37599 12.64676 14.25883 3.15840 ## [211] -10.76935 2.19830 -8.73947 4.00850 3.62611 -17.81072 ## [217] -12.41588 -6.69386 7.75646 -11.16493 1.91316 7.78288 ## [223] -3.37308 -5.22255 -5.94743 3.41192 -1.44003 -14.69403 ## [229] -16.01522 0.80224 4.83779 -0.27959 -1.94487 -11.37567 ## [235] -4.60579 -8.93297 -3.71380 -18.06480 -15.37193 4.00617 ## [241] -2.06777 0.33373 2.89201 2.44175 8.08600 3.15840 ## [247] 7.10499 13.94798 12.48005 13.55371 9.32552 12.83681 ## [253] 10.97984 9.70385 8.72935 7.20409 -3.48808 3.81904 ## [259] 3.40685 4.05351 4.66233 5.79823 -5.23964 4.00617 ## [265] -4.31736 -2.00880 -10.24199 -18.64861 -8.88931 -12.69695 ## [271] -15.08385 -7.69617 -9.85629 -13.29220 -13.47213 -4.30503 ## [277] -25.01357 -2.00880 -25.32177 -0.27959 -0.40166 18.35311 ## [283] 1.55778 1.57990 9.12960 3.41192 20.72435 20.62856 ## [289] 18.17852 10.64086 12.62090 7.20409 4.22797 5.80096 ## [295] -3.37308 -2.74893 7.45437 4.00850 1.09304 -1.70778 ## [301] 14.12926 -4.35133 -9.26910 4.48275 0.36995 -0.80544 ## [307] -6.45486 -7.07776 -9.29788 1.02562 -5.23964 -6.90228 ## [313] 1.08165 -10.91041 10.67512 5.84342 11.17242 15.04992 ## [319] 7.72134 5.90872 -0.36335 9.37768 -5.87291 1.40892 ## [325] 7.83042 7.36132 3.86490 9.92542 15.80204 15.04992 ## [331] 16.96669 9.00075 10.80482 6.99138 9.32552 -5.34393 ## [337] 17.72860 9.70385 7.27001 7.88442 -7.34610 1.83712 ## [343] 4.02320 -11.40659 -1.48016 3.41192 -7.13945 0.37002 ## [349] 3.78116 5.95580 -0.51310 8.56475 3.45637 -6.75119 ## [355] -7.68757 -8.93297 -4.83062 -6.13329 -5.87291 -9.49953 ## [361] 4.23107 8.29833 -0.99954 -0.95992 -31.26585 -19.30335 ## [367] -13.85114 -14.49861 -13.20674 -12.69562 -5.87291 -11.57733 ## [373] -11.96596 -7.63087 -15.10643 -14.56660 -17.37696 -0.80544 ## [379] -1.52401 2.19830 0.19506 2.21877 -8.40598 8.68122 ## [385] 5.58083 5.01879 16.99890 8.56475 -1.17326 -1.46608 ## [391] 5.87227 8.38234 8.57119 4.60507 2.35957 2.44782 ## [397] 10.97984 7.36132 11.64801 2.44175 -5.80289 -7.41183 ## [403] -5.83850 -8.93297 -2.03857 6.39480 14.39167 -4.82448 ## [409] -7.46678 -7.16236 -23.37599 -15.24694 -6.57450 20.33503 ## [415] 9.57041 7.76394 -15.44038 -17.46823 -6.50618 -6.90228 ## [421] -1.61785 2.20775 8.72935 3.80242 0.36995 15.71056 ## [427] 13.26855 10.23755 -14.88197 -21.64426 -9.67252 -3.26613 ## [433] 15.02910 -2.94581 10.18868 5.84342 8.85760 2.49776 ## [439] 0.32506 -7.69617 -16.55719 -7.32644 -4.60637 -4.30503 ## [445] -9.71637 -11.37892 15.05312 15.36809 9.62921 -0.14480 ## [451] 4.63956 0.96149 7.45437 5.20165 1.09304 14.39516 ## [457] 11.87967 14.85741 8.24290 25.57311 9.62921 -6.75119 ## [463] -9.53664 0.96149 0.75347 0.42904 -16.63847 0.37002 ## [469] -9.26645 -17.93800 -26.78110 -34.29629 -35.12388 -24.58847 ## [475] -19.39835 -7.69617 -6.50584 -6.13329 -17.90500 -15.73292 ## [481] -4.76728 -8.56788 -3.91821 -7.08293 -7.34610 -2.78736 ## [487] 5.25592 15.18479 16.38891 8.78110 19.45781 1.92837 ## [493] 4.68099 -2.00880 -0.99954 -5.72226 -4.25968 3.81904 ## [499] -0.29129 -7.69617 -0.36335 -3.15041 -3.97311 4.52562 ## [505] -8.36662 7.82982 -3.43177 8.56475 1.91316 -3.44800 ## [511] 12.03584 4.05351 16.94731 7.58795 2.35957 -6.38283 ## [517] -2.51769 -2.94581 -0.51310 -11.16493 2.68476 17.69247 ## [523] 10.18677 9.61915 0.19506 3.41192 6.79245 3.48672 ## [529] 4.23107 14.85741 11.16157 3.80242 1.91316 18.35311 ## [535] 19.43212 10.23755 13.03846 9.37768 18.19127 14.91461 ## [541] 1.53157 1.73925 2.89201 1.08108 -6.57450 -10.71503 ## [547] -12.61842 -10.78819 -8.18107 -12.09905 -8.40598 -3.78558 ## [553] -10.61621 -17.46949 -8.29621 -18.64861 -11.20413 -4.76928 ## [559] -1.52401 -7.69617 -4.27221 -6.72986 1.72631 -10.53843 ## [565] -3.86744 -5.28834 0.45979 -5.72226 -6.57450 -4.10864 ## [571] -2.14036 -7.69617 -9.85629 -9.11617 -11.57232 -8.98008 ## [577] -6.11703 1.27074 8.24290 -14.56660 -8.88931 -16.00015 ## [583] -11.38571 -13.26180 -13.76515 -7.92302 -5.23964 0.37002 ## [589] -25.01357 -5.28834 -19.48443 -12.52560 10.40081 19.67439 ## [595] 22.51390 13.32957 12.48005 9.37768 -14.73866 8.16177 ## [601] 3.78116 -3.41432 -6.35043 -0.27959 6.54279 7.78288 ## [607] 23.74662 15.18479 3.54551 6.39480 9.32552 13.35626 ## [613] 3.78116 10.17235 16.02601 18.08943 21.97488 17.69247 ## [619] 24.97933 33.11851 27.55708 15.34344 20.09108 21.14801 ## [625] 12.32959 15.79442 -2.45888 14.68776 11.17242 -3.44800 ## [631] -14.46749 -3.98574 -7.62266 -4.94014 -0.17350 -9.49953 ## [637] -15.56530 -11.84742 -11.70132 -1.64026 -7.34610 -8.73311 ## [643] -8.92028 -7.69617 -4.83062 -10.30932 -11.57232 -11.05788 ## [649] 0.63173 -2.94581 -4.40465 -0.27959 -1.94487 5.14032 ## [655] 8.33770 4.67192 5.22074 -4.34356 -0.17350 7.64232 ## [661] 10.08000 5.48729 10.18868 -7.76326 4.22797 27.60207 ## [667] 2.79049 -5.84095 3.54551 9.97426 6.79245 6.08397 ## [673] 6.03075 8.29833 4.83779 -3.68126 -5.03129 3.81904 ## [679] 0.94142 -0.89372 4.10392 5.20165 8.05899 9.72012 ## [685] 4.23107 3.14477 7.75646 7.20409 2.68476 -0.14480 ## [691] -11.38571 -17.59063 -13.76515 -5.53671 -12.83886 0.88947 ## [697] -5.66711 0.33373 2.89201 -11.84527 11.94402 5.80096 ## [703] 6.48863 10.85596 6.33755 7.58795 11.22533 3.48672 ## [709] 4.23107 13.92040 12.62090 11.28609 -0.40166 -3.44800 ## [715] -13.23478 -9.55138 -5.94743 -13.29220 -10.30579 -10.01898 ## [721] -6.56695 -15.12696 -11.70132 -5.72226 -1.17326 -5.42992 ## [727] -1.52401 -7.07776 1.31188 4.00850 0.45977 4.00617 ## [733] 3.78116 -0.13478 16.02601 7.20409 7.31439 7.78288 ## [739] 1.55778 -5.84095 -0.36335 -4.34356 -8.40598 6.08397 ## [745] -9.71637 9.70385 11.64801 10.60576 -1.94487 -6.09056 ## [751] 0.32506 9.61915 8.01278 9.97426 1.09304 10.75902 ## [757] 13.67935 17.66845 4.83779 3.80242 -11.20413 1.17648 ## [763] 2.17413 0.34309 -0.92176 7.58795 10.59206 4.52562 ## [769] 10.52992 5.95580 0.94623 -1.64026 1.91316 -4.76928 ## [775] -4.60579 -2.13053 -7.06425 4.60507 1.09304 -3.78558 ## [781] 6.03075 -14.65846 4.35135 7.88442 14.25883 13.72864 ## [787] 8.33770 -3.36734 2.98710 -7.92302 -3.33984 6.60342 ## [793] 1.53157 -2.94581 -5.37754 -18.64861 2.68476 -0.80544 ## [799] -15.08385 -6.45936 -9.85629 -4.34356 -4.60637 1.40892 ## [805] -11.51604 3.61327 -1.48599 3.12208 6.54279 -3.44800 ## [811] 2.17413 -0.89372 2.42869 3.41192 6.15918 6.60342 ## [817] 4.68099 0.33373 0.45979 5.16308 -5.80289 9.76480 ## [823] 2.17413 -21.91946 -6.50584 3.41192 11.22533 -3.78558 ## [829] 5.13091 -4.81983 11.16157 3.12208 -4.25968 -1.46608 ## [835] 0.32506 2.19830 5.22074 7.58795 2.35957 -5.86338 ## [841] -2.51769 -0.60328 0.94623 5.16308 5.77118 -2.12672 ## [847] 3.40685 -7.69617 1.31188 0.42904 -13.47213 -1.18833 ## [853] 7.38050 -5.75685 5.81068 1.08108 -0.40166 2.49776 ## [859] -5.22214 2.81670 0.19506 -3.15041 -3.97311 4.52562 ## [865] 5.58083 3.61327 -2.94532 3.80242 13.48723 15.04992 ## [871] 17.58305 11.47436 7.45437 12.95713 1.09304 12.31736 ## [877] 14.57918 6.42431 9.70223 7.88442 -0.40166 -0.80544 ## [883] 7.72134 -0.89372 -9.85629 -15.08192 -9.03925 -11.05788 ## [889] -0.71802 -5.28834 -9.26910 -7.08293 -13.51894 -16.66079 ## [895] -15.70021 -16.97223 -18.23242 -17.46823 -7.77271 -21.96632 ## [901] -10.16629 -16.06397 -9.75554 -9.12393 -5.80289 -12.03631 ## [907] -20.01470 -18.20904 -8.18107 -19.85453 -9.67252 -10.53843 ## [913] -13.31571 -4.35133 -5.37754 -4.36159 5.77118 6.46160 ## [919] -1.52401 4.05351 7.45437 8.18453 -4.60637 12.31736 ## [925] 7.83042 9.70385 3.37846 11.28609 3.45637 -6.75119 ## [931] 2.17413 7.76394 4.66233 8.78110 2.99284 -10.53843 ## [937] -8.81654 3.14477 4.83779 -4.36159 -0.40166 -3.44800 ## [943] -8.92028 3.43511 1.87029 -8.51959 -7.13945 2.44782 ## [949] 3.78116 5.48729 11.64801 -5.72226 0.36995 6.46160 ## [955] 2.79049 11.47436 7.45437 2.21877 -3.97311 5.04507 ## [961] -3.41753 3.61327 -1.97243 -7.76326 -10.43252 -8.07247 ## [967] -10.15300 -3.98574 -15.44038 -12.09905 -2.07330 -7.42173 ## [973] -5.66711 -4.81983 10.18868 -3.68126 -12.74733 1.83712 ## [979] 8.33770 7.14553 10.24641 -1.95726 8.05899 3.48672 ## [985] -3.86744 -3.41432 -3.91821 2.44175 4.22797 -4.10864 ## [991] 0.32506 -3.98574 -7.62266 5.79823 3.62611 -7.42173 ## [997] -7.91670 1.27074 4.83779 1.08108 -5.03129 -10.05439 ## [1003] -2.14036 -3.98574 1.31188 5.20165 -6.50618 6.08397 ## [1009] 3.78116 16.26293 15.53957 13.32709 8.08600 21.65631 ## [1015] 19.43212 12.09277 5.77915 16.53659 11.22533 21.66746 ## [1021] 11.42976 5.48729 -4.89110 -6.40260 -11.97573 -1.46608 ## [1027] -13.23478 -5.22255 -11.53152 -12.09905 -9.03925 -8.98008 ## [1033] -8.81654 -7.63087 0.94623 -5.72226 8.08600 5.80096 ## [1039] 4.02320 -2.13053 -2.59698 2.81535 -0.17350 -13.65513 ## [1045] 1.08165 -6.69386 -5.37754 -0.27959 -26.63622 -2.12672 ## [1051] -1.52401 -16.35382 -13.20674 -0.76411 -3.33984 -6.38283 ## [1057] 4.23107 -4.81983 -1.48599 5.16308 4.22797 -2.12672 ## [1063] 2.79049 9.61915 1.31188 1.62220 1.72631 11.27846 ## [1069] -5.66711 13.92040 1.91912 -4.36159 11.94402 12.40736 ## [1075] 16.35034 13.94798 13.03846 15.34344 14.39167 7.64232 ## [1081] 9.63009 11.57787 8.24290 8.56475 17.34525 9.76480 ## [1087] 10.80313 10.85596 0.19506 18.32632 1.72631 8.16177 ## [1093] 3.78116 -2.94581 -1.97243 0.40074 -11.97573 -1.46608 ## [1099] 0.94142 -8.31457 -9.85629 -0.76411 -11.57232 -4.82448 ## [1105] -16.91505 -6.69386 -1.48599 1.76141 -9.66092 -0.14480 ## [1111] 2.79049 7.76394 -3.71380 -0.76411 -5.23964 -7.94118 ## [1117] -8.81654 -14.18995 7.75646 2.44175 4.22797 -4.76928 ## [1123] -4.60579 -3.98574 -0.36335 4.00850 10.59206 6.60342 ## [1129] -0.26810 -1.54029 5.32423 7.20409 4.99958 -2.12672 ## [1135] 4.02320 6.52713 6.33755 -0.16753 6.15918 5.56452 ## [1141] -3.41753 5.48729 -19.97088 7.20409 1.91316 -2.12672 ## [1147] 0.94142 -0.89372 0.19506 -4.94014 6.15918 12.31736 ## [1153] 15.02910 -4.81983 -2.45888 -0.27959 12.71562 5.14032 ## [1159] -0.90765 4.67192 4.10392 9.37768 3.62611 -0.14943 ## [1165] 7.38050 2.67626 6.29712 7.20409 2.68476 -9.39375 ## [1171] -2.75672 -2.74893 -6.50584 -14.48535 -9.67252 -13.65513 ## [1177] -4.31736 -2.00880 6.29712 11.28609 0.36995 5.80096 ## [1183] 4.63956 14.56638 12.48005 12.95713 1.72631 -4.82448 ## [1189] -4.76728 1.27074 2.89201 -9.12393 -0.40166 -10.05439 ## [1195] -21.86377 -10.16978 -13.76515 -10.90589 -17.27173 -0.14943 ## [1201] -12.41588 -4.35133 -12.18777 -5.04193 -8.11771 0.51584 ## [1207] -0.90765 4.67192 -2.59698 -2.55383 -0.17350 -5.34393 ## [1213] 12.77951 12.04638 7.27001 -3.00093 -2.71647 4.47968 ## [1219] 0.32506 5.90872 4.66233 -0.76411 -5.23964 -7.42173 ## [1225] 2.88132 8.76683 -2.94532 -2.32059 15.03044 9.76480 ## [1231] 5.87227 0.34309 -3.15539 -1.95726 -3.97311 -0.14943 ## [1237] -13.76563 -7.16236 -0.02665 -7.08293 -5.80289 -0.80544 ## [1243] -7.07121 -3.98574 -10.41470 -11.50247 -0.80677 1.92837 ## [1249] -10.61621 -11.37892 1.43268 -3.68126 2.68476 8.44352 ## [1255] -5.83850 4.05351 11.36323 9.97426 18.19127 16.47296 ## [1261] 1.53157 14.85741 16.99890 18.76976 9.62921 1.83712 ## [1267] 0.94142 14.56638 14.15527 17.13316 6.15918 0.37002 ## [1273] 2.43140 6.89281 1.91912 -4.36159 -25.09301 -10.05439 ## [1279] -18.16563 -8.31457 -14.32356 -10.90589 -3.97311 -14.17458 ## [1285] -3.86744 -14.18995 2.40557 -1.64026 3.45637 9.76480 ## [1291] 5.25592 12.09277 11.92164 6.99138 30.22337 15.43406 ## [1297] 18.62844 15.32592 17.97179 9.24509 11.17242 1.83712 ## [1303] 11.41948 6.52713 10.80482 8.78110 -1.44003 -0.14943 ## [1309] -5.66711 4.55028 10.18868 11.96642 6.54279 12.40736 ## [1315] 20.04848 20.13202 20.85617 17.13316 13.75840 18.03131 ## [1321] 10.08000 12.04638 11.64801 1.08108 2.68476 -0.14480 ## [1327] -11.38571 -11.40659 -12.08993 3.41192 8.05899 -5.34393 ## [1333] -4.76728 8.29833 -8.78265 -7.76326 -9.66092 -14.67887 ## [1339] -14.46749 -11.40659 -8.73947 -12.69562 -14.73866 -12.09678 ## [1345] -3.86744 -25.43409 -5.37754 -6.40260 16.57365 4.47968 ## [1351] 4.63956 2.19830 0.19506 -5.53671 -0.17350 -2.74668 ## [1357] -4.76728 5.95580 -2.94532 -4.36159 4.22797 4.47968 ## [1363] 13.26855 -4.60415 0.75347 -2.55383 -5.23964 -8.98008 ## [1369] 1.98149 0.33373 -7.80977 -10.48460 -2.71647 -3.44800 ## [1375] -2.14036 1.57990 -5.38902 -1.36068 -3.97311 -2.74668 ## [1381] 1.53157 -3.88282 -15.10643 -4.36159 8.85760 11.74672 ## [1387] 7.72134 5.29032 2.98710 -5.53671 1.72631 2.96727 ## [1393] 8.73025 -0.13478 1.43268 -1.64026 -7.34610 -16.00015 ## [1399] -18.78199 -22.53787 -18.23242 -19.85453 -30.57036 -23.52467 ## [1405] -30.41258 -33.86719 -25.32177 -13.88627 5.77118 -2.78736 ## [1411] -8.30393 -0.89372 8.57119 3.41192 -1.44003 -1.70778 ## [1417] 0.18182 5.01879 -5.86399 1.76141 0.36995 -7.41183 ## [1423] 0.32506 2.19830 0.75347 -4.94014 2.99284 -3.26613 ## [1429] -3.86744 8.29833 2.89201 11.96642 3.45637 -8.07247 ## [1435] -2.75672 8.38234 -0.36335 -5.53671 -2.07330 0.37002 ## [1441] 6.48066 -11.37892 0.45979 1.08108 -5.80289 8.44352 ## [1447] 1.55778 -7.07776 -5.38902 5.79823 -14.10539 -15.73292 ## [1453] -7.01687 -13.25294 -15.10643 -20.00928 -20.46338 -16.66079 ## [1459] -16.93292 -13.26180 -10.97311 -5.53671 -0.17350 -6.38283 ## [1465] -2.06777 -5.75685 1.91912 -0.95992 10.40081 -2.12672 ## [1471] 10.18677 14.56638 18.06413 14.15029 19.45781 9.72012 ## [1477] 11.42976 8.29833 5.81068 17.40910 15.03044 5.80096 ## [1483] 8.33770 -5.84095 4.66233 6.99138 -2.70657 -6.90228 ## [1489] -2.06777 -17.93800 -7.80977 -0.27959 14.25883 0.51584 ## [1495] 4.63956 -4.60415 -7.62266 1.62220 -5.87291 -4.30503 ## [1501] 3.78116 -0.13478 -9.75554 -10.48460 -17.37696 -4.76928 ## [1507] -1.52401 -7.07776 -15.99878 -12.69562 -7.77271 -18.84962 ## [1513] -25.01357 -9.97340 -21.91666 -15.24694 1.14155 -10.71503 ## [1519] -8.30393 1.57990 0.19506 -16.87165 -7.77271 -7.42173 ## [1525] -8.81654 -8.56788 -7.80977 -17.28794 -7.34610 -13.35759 ## [1531] -10.76935 -13.26180 -7.62266 -13.29220 -0.80677 -0.14943 ## [1537] -2.51769 -0.13478 -9.75554 -18.64861 -11.97573 -8.73311 ## [1543] -16.93292 -16.35382 -16.55719 -13.88877 -6.50618 -13.65513 ## [1549] -4.76728 -3.41432 2.89201 -11.16493 -8.11771 -1.46608 ## [1555] 3.40685 0.96149 2.42869 -1.36068 0.45977 -6.90228 ## [1561] 6.93058 -0.13478 6.29712 7.20409 1.91316 11.08608 ## [1567] 5.87227 5.90872 5.77915 5.20165 -0.80677 6.60342 ## [1573] 2.88132 12.51488 -7.80977 -9.80427 -21.23499 -19.96399 ## [1579] -10.15300 -18.82744 -14.88197 -17.46823 -14.10539 -11.05788 ## [1585] -23.66382 -19.81202 -25.80821 -20.68961 1.14155 8.44352 ## [1591] 12.03584 10.85596 10.24641 10.57083 11.85859 11.27846 ## [1597] 13.67935 5.95580 7.75646 14.68776 1.14155 -0.14480 ## [1603] 4.63956 1.57990 -0.92176 9.37768 11.85859 11.79791 ## [1609] 4.23107 12.04638 8.24290 12.64676 2.68476 -5.42992 ## [1615] -3.98943 5.29032 8.57119 10.57083 21.35761 5.56452 ## [1621] 6.93058 12.04638 4.35135 -0.27959 -9.66092 0.51584 ## [1627] -3.98943 -7.69617 1.31188 -1.95726 8.05899 -10.01898 ## [1633] 2.43140 6.89281 -5.86399 -3.68126 -14.29054 -6.09056 ## [1639] -7.07121 -13.88021 -6.50584 -5.53671 -5.23964 -12.61623 ## [1645] -2.96761 -8.09937 -6.35043 -4.36159 -6.57450 -10.71503 ## [1651] 2.17413 -1.51213 -1.48016 -1.36068 -3.33984 8.68122 ## [1657] -12.41588 8.29833 -0.02665 -13.88627 11.94402 -13.35759 ## [1663] -7.07121 -5.84095 -3.71380 -1.95726 -7.77271 -9.49953 ## [1669] 1.53157 -29.65064 -0.99954 -10.48460 -13.51894 1.83712 ## [1675] 0.32506 -6.45936 3.54551 11.76398 -2.07330 -2.22723 ## [1681] 12.32959 -1.07179 12.13446 10.60576 -9.66092 -6.09056 ## [1687] -8.30393 -14.49861 -6.50584 -15.08192 1.09304 -4.82448 ## [1693] -7.91670 -2.94581 -1.97243 -4.36159 -4.25968 3.15840 ## [1699] -5.22214 2.19830 1.31188 6.39480 9.95879 13.87571 ## [1705] 13.22943 20.47948 9.21579 5.84342 -3.48808 3.81904 ## [1711] 1.55778 9.00075 12.48005 12.95713 15.65820 11.79791 ## [1717] 7.83042 13.92040 -0.02665 7.20409 -4.25968 -2.12672 ## [1723] 0.94142 -4.60415 2.98710 -14.48535 -7.13945 -8.98008 ## [1729] -10.61621 -15.12696 -9.75554 12.64676 10.40081 1.17648 ## [1735] 18.81576 18.27681 22.53140 17.72974 15.02493 25.82306 ## [1741] 17.72860 21.41649 18.94468 20.81077 2.68476 -0.14480 ## [1747] 9.57041 2.19830 10.24641 6.99138 12.49186 21.66746

# fast plot to check data plot(Index, type = "h", xlim = c(0, 1752), ylim = c(-35, 35))

Custom Plotting

The first full plot of our data has confirmed a viable index has been created but this graphic is not likely to be used in this form in a presentation or paper.

Some additional steps one may take to customize the graph include:

i) subsetting into a smaller time-frame i.e. the last 20 years

ii) adding colour, labels, text and lines

iii) smoothing and equation-fitting

Creating masks using logical vectors

The firts step it to subset the data and also allocate a mask to positive and negative value to be used in our colour plotting scheme.

Nb. When subsetting its always a good idea to replicate the original object/data.

# replicate data Index2 <- Index # locigal expression to show pos/neg split Index2 > 0

## [1] FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE TRUE TRUE ## [12] FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE ## [23] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [34] FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE ## [45] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE ## [56] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [67] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE ## [78] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [89] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE ## [100] TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [111] FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE TRUE ## [122] TRUE FALSE TRUE TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE ## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [144] FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE ## [155] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [166] TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [177] TRUE TRUE TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE ## [188] FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE TRUE FALSE ## [199] FALSE FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE ## [210] TRUE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE TRUE FALSE ## [221] TRUE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE ## [232] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE ## [243] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [254] TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE ## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [276] FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE ## [287] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE ## [298] TRUE TRUE FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE ## [309] FALSE TRUE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [320] TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [331] TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE ## [342] TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE TRUE FALSE TRUE ## [353] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE ## [364] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [375] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE TRUE TRUE ## [386] TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE ## [397] TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE ## [408] FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE FALSE ## [419] FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [430] FALSE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE ## [441] FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE TRUE ## [452] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [463] FALSE TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [474] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [485] FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE ## [496] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE ## [507] FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE ## [518] FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [529] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [540] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE ## [551] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [562] FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [573] FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE ## [584] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE TRUE ## [595] TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE TRUE ## [606] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [617] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [628] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [639] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [650] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE TRUE ## [661] TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE ## [672] TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE TRUE TRUE ## [683] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE ## [694] FALSE FALSE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE ## [705] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE ## [716] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [727] FALSE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE ## [738] TRUE TRUE FALSE FALSE FALSE FALSE TRUE FALSE TRUE TRUE TRUE ## [749] FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [760] TRUE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [771] TRUE FALSE TRUE FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE ## [782] FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE ## [793] TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE ## [804] TRUE FALSE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE TRUE ## [815] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE FALSE ## [826] TRUE TRUE FALSE TRUE FALSE TRUE TRUE FALSE FALSE TRUE TRUE ## [837] TRUE TRUE TRUE FALSE FALSE FALSE TRUE TRUE TRUE FALSE TRUE ## [848] FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE FALSE TRUE ## [859] FALSE TRUE TRUE FALSE FALSE TRUE TRUE TRUE FALSE TRUE TRUE ## [870] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [881] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [892] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [903] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [914] FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE ## [925] TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [936] FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE TRUE TRUE FALSE ## [947] FALSE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [958] TRUE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE ## [969] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE TRUE ## [980] TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE TRUE TRUE FALSE ## [991] TRUE FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE TRUE FALSE ## [1002] FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [1013] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [1024] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1035] TRUE FALSE TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE TRUE ## [1046] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1057] TRUE FALSE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [1068] TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE ## [1079] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1090] TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE TRUE FALSE ## [1101] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE ## [1112] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE ## [1123] FALSE FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE ## [1134] FALSE TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE ## [1145] TRUE FALSE TRUE FALSE TRUE FALSE TRUE TRUE TRUE FALSE FALSE ## [1156] FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE TRUE TRUE ## [1167] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1178] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [1189] FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1200] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE ## [1211] FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE ## [1222] FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE ## [1233] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1244] FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE TRUE TRUE ## [1255] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1266] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [1277] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [1288] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1299] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE ## [1310] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1321] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE TRUE ## [1332] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1343] FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE ## [1354] FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE TRUE TRUE FALSE ## [1365] TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE ## [1376] TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE ## [1387] TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE FALSE FALSE ## [1398] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1409] TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE FALSE ## [1420] TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE ## [1431] TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE TRUE TRUE ## [1442] FALSE TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE FALSE ## [1453] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1464] FALSE FALSE FALSE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE ## [1475] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE ## [1486] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE ## [1497] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE ## [1508] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE ## [1519] FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1530] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1541] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [1552] FALSE FALSE FALSE TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE ## [1563] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE ## [1574] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1585] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1596] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE ## [1607] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE ## [1618] TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE ## [1629] TRUE FALSE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE ## [1640] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1651] TRUE FALSE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE ## [1662] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE ## [1673] FALSE TRUE TRUE FALSE TRUE TRUE FALSE FALSE TRUE FALSE TRUE ## [1684] TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE ## [1695] FALSE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE ## [1706] TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1717] TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE ## [1728] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1739] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE ## [1750] TRUE TRUE TRUE

# assign all values in Index2 that are less than zero to the value 0 Index2[Index2 < 0] <- 0 # check assignment head(Index2)

## [1] 0.000 0.000 0.000 0.000 1.913 0.000

What have we achieved?

In technical terms we have masked all negative values in our data set by changing their original values to an arbitrary one - in this case the value 0.

We can see how this changes *the presentation of our data below. Plotting our entire index yields the following:

plot(Index2, type = "h", col = "blue")

Notice 'the gaps' in betweent the 'high-rises' - this where the negative indice values used to to be; now they have assigned to the value zero.

We will repeat this process for the positive values and then we can plot these data as the same time… assigning a different colour to each range of data.

# replicate original data Index3 <- Index # logical expression to show neg/pos split Index < 0

## [1] TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE FALSE FALSE ## [12] TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE ## [23] TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [34] TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE ## [45] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE TRUE FALSE ## [56] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE ## [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [89] FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE ## [100] FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [111] TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE TRUE FALSE ## [122] FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE TRUE TRUE TRUE ## [133] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [144] TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE FALSE ## [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [166] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [177] FALSE FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE ## [188] TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE ## [199] TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE ## [210] FALSE TRUE FALSE TRUE FALSE FALSE TRUE TRUE TRUE FALSE TRUE ## [221] FALSE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE ## [232] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE ## [243] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [254] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE ## [265] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [276] TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE ## [287] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE ## [298] FALSE FALSE TRUE FALSE TRUE TRUE FALSE FALSE TRUE TRUE TRUE ## [309] TRUE FALSE TRUE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [320] FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [331] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE ## [342] FALSE FALSE TRUE TRUE FALSE TRUE FALSE FALSE FALSE TRUE FALSE ## [353] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE ## [364] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [375] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE FALSE FALSE ## [386] FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE ## [397] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE ## [408] TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE TRUE ## [419] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [430] TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE ## [441] TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE FALSE ## [452] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [463] TRUE FALSE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [474] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [485] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE ## [496] TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE ## [507] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE ## [518] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [529] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [540] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE ## [551] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [562] TRUE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE ## [573] TRUE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE ## [584] TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE ## [595] FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE TRUE TRUE FALSE ## [606] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [617] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [628] FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [639] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [650] TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE TRUE FALSE ## [661] FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE ## [672] FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE FALSE FALSE ## [683] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE ## [694] TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE ## [705] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE ## [716] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [727] TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE ## [738] FALSE FALSE TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE ## [749] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [760] FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [771] FALSE TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE TRUE FALSE ## [782] TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE TRUE FALSE ## [793] FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE ## [804] FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE ## [815] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE TRUE ## [826] FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE ## [837] FALSE FALSE FALSE TRUE TRUE TRUE FALSE FALSE FALSE TRUE FALSE ## [848] TRUE FALSE FALSE TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE ## [859] TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE FALSE FALSE ## [870] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [881] TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [892] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [903] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [914] TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE ## [925] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [936] TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE TRUE ## [947] TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [958] FALSE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE ## [969] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE FALSE ## [980] FALSE FALSE TRUE FALSE FALSE TRUE TRUE TRUE FALSE FALSE TRUE ## [991] FALSE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE ## [1002] TRUE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [1013] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [1024] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1035] FALSE TRUE FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE FALSE ## [1046] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1057] FALSE TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE ## [1068] FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE ## [1079] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1090] FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE ## [1101] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE ## [1112] FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE ## [1123] TRUE TRUE TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE ## [1134] TRUE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE TRUE FALSE ## [1145] FALSE TRUE FALSE TRUE FALSE TRUE FALSE FALSE FALSE TRUE TRUE ## [1156] TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE ## [1167] FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1178] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [1189] TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1200] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE ## [1211] TRUE TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE ## [1222] TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE ## [1233] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1244] TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE TRUE FALSE FALSE ## [1255] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1266] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE ## [1277] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [1288] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1299] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE ## [1310] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1321] FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE ## [1332] TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1343] TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE ## [1354] TRUE TRUE TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE TRUE ## [1365] FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE ## [1376] FALSE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE ## [1387] FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE TRUE TRUE ## [1398] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1409] FALSE TRUE TRUE TRUE FALSE FALSE TRUE TRUE FALSE FALSE TRUE ## [1420] FALSE FALSE TRUE FALSE FALSE FALSE TRUE FALSE TRUE TRUE FALSE ## [1431] FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE ## [1442] TRUE FALSE FALSE TRUE FALSE FALSE TRUE TRUE FALSE TRUE TRUE ## [1453] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1464] TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE ## [1475] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE ## [1486] FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE ## [1497] TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE ## [1508] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE ## [1519] TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1530] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1541] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE ## [1552] TRUE TRUE TRUE FALSE FALSE FALSE TRUE FALSE TRUE FALSE TRUE ## [1563] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE ## [1574] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1585] TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1596] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE ## [1607] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE ## [1618] FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE TRUE ## [1629] FALSE TRUE FALSE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE ## [1640] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [1651] FALSE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE FALSE ## [1662] TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE ## [1673] TRUE FALSE FALSE TRUE FALSE FALSE TRUE TRUE FALSE TRUE FALSE ## [1684] FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE ## [1695] TRUE TRUE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE ## [1706] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1717] FALSE FALSE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE TRUE ## [1728] TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE ## [1739] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE ## [1750] FALSE FALSE FALSE

# assign all values in Index2 that are greater than zero to the value 0 Index3[Index3 > 0] <- 0 # check assignment head(Index3)

## [1] -5.6671 -0.6033 -3.9182 -4.3616 0.0000 -4.7693

plot(Index3, type = "h", col = "red")

Choosing a time-frame

Now that we have two sets of masked values for our data. We can suitably subset them for a shorter time period.

Let's look at data for the last 12 years (2000-2011).

# data of pos values Short_Pos <- Index2[1609:1752] Short_Pos

## [1] 4.2311 12.0464 8.2429 12.6468 2.6848 0.0000 0.0000 5.2903 ## [9] 8.5712 10.5708 21.3576 5.5645 6.9306 12.0464 4.3513 0.0000 ## [17] 0.0000 0.5158 0.0000 0.0000 1.3119 0.0000 8.0590 0.0000 ## [25] 2.4314 6.8928 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 ## [33] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 ## [41] 0.0000 0.0000 2.1741 0.0000 0.0000 0.0000 0.0000 8.6812 ## [49] 0.0000 8.2983 0.0000 0.0000 11.9440 0.0000 0.0000 0.0000 ## [57] 0.0000 0.0000 0.0000 0.0000 1.5316 0.0000 0.0000 0.0000 ## [65] 0.0000 1.8371 0.3251 0.0000 3.5455 11.7640 0.0000 0.0000 ## [73] 12.3296 0.0000 12.1345 10.6058 0.0000 0.0000 0.0000 0.0000 ## [81] 0.0000 0.0000 1.0930 0.0000 0.0000 0.0000 0.0000 0.0000 ## [89] 0.0000 3.1584 0.0000 2.1983 1.3119 6.3948 9.9588 13.8757 ## [97] 13.2294 20.4795 9.2158 5.8434 0.0000 3.8190 1.5578 9.0007 ## [105] 12.4800 12.9571 15.6582 11.7979 7.8304 13.9204 0.0000 7.2041 ## [113] 0.0000 0.0000 0.9414 0.0000 2.9871 0.0000 0.0000 0.0000 ## [121] 0.0000 0.0000 0.0000 12.6468 10.4008 1.1765 18.8158 18.2768 ## [129] 22.5314 17.7297 15.0249 25.8231 17.7286 21.4165 18.9447 20.8108 ## [137] 2.6848 0.0000 9.5704 2.1983 10.2464 6.9914 12.4919 21.6675

# data of neg values Short_Neg <- Index3[1609:1752] Short_Neg

## [1] 0.00000 0.00000 0.00000 0.00000 0.00000 -5.42992 -3.98943 ## [8] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 ## [15] 0.00000 -0.27959 -9.66092 0.00000 -3.98943 -7.69617 0.00000 ## [22] -1.95726 0.00000 -10.01898 0.00000 0.00000 -5.86399 -3.68126 ## [29] -14.29054 -6.09056 -7.07121 -13.88021 -6.50584 -5.53671 -5.23964 ## [36] -12.61623 -2.96761 -8.09937 -6.35043 -4.36159 -6.57450 -10.71503 ## [43] 0.00000 -1.51213 -1.48016 -1.36068 -3.33984 0.00000 -12.41588 ## [50] 0.00000 -0.02665 -13.88627 0.00000 -13.35759 -7.07121 -5.84095 ## [57] -3.71380 -1.95726 -7.77271 -9.49953 0.00000 -29.65064 -0.99954 ## [64] -10.48460 -13.51894 0.00000 0.00000 -6.45936 0.00000 0.00000 ## [71] -2.07330 -2.22723 0.00000 -1.07179 0.00000 0.00000 -9.66092 ## [78] -6.09056 -8.30393 -14.49861 -6.50584 -15.08192 0.00000 -4.82448 ## [85] -7.91670 -2.94581 -1.97243 -4.36159 -4.25968 0.00000 -5.22214 ## [92] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 ## [99] 0.00000 0.00000 -3.48808 0.00000 0.00000 0.00000 0.00000 ## [106] 0.00000 0.00000 0.00000 0.00000 0.00000 -0.02665 0.00000 ## [113] -4.25968 -2.12672 0.00000 -4.60415 0.00000 -14.48535 -7.13945 ## [120] -8.98008 -10.61621 -15.12696 -9.75554 0.00000 0.00000 0.00000 ## [127] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 ## [134] 0.00000 0.00000 0.00000 0.00000 -0.14480 0.00000 0.00000 ## [141] 0.00000 0.00000 0.00000 0.00000

Overplotting two sets of data using par

Now we can plot our two sets of data using an overplotting technique:

# set axis limits and type of plot plot(Short_Pos, ylim = c(-30, 30), type = "h", lwd = 2.9, col = "blue") # set overplotting to 'on' par(new = T) plot(Short_Neg, ylim = c(-30, 30), type = "h", lwd = 2.9, col = "red")

Cleaning up the plot

Without any alteration our overplotting produces two sets of labels on the y-axis.

Lets clean up the plot, by removing default annotations and adding our own labels.

# set axis limits and type of plot turn off axes and labels for this plot plot(Short_Pos, ylim = c(-32, 32), type = "h", lwd = 2.9, col = "blue", ylab = "", axes = FALSE, ann = F) # set overplotting to 'on' par(new = T) # turn off defualt axes set new ylab set new xlab plot(Short_Neg, ylim = c(-32, 32), type = "h", lwd = 2.9, col = "red", font.lab = 2, axes = FALSE, ylab = "Southern Oscillation Index", xlab = "Years") # custom x axis at beginning of each year axis(1, at = c(0, 13, 25, 37, 49, 61, 73, 85, 97, 109, 121, 133, 145), label = c("2000", "2001", "2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011", "2012"), font = 2) # default y axis with bold horizontal labels axis(2, font = 2, las = 1) # keep overploting 'on' par(new = T) # add zero line across xaxis abline(h = 0, col = "black") # add horizontal dotted lines to plot abline(h = 15, col = "16", lty = 3) abline(h = -15, col = "16", lty = 3) # add some text to plot text(37, -32.5, "Source: http://www.cgd.ucar.edu/cas/catalog/climind/soiAnnual.html", cex = 0.75) # add box around plot box()

# turn overplotting 'off' par(new = F)

[r-module-geo-blog]

Chapter Seven: Sea Ice Decline

In this chapter we will learn how to access data-sets directly from the Internet. Then we will use several looping functions to make some calculations and produce our final output. Lets get started! :D

#url<-“http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv”

# read in data associating missing values with numeric -9999 all_data <- read.csv("plot.csv", na.strings = c(-9999)) # inspect data structure head(all_data)

## Month Day X1980.s.Average X1990.s.Average X2000.s.Average X2002 X2003 ## 1 1 1 13177778 13025500 12486531 NA 12787813 ## 2 1 2 13208924 13044813 12536141 NA 12856875 ## 3 1 3 13279219 13080563 12594219 NA 12936250 ## 4 1 4 13357366 13129188 12646688 NA 13063281 ## 5 1 5 13404196 13160125 12689953 NA 13182813 ## 6 1 6 13445714 13212563 12730125 NA 13249063 ## X2004 X2005 X2006 X2007 X2008 X2009 X2010 X2011 ## 1 12708906 12287969 12325625 12344375 12322031 12286250 12255313 12008906 ## 2 12743750 12335625 12381250 12395156 12383281 12373125 12279063 12082813 ## 3 12760781 12395938 12458438 12440469 12448281 12465938 12327656 12121563 ## 4 12847969 12539063 12550938 12520469 12520000 12548750 12407031 12243594 ## 5 12918125 12596406 12592344 12586563 12561250 12589531 12477344 12369531 ## 6 12902188 12595313 12560781 12578750 12611250 12637031 12537031 12410938 ## X2012 ## 1 12524531 ## 2 12562031 ## 3 12553281 ## 4 12536094 ## 5 12557813 ## 6 12632969

tail(all_data)

## Month Day X1980.s.Average X1990.s.Average X2000.s.Average X2002 ## 361 12 26 12899102 12729931 12230453 12436875 ## 362 12 27 12923164 12750375 12264766 12477813 ## 363 12 28 12955898 12797000 12298328 12516719 ## 364 12 29 12993090 12833938 12335969 12550938 ## 365 12 30 13033958 12871938 12383063 12635313 ## 366 12 31 13099479 12913563 12415532 12720625 ## X2003 X2004 X2005 X2006 X2007 X2008 X2009 ## 361 12495938 12150313 12091094 12050313 12055625 12141719 12087344 ## 362 12504219 12219688 12149375 12075625 12101406 12116719 12153125 ## 363 12554219 12324531 12183594 12102344 12170781 12133281 12204688 ## 364 12586406 12320938 12194688 12168125 12236563 12181250 12232969 ## 365 12618750 12285156 12236094 12245469 12304531 12207656 12252656 ## 366 12655938 12265000 12282188 12290000 12369688 12230156 12255469 ## X2010 X2011 X2012 ## 361 11631406 12036875 NA ## 362 11653125 12076719 NA ## 363 11717656 12108438 NA ## 364 11799375 12295156 NA ## 365 11865781 12387656 NA ## 366 11914688 12455156 NA

The data is arranged by each successive date in each month. These data include averages from the 1980's/1990's/2000's and yearly data from 2002-2012. Lets firstly separate these data into two distinct files for clarity.

# sub-select all rows, columns 1 to 5 avg_data <- all_data[, c(1:5)] # sub-select all rows, columns 1,2,6 to 16 yearly_data <- all_data[, c(1, 2, 6:16)]

Lets add a count of the days in the year for the yearly data, and additionally convert the sea-ice extent into the units of millions of km2.

# count of days added as additional column in yearly_data for (i in 1:366) yearly_data$days[i] <- i # convert all data to mil km2 avg_data[, c(3:5)] <- avg_data[, c(3:5)]/1e+06 yearly_data[, c(3:13)] <- yearly_data[, c(3:13)]/1e+06

We will now calculate the minimum sea-ice extent for each individual year (2002-2012), and the specific time in the year at which the minimum extent occurred.

We use a loop structure to execute our commands for each year. The output will be stored in individual vectors.

# create vectors for output of type numeric and character min_asie <- numeric(11) min_lab <- numeric(11) min_r <- numeric(11) min_d <- numeric(11) min_m <- numeric(11) min_date <- character(11)

# loop will be executed 11 times for (j in 2002:2012) { # assign counters for column selection and vector assignment which_yr <- j col_sel <- which_yr - 1999 cnt <- j - 2001 # col_sel selects column to compute minimum on na.rm = TRUE removes na's # from data min_asie[cnt] <- min(yearly_data[, col_sel], na.rm = TRUE) # round asie result to 2 digits be used later as a label in plot min_lab[cnt] <- round(min_asie[cnt], 2) # find the row number of min extent min_r[cnt] <- which(yearly_data[, col_sel] == min_asie[cnt]) # find the day of minimum using created min_r index min_d[cnt] <- (yearly_data[min_r[cnt], 2]) # find month of minimum using min_r index min_m[cnt] <- (yearly_data[min_r[cnt], 1]) # paste date results in a character string min_date[cnt] <- paste(min_m[cnt], "/", min_d[cnt], "/", j, sep = "") } # reset counters for later use in plots! j <- 2002 which_yr <- 2002 col_sel <- which_yr - 1999 cnt <- j - 2001

We have achieved alot in the loop :)

Lets review the implementation:

The counters enabled us to

i) select the appropriate column from our yearly data to make each minimum calculation (using col_sel)

ii) assign our results to a specific position in our vectors (using cnt)

Indexing i.e. [ ] and select functions which and == enabled us to

iii) 'find' and select certain portions of our data… i.e. the row of the minimum ice extent

iv) and were then used to extract the associated dates when this occurred

head(yearly_data)

## Month Day X2002 X2003 X2004 X2005 X2006 X2007 X2008 X2009 X2010 X2011 ## 1 1 1 NA 12.79 12.71 12.29 12.33 12.34 12.32 12.29 12.26 12.01 ## 2 1 2 NA 12.86 12.74 12.34 12.38 12.40 12.38 12.37 12.28 12.08 ## 3 1 3 NA 12.94 12.76 12.40 12.46 12.44 12.45 12.47 12.33 12.12 ## 4 1 4 NA 13.06 12.85 12.54 12.55 12.52 12.52 12.55 12.41 12.24 ## 5 1 5 NA 13.18 12.92 12.60 12.59 12.59 12.56 12.59 12.48 12.37 ## 6 1 6 NA 13.25 12.90 12.60 12.56 12.58 12.61 12.64 12.54 12.41 ## X2012 days ## 1 12.52 1 ## 2 12.56 2 ## 3 12.55 3 ## 4 12.54 4 ## 5 12.56 5 ## 6 12.63 6

Lets construct some preliminary axes for our plots.

# names of the months month_names <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Se", "Oct", "Nov", "Dec") # tick-marks for beginning of month month_beg <- c(1, 32, 60, 91, 121, 151, 182, 213, 244, 274, 305, 335, 366) # label-marks for middle of month month_pos <- c(16, 46, 75, 106, 135, 165, 200, 228, 255, 289, 320, 355)

Lets lay out the base plot:

# no plot of data, set up yaxis + labels only yy <- "Historical Se Ice Extent - 10e6sq kilometrere" plot(yearly_data[, 14], type = "n", col = "blue", axes = F, xlab = "", ylab = yy, ylim = c(0, 20), xaxs = "i", yaxs = "i") # designate labels at intervals along x axis axis(side = 1, at = month_beg, labels = F, xaxs = "i") axis(side = 1, at = month_pos, labels = month_names, tick = F, line = F, xaxs = "i") # y-axis labels at default positions axis(side = 2, at = NULL, labels = T, yaxs = "i", las = 1) # add minimum markers for each year points(yearly_data[min_r, 14], min_asie, col = "orange", pch = 3, cex = 1)

The last requirement is to add the data to our plot. Lets add one year to illustrate:

# no plot of data, set up y-axis + labels plot(yearly_data[, 14], type = "n", col = "blue", axes = F, xlab = "", ylab = "Historical Arctic Sea ice Extent - 10e6 Sq Km", ylim = c(0, 20), xaxs = "i", yaxs = "i") # designate labels at intervals along x-axis axis(side = 1, at = month_beg, labels = F, xaxs = "i") axis(side = 1, at = month_pos, labels = month_names, tick = F, line = F, xaxs = "i") # y-axis labels at defualt positions axis(side = 2, at = NULL, labels = T, yaxs = "i", las = 1) # add minimum markers for each year points(yearly_data[min_r, 14], min_asie, col = "orange", pch = 3, cex = 1) legend("bottomleft", c("minima", "2002"), cex = 0.8, col = c("orange", "blue"), pch = 3, lty = 2) ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### Re-estlablish the 'base plot' as above ### ### ### ### ### # add data using points plotting function points(yearly_data[, 14], yearly_data[, 3], col = "blue", type = "l", lwd = 1.5)

Above we re-plotted the entire 'base plot' and overlayed one year of data, 2002. We can see there are gaps in this data (no coverage), and so a straightforward continuous line plot would be a good choice, as its easy to spot these data gaps.

We can repeat this process for all years, by using a looping sequence similar to that which calculated our minimum sea-ice extent and minimum dates.

# no plot of data, set up yaxis + labels only plot(yearly_data[, 14], type = "n", col = "blue", axes = F, xlab = "", ylab = "Sea Ice Extent - 10e6 Sq Km", main = "Recent Historical Arctic Sea Ice Extent, 2002-2012", cex.main = 1.2, font.main = 3, ylim = c(0, 20), xaxs = "i", yaxs = "i", font.lab = 2) # designate labels at intervals along x axis axis(side = 1, at = month_beg, labels = F, xaxs = "i") axis(side = 1, at = month_pos, labels = month_names, tick = F, line = F, xaxs = "i") # y-axis labels at default positions axis(side = 2, at = NULL, labels = T, yaxs = "i", las = 1) # add minimum markers for each year # points(yearly_data[min_r,14],min_asie,col='orange',pch=15,cex=2) my.col <- rainbow(11) new_sm <- round(min_asie, digits = 2) new_sm <- as.character(new_sm) legend("bottomleft", title = "Minimum Extent", c(paste("2002", new_sm[1]), paste("2003", new_sm[2]), paste("2004", new_sm[3]), paste("2005", new_sm[4]), paste("2006", new_sm[5]), paste("2007", new_sm[6]), paste("2008", new_sm[7]), paste("2009", new_sm[8]), paste("2010", new_sm[9]), paste("2011", new_sm[10]), paste("2012", new_sm[11])), cex = 0.8, col = c(my.col), lty = 1) ############################## my.col

## [1] "#FF0000FF" "#FF8B00FF" "#E8FF00FF" "#5DFF00FF" "#00FF2EFF" ## [6] "#00FFB9FF" "#00B9FFFF" "#002EFFFF" "#5D00FFFF" "#E800FFFF" ## [11] "#FF008BFF"

# loop will be executed 11 times for (j in 2002:2012) { # assign counters for column selection and vector assignment which_yr <- j col_sel <- which_yr - 1999 cnt <- j - 2001 # plot data as lines with individual colours points(yearly_data[, col_sel], col = my.col[cnt], type = "o", pch = "4", cex = 0.5) }

Thats it!

[r-module-geo-blog]

Statistics and Spatial Data

R' real strength, arguably is its immense focus on statistical functionality. Statistics and applied statistics remain the most popular and documented aspect of the R environment. Several comprehensive add-on packages provide additional functionality to base stats, and many more exist for specific statistics requirements.

Here, is just a short introduction to some 'elementary statistics'.

Correlation and Simple Linear Fitting

Lets load in the required MASS (link) package, which contains sample data.

# load MASS package library(MASS) # load sample data data(cats) # check structure of data-set str(cats)

## 'data.frame': 144 obs. of 3 variables: ## $ Sex: Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ... ## $ Bwt: num 2 2 2 2.1 2.1 2.1 2.1 2.1 2.1 2.1 ... ## $ Hwt: num 7 7.4 9.5 7.2 7.3 7.6 8.1 8.2 8.3 8.5 ...

attach(cats)

## The following object(s) are masked from 'cats (position 4)': ## ## Bwt, Hwt, Sex

# plot data bodyweight vs height plot(Bwt, Hwt)

Correlations

We can compute a correlation of the data simply:

# compute correlation cor(Bwt, Hwt)

## [1] 0.8041

To run a complete correlation computation, i.e. including a t-test, and confidence intervals we can:

# complete correlation test my_cor <- cor.test(Bwt, Hwt)

The confidence level is set by default to 95%, we can change this using:

# show structure of object str(my_cor)

## List of 9 ## $ statistic : Named num 16.1 ## ..- attr(*, "names")= chr "t" ## $ parameter : Named int 142 ## ..- attr(*, "names")= chr "df" ## $ p.value : num 0 ## $ estimate : Named num 0.804 ## ..- attr(*, "names")= chr "cor" ## $ null.value : Named num 0 ## ..- attr(*, "names")= chr "correlation" ## $ alternative: chr "two.sided" ## $ method : chr "Pearson's product-moment correlation" ## $ data.name : chr "Bwt and Hwt" ## $ conf.int : atomic [1:2] 0.738 0.855 ## ..- attr(*, "conf.level")= num 0.95 ## - attr(*, "class")= chr "htest"

# set new conf.level cor.test(Bwt, Hwt, conf.level = 0.99)

## ## Pearson's product-moment correlation ## ## data: Bwt and Hwt ## t = 16.12, df = 142, p-value < 2.2e-16 ## alternative hypothesis: true correlation is not equal to 0 ## 99 percent confidence interval: ## 0.7130 0.8685 ## sample estimates: ## cor ## 0.8041

NB. To investigate subsets within our data, we need to arrange our arguments (Bwt and Hwt) slightly differently. We use the formula that is applied when calculating Linear Models (see below), BUT with both 'terms' on the RHS of the tidle (~).

# compute correlation test of subset of data NOTE: both terms on the RHS # of the ~ cor.test(~Bwt + Hwt, subset = (Sex == "F"))

## ## Pearson's product-moment correlation ## ## data: Bwt and Hwt ## t = 4.215, df = 45, p-value = 0.0001186 ## alternative hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## 0.2890 0.7106 ## sample estimates: ## cor ## 0.532

Introduction to Linear Models

Creating linear models with R is simple.

Specifically a linear model will be formulated in the form below:

response ~ term(s)

where the response is modelled by the linear predictor term(s).

The response is on the LHS of the ~. The term(s) are on the RHS of the ~.

I.e. 'y' modelled by 'x'.

# compute linear model my_lm <- lm(Hwt ~ Bwt) # print results my_lm

## ## Call: ## lm(formula = Hwt ~ Bwt) ## ## Coefficients: ## (Intercept) Bwt ## -0.357 4.034

Note: The LHS is the response (Height!) and the RHS, the predictor term (Body Weight!).

we stored (assigned) our linear model and correlation test results to the data objects, my_lm and my_cor respectively. We can exploit this data object in further plots/analysis.

# plot data again plot(Bwt, Hwt) # add linear model result abline(my_lm, col = "red", lty = 2) # add text of correlation statistic using data object my_cor text(2.5, 18.5, "Correlation =", font = 2) text(2.97, 18.5, my_cor$estimate, font = 2)

Multiple Regressions

Computing multiple regressions in R is straightforward.

They share the same general formula,i.e the arguments are arranged in the same way as a linear model.

Specifically:

y ~ x => y is modelled by x (linear)

y ~ x + x2 => is modelled by x and x2 (multiple)

y ~ x + x2 + x3 => y is modelled by x, x2, and x3 (multiple)

y ~ x * x2 => y is modelled by x and x2 AND the interaction between them (multiple)

# load in sample data new <- read.table("mlr04.csv", sep = ",", header = T) head(new)

## X1 X2 X3 X4 ## 1 85.1 8.5 5.1 4.7 ## 2 106.3 12.9 5.8 8.8 ## 3 50.2 5.2 2.1 15.1 ## 4 130.6 10.7 8.4 12.2 ## 5 54.8 3.1 2.9 10.6 ## 6 30.3 3.5 1.2 3.5

attach(new)

## The following object(s) are masked from 'new (position 4)': ## ## X1, X2, X3, X4

linear <- lm(X1 ~ X2) plot(X2, X1)

summary(linear)

## ## Call: ## lm(formula = X1 ~ X2) ## ## Residuals: ## Min 1Q Median 3Q Max ## -13.14 -9.03 -3.69 3.21 29.72 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 15.51 11.60 1.34 0.21799 ## X2 7.98 1.22 6.52 0.00018 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 14.3 on 8 degrees of freedom ## Multiple R-squared: 0.842, Adjusted R-squared: 0.822 ## F-statistic: 42.5 on 1 and 8 DF, p-value: 0.000184

# create multiple regression multiple <- lm(X1 ~ X2 + X3) summary(multiple)

## ## Call: ## lm(formula = X1 ~ X2 + X3) ## ## Residuals: ## Min 1Q Median 3Q Max ## -15.417 -2.570 0.805 2.120 11.046 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 11.85 6.76 1.75 0.1233 ## X2 4.23 1.15 3.67 0.0080 ** ## X3 7.44 1.81 4.12 0.0045 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 8.24 on 7 degrees of freedom ## Multiple R-squared: 0.954, Adjusted R-squared: 0.941 ## F-statistic: 72.1 on 2 and 7 DF, p-value: 2.13e-05

plot(multiple)

multiple$coeff["X2"]

## X2 ## 4.228

[r-module-geo-blog]

Sources, Credit and Appendix of Commands

Many Thanks…

Without the reference and insight of several key contributors/sources this short guide would not be possible.

My sincere thanks go out to all those who have contributed by sharing their experiences, advice and insight into this great tool, R!!

Important sources include:

Quick R

Revolutions

Using R

Rtips Revival 2012!

R Tutorials: Carolina Coastal University

A special thanks to Kelly O'Day at Climate Charts and Graphs. From which specific inspiration and large parts of the implementation of Chapter 7 were originally derived.

Appendix of Commands

Some commands used within this module are include here, in alphabetical order.

Each with a brief description and the arguments required.

abline(h=15,col=“red”,lty=“3”) - add a horizontal line

apropos(“plot”) - list all functions with 'plot' in title

c(i,ii,iii) - concatenate a vector - i can be numeric, character or LOGICAL

cwd() - show the current working directory

dim(my.array) <-c(x,y) - assign my.array dimensions of x and y

getwd() - get current working directory

head(my.data) - show first rows of object my.data

hist(new.data$temp) - histogram of new.data, column 'temp'

history() - display last 25 commands

length(my.vector) - show length of my.vector

load(file=“some.data.Rdata”) - load workspace objects

load.history(file=“myfile”) - load command history

ls() - list all objects in the workspace

mean(my.vector) - calculate sum of my.vector

nrow(x) - show no of rows in array x

ncol(X) - show no of cols in array x

numeric(x) - create numeric vector of length x

par(new=T) - set overplotting to on, use after plot() command

par(new=F) - set overplotting to off

plot(one,two) - simple x.y plot of object one and two

points(mydata_x,mydata_y,col=“red”) - add points to a xy plot

read.table(“mydata.csv”,sep=“,”,header=T) - load csv file, with header and sep “,”

rep(i:ii,iii) - replicate vector seq i:ii, iii times

save.history(file=“myfile”) - save command history

save.image() - save all workspace objects to file .Rdata

scan(mytextfile.ascii) - scan in .ascii file one item at a time

scan(mytextfile.txt) - scan in .txt file one item at a time

sd(my.vector) - calculate std dev of my.vector

seq(i,ii,by=iii) - create a vector from i to ii, by iii

str(my.object) - list attributes of object my.object

sum(my.vector) - calculate sum of my.vector

tail(my.data) - show last rows of object my.data

text(x,y,“my.message”“) - add text plot at pos x and y, text in ”“