第九届生物多样性会议(2010-11-4, 福建, 厦门)
An brief introduction to R , Package and
Graphics used in biodiversity research
Kechang NIU (牛克昌)
Department of Ecology , Peking University
Outline
What this is
Section 1 >-A short, highly incomplete tour around some
of the basic concepts of R as a programming language
Section 2 >-Some hints on how to obtain documentation
on the many library functions (packages)
Section 3>-A brief introduction of R graphics
Section 1
R, S and S-plus
S: an interactive environment for data analysis
developed at Bell Laboratories since 1976
1988 - S2: RA Becker, JM Chambers, A Wilks
1992 - S3: JM Chambers, TJ Hastie
1998 - S4: JM Chambers
Exclusively licensed by AT&T/Lucent to Insightful
Corporation, Seattle WA. Product name: “S-plus”.
Implementation languages C, Fortran.
See:
http://cm.bell-labs.com/cm/ms/departments/sia/S/history.html
Section 1
R, S and S-plus
R: initially written by Ross Ihaka and Robert Gentleman
at Dep. of Statistics of U of Auckland, New Zealand
during 1990s.
Since 1997: international “R-core” team of ca. 15 people
with access to common CVS archive.
GNU General Public License (GPL)
- can be used by anyone for any purpose
- contagious
Open Source
-quality control!
-efficient bug tracking and fixing system supported by
the user community
Section 1
What R is
What R is not
o data handling and
storage: numeric, textual
o is not a database, but
connects to DBMSs
o matrix algebra
o has no graphical user
interfaces, but connects
to Java, TclTk
o hash tables and regular
expressions
o high-level data analytic
and statistical functions
o classes (“OO”)
o graphics
o programming language:
loops, branching,
subroutines
o language interpreter
can be very slow, but
allows to call own C/C++
code
o no spreadsheet view of
data, but connects to
Excel/MsOffice
o no professional /
commercial support
Section 1
R Packaging and statistics
o Packaging: a crucial infrastructure to efficiently
produce, load and keep consistent software libraries
from (many) different sources / authors
o Statistics: most packages deal with statistics and data
analysis
o State of the art: many statistical researchers provide
their methods as R packages
Section 1
R as a calculator
> seq(0, 5, length=6)
0.0
-0.5
[1] 1.414214
-1.0
> sqrt(2)
sin(seq(0, 2 * pi, length = 100))
[1] 5
0.5
1.0
> log2(32)
0
20
40
[1] 0 1 2 3 4 5
> plot(sin(seq(0, 2*pi, length=100)))
60
Index
80
100
Section 1
Variables
> a = 49
> sqrt(a)
[1] 7
> a = "The dog ate my homework"
> sub("dog","cat",a)
[1] "The cat ate my homework“
> a = (1+1==3)
>a
[1] FALSE
Numeric
Character
string
Logical
Section 1
Missing values
Variables of each data type (numeric, character, logical)
can also take the value NA: not available.
>- NA is not the same as 0
>- NA is not the same as “”
>- NA is not the same as FALSE
Any operations (calculations, comparisons) that involve NA
> NA | TRUE
may or may not produce NA:
[1] TRUE
> NA==1
> NA & TRUE
[1] NA
[1] NA
> 1+NA
[1] NA
> max(c(NA, 4, 7))
[1] NA
Section 1
Functions and Operators
Functions do things with data
“Input”: function arguments (0,1,2,…)
“Output”: function result (exactly one)
Example:
add = function(a,b)
{ result = a+b
return(result) }
Operators:
Short-cut writing for frequently used functions
of one or two arguments.
Examples: + - * / ! & | %%
Section 1
Vectors and Matrices
Vector: an ordered collection of data of the
same type
> a = c(1,2,3)
> a*2
[1] 2 4 6
Matrix: a rectangular table of data of the same
type
Example: the expression values for 10000 genes
for 30 tissue biopsies: a matrix with 10000
rows and 30 columns.
Section 1
Data frames
Data frame: is supposed to represent the typical data
table that researchers come up with – like a
spreadsheet.
It is a rectangular table with rows and columns; data
within each column has the same type (e.g. number, text,
logical), but different columns may have different
types. And each row on one case
Example:
>a
localisation tumorsize progress
XX348 proximal
6.3 FALSE
XX234
distal
8.0 TRUE
XX987 proximal
10.0 FALSE
Section 1
Loops
When the same or similar tasks need to be
performed multiple times; for all elements of a
list; for all columns of an array; etc.
for(i in 1:10) {
print(i*i)
}
i=1
while(i<=10) {
print(i*i)
i=i+sqrt(i)
}
Section 1
lapply, sapply, apply
When the same or similar tasks need to be performed multiple
times for all elements of a list or for all columns of an array. May be
easier and faster than “for” loops
lapply( li, fct )
To each element of the list li, the function fct is applied. The
result is a list whose elements are the individual fct results.
>
>
>
>
>
>
>
>
li = list("klaus","martin","georg")
lapply(li, toupper)
[[1]]
[1] "KLAUS"
[[2]]
[1] "MARTIN"
[[3]]
[1] "GEORG"
Section 1
Regular expressions
A tool for text matching and replacement which is available in similar
forms in many programming languages (Perl, Unix shells, Java)
> a = c("CENP-F","Ly-9", "MLN50", "ZNF191", "CLH-17")
> grep("L", a)
[1] 2 3 5
> grep("L", a, value=T)
[1] "Ly-9"
"MLN50" "CLH-17"
> grep("^L", a, value=T)
[1] "Ly-9"
> grep("[0-9]", a, value=T)
[1] "Ly-9"
"MLN50" "ZNF191" "CLH-17"
> gsub("[0-9]", "X", a)
[1] "CENP-F" "Ly-X"
"MLNXX"
"ZNFXXX" "CLH-XX"
Section 1
Storing data
Every R object can be stored into and restored
from a file with the commands
“save” and “load”.
This uses the XDR (external data representation)
standard of Sun Microsystems and others, and is
portable between MS-Windows, Unix, Mac.
> save(x, file=“x.Rdata”)
> load(“x.Rdata”)
Section 1
Importing and exporting data
There are many ways to get data into R and out of R.
Most programs (e.g. Excel), as well as humans, know how
to deal with rectangular tables in the form of tabdelimited text files.
> x = read.delim(“filename.txt”)
also: read.table, read.csv
> write.table(x, file=“x.txt”, sep=“\t”)
Section 1
Statistical functions
t.test()
Student (test t), determines if the averages of two
populations are statistically different.
prop.test()
hypothesis tests with proportions
Non-parametrical tests
kruskal.test()
Kruskal-Wallis test (variance analysis)
chisq.test() 2 test for convergence
ks.test()
Kolmogorov-Smirnov test
Section 1
lm()
> pisa.lm <- lm(inclin ~ year)
> pisa.lm
> summary(pisa.lm)
Call:
lm(formula = inclin ~ year)
Residuals:
Min
1Q Median
3Q Max
-5.9670 -3.0989 0.6703 2.3077 7.3956
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -61.1209 25.1298 -2.432 0.0333 *
year
9.3187 0.3099 30.069 6.5e-12 ***
--Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 4.181 on 11 degrees of freedom
Multiple R-Squared: 0.988,
Adjusted R-squared: 0.9869
F-statistic: 904.1 on 1 and 11 DF, p-value: 6.503e-012
Section 1
Getting help
Details about a specific command whose name you
know (input arguments, options, algorithm, results):
>? t.test
or
>help(t.test)
Section 1
Getting help
o
HTML search
engine
o
search for topics
with regular
expressions:
“help.search”
Web sites
Section 1
www.r-project.org
cran.r-project.org
www.bioconductor.org
Full text search:
www.r-project.org
or
www.google.com
with ‘… site:.r-project.org’ or other R-specific
keywords
Section 2
R Packages
• There are many contributed
packages that can be used to
extend R.
• These libraries are created and
maintained by the authors.
Section 2
Installing R packages from CRAN
• install.packages(‘packageName’)
• OR in Rgui:
– select a local repository (if needed)
– select package(s) from list
Section 2
Ade4>-Analysis of Ecological Data : Exploratory
and Euclidean methods in Environmental sciences
http://pbil.univlyon1.fr/ADE4/home.php?lang=eng
ade4 is characterized by :
-the implementation of graphical and statistical functions
- the availability of numerical data
- the redaction of technical and thematic documentation
- the inclusion of bibliographic references
Section 2
Vegan: R functions for vegetation ecologists
Vegan: R Labs for Vegetation Ecologists
http://ecology.msu.montana.edu/labdsv/R/labs/
•Familiarization with Data
Lab 1 Loading Vegetation Data and Simple Graphical Data Summaries
Lab 2 Loading Site/Environment Data and Simple Graphical Summaries
Lab 3 Vegetation Tables and Summaries
•Modeling Species Distributions
Lab 4 Modeling Species Distributions with Generalized Linear Models
Lab 5 Modeling Species Distributions with Generalized Additive Models
Lab 6 Modeling Species Distributions with Classification Trees
•Ordination
Lab 7 Principal Components Analysis
Lab 8 Principal Coordinates Analysis
Lab 9 Nonmetric Multi-Dimensional Scaling
Lab 10 Correspondence Analysis and Detrended Corresponence Analysis
Lab 11 Fuzzy Set Ordination
Lab 12 Canonical Correspondence Analysis
•Cluster Analysis
Lab 13 Cluster Analysis
Lab 14 Discriminant Analysis with Tree Classifiers
•Introduction
R for Ecologists, a primer on the S language and available software
Analysis of Phylogenetics
and Evolution
http://ape.mpl.ird.fr/
•ape is not a classical software (i.e., a "black box" though it can be used as
one) but an environment for data analyses and development of new methods
where interactivity and user's decisions are left an important place.
•ape is written in R, but some functions for computer-intensive tasks are
written in C.
•ape is available for all main computer operating systems (Linux, Unix,
Windows, MacOS X 10.4 and later).
•ape is distributed under the terms of the GNU General Public Licence,
meaning that it is free and can be freely modified and redistribued (under
some conditions).
smatr:(Standardised) Major Axis Estimation and Testing Routines
This package provides methods of fitting bivariate lines in allometry using
the major axis (MA) or standardised major axis (SMA), and for making
inferences about such lines.
BiodiversityR: GUI for biodiversity and community
ecology analysis
This package provides a GUI and some utility functions
(often based on the vegan package) for statistical analysis
of biodiversity and ecological communities,
• species accumulation curves,
• diversity indices, Renyi profiles,
• GLMs for analysis of species abundance and presenceabsence
• distance matrices,
• Mantel tests, and cluster, constrained and unconstrained
ordination analysis. A book on biodiversity and community
ecology analysis is available for free download from the
website.
untb: ecological drift under the UNTB
A collection of utilities for biodiversity data. Includes the simulation of
ecological drift under Hubbell's Unified Neutral Theory of Biodiversity, and the
calculation of various diagnostics such as Preston curves.
Functions (34)
alonso Various functions from Alonso and McKane
bci Barro Colorado Island (BCI) dataset
butterflies
abundance data for butterflies
etienne Etienne's sampling formula
expected.abundance
Expected abundances under the neutral model
plot.count Abundance curves
rand.neutral
Random neutral ecosystem
volkov Expected frequency of species
zsm Zero sum multinomial distribution as derived by McKane 2004
sem:
Structural Equation Models
sem
is an R package for fitting structural-equation models. The package supports general
structural equation models with latent varibles, fit by maximum likelihood assuming
multinormality, and single-equation estimation for observed-variable models by two-stage
least.squares
1
2
x1
x2
x11
1
x 21
 31
3
4
x3
x4
5
x5
6
x6
x32
 21
 12
2
x42
x53
x63
1
 11
 22
 32
 23
3
1

2
21

y11
y1
1
y21
y2
2
y32
y3
3
y42
y4
4
2
R2WinBUGS: Running WinBUGS and
OpenBUGS from R / S-PLUS
The R2WinBUGS package provides convenient
functions to call WinBUGS from R.
It automatically writes the data and scripts in a
format readable by WinBUGS for processing in
batch mode, which is possible since version 1.4.
After the WinBUGS process has nished, it is
possible either to read the resulting data into R
by the package itself which gives a compact
graphical summary of inference and
convergence diagnostics or to use the facilities
of the coda package for further analyses of the
output.
Others
bio.infer: Maximum Likelihood Method for Predicting
Environmental Conditions from Assemblage Composition
popbio:
Construction and analysis of matrix population models
Construct and analyze projection matrix models from a demography study
of marked individuals classified by age or stage. The package covers
methods described in Matrix Population Models by Caswell (2001) and
Quantitative Conservation Biology by Morris and Doak (2002).
simecol
allows to implement ecological models (ODEs, IBMs, ...) using a
template-like object-oriented stucture. It helps to organize scenarios and may
also be useful for other areas.
demogR: Analysis of age-structured demographic models
Section 3
Making Elegant Graphics in R
Graphical capabilities
100
120
Simple, exploratory graphics are
easy to produce.
0
20
40
60
dist
80
Publication quality graphics can be created.
Several device drivers are available including:
5
15
20
25
speed
8
6
4
10
2
0
Postscript,
5
-2
-10
0
-5
0
-5
X
5
10
-10
Y
Z
On-screen graphics
10
Section 3
Types of graphics functions
High level functions
>-Plot graphs from data
Low level functions
>-Add to/Modify existing graphs
Graphics parameters
>-Control how plots are displayed
Section 3
Some and high level functions
•
•
•
•
•
•
•
•
Low
points()
lines()
abline()
text()
mtext()
title()
legend()
etc
High
•
•
•
•
plot()
hist()
barplot()
boxplot()
Section 3
Graphical parameters
• list of parameters held in par()
• controls how plots are displayed
• values in par() can be over-ridden by
calls to high-level functions
– values are changed temporarily for plot
– no permanent change to par() values
• can be modified directly
Section 3
Plotting single variables: histogram
• First generate data
rrnorm takes n random samples from a
normal distribution
Can specify mean and sd, defaults are
(0,1)
x<-rnorm(100)
• Plot histogram:
Histogram of x
hist(x)
10
5
0
Frequency
15
20
hist is a generic function to plot
histograms of a dataset
-2
-1
0
x
1
2
Section 3
Plotting single variables: boxplot
-1
-2
• Add title
title(main="Boxplot
of x1")
0
1
• boxplot function
boxplot(x1)
2
Boxplot of x1
Section 3
3
0
-1
-2
-3
boxplot(x1,x2,x3)
title(main=“three
boxplots”)
OR
bp3<data.frame(x1,x2,x3)
boxplot(bp3)
1
2
Generate more data
x2<-rnorm(100)
x3<-rnorm(100,2)
4
Plotting single variables: boxplot
1
2
3
Section 3
• barplot(tN,col=rainbo
w(20)
15
10
5
• barplot(tN)
0
• tN <- table(Ni <rpois(100, lambda=5))
20
Plotting single variables: barplot
0
1
2
3
4
5
6
7
8
9
11
12
Section 3
Barplots for grouped data
• Using the VADeaths table from the datasets package:
•
>-library(datasets); data(VADeaths); VADeaths
• Rural Male Rural Female Urban Male Urban Female
50-54
11.7
8.7
15.4
8.4
55-59
18.1
11.7
24.3
13.6
60-64
26.9
20.3
37.0
19.3
65-69
41.0
30.9
54.6
35.1
70-74
66.0
54.3
71.1
50.0
Section 3
Barplots for grouped data
60
50-54
55-59
60-64
65-69
70-74
0
0
10
20
50
30
40
100
50
150
70-74
65-69
60-64
55-59
50-54
70
200
barplot(VADeaths, col=rainbow(5
barplot(VADeaths,
col=rainbow(5),legend=T) legend=T,beside=T)
Rural Male
Rural Female
Urban Male
Urban Female
Rural Male
Rural Female
Urban Male
Urban Female
Section 3
Plotting single variables: scatterplot
function (n=20)
{x1<-rnorm(n)
par(mfrow=c(3,2))
plot(x1,type="p",main="type p:
points")
plot(x1,type="l",main="type l:
lines")
plot(x1,type="b",main="type b:
both")
plot(x1,type="o",main="type o:
overplot")
plot(x1,type="h",main="type h:
histogram")
plot(x1,type="s",main="type s:
steps")
par(mfrow=c(1,1))}
Section 3
Plotting two variables: scatterplot
• Example data: “women” from the datasets package .
Heights vs weights for America women
To get data into workspace: >library(datasets);data(women); attach(women)
140
130
plot(women$height,women$weight)
120
plot(weight~height)
weight
150
160
plot(x=height,y=weight, data=women)
58
60
62
64
66
height
68
70
72
Section 3
Adding to plots: lines
• Straight lines added using “abline” function
abline(v=x) >-draws a vertical line at value x
abline(v=y) >-draws a horizontal line at value y
abline(a=x,b=m)>- draws a line with intercept x and
gradient m
• lines() adds a series of connected lines
– takes vectors of x and y co-ordinates and draws lines
between them
• segments() draws (disjoint) lines segments between
pairs of points
Section 3
Adding to plots: lines
100
80
60
40
20
abline(cars.lm,
col="red“,lty=2)
0
cars.lm<-lm(dist~speed)
dist
library(datasets)
data(cars)
attach(cars)
plot(dist~speed)
120
• “cars” from the datasets package has speed
and stopping distances
• draw scatter plot and add regression line
5
10
15
speed
20
25
Section 3
Adding to plots: lines and text
•
•
•
•
•
•
•
•
•
•
•
•
•
•
function (df=1,p=0.05)
{x<-seq(1,50,0.01)
y<-dchisq(x,df)
plot(x,y,type="l",xlim=c(0,10),ann=F)
title(xlab=substitute(paste(chi^2,"stati
stic")),ylab="Frequency",main=(substitut
e(paste(chi^2,"distribution with",df,"df,
p=",p))))
critx<-qchisq(p,df,lower.tail=F)
abline(v=critx,col="red")
xr<-c(critx,x[x>critx])
yr<-c(dchisq(critx,df),y[x>critx])
polygon(x=c(critx,critx,xr,10),
y=c(0,dchisq(critx,df),yr,0),col="red”)
xr<-round(critx,0)+2
yr<-max(y)/4
rx<-round(critx,3)
text(xr,yr,substitute(paste("Critical
",chi^2,"=",rx)),col="red")}
Section 3
Adding to plots: margin text
• mtext() adds text in the
margins of figures
• Margins are numbered
1:4, clockwise from the
bottom
• Text is printed a
specified number of
lines away from the edge
of the plot area
• Can print in outer
margins if available
Section 3
Adding to plots: legends
• “Legend” function adds labels to plots
• Needs x and y co-ordinates, and a list of labels
• Can be argument to high-level function, or called separately
legend(x=30,y=0.06,
legend=c("Null","Alternative"
,"Threshold"), lty=c(1,1,2),
col=c("blue","black","red"),lw
d=2)
Section 3
Plotting symbols, pch
• 25 default plotting characters
• Controlled by pch parameter
• Add pch = n argument to a high-level function
plots points of type n
• Argument to pch can be a single character, so
pch = “x” will plot letter x
• pch can take vectors as arguments, so
different points will have different symbols
function () {x<-rep(1:5,times=5) ,y<-rep(1:5,each=5)
par(xpd=T),
plot(x,y,type="p",pch=1:25,cex=3,axes=F,ann=F,bty="n",bg="red")
title(main="Data symbols 1:25")
text(x,y,labels=1:25,pos=1,offset=1)}
Section 3
Plotting symbols: size and width
• symbol size controlled by cex
– cex=n plots a figure n times
the default size
• Line width controlled by lwd
– lwd=n gives the width of the
line
function ()
{x<-rep(1:5,times=5),y<-rep(1:5,each=5)
plot(x,y,type="p",pch=3,lwd=x,cex=y,xli
m=c(0,5),ylim=c(0,5),bty="n", ann=F)
title(main="R plotting symbols: size and
width", xlab="Width (lwd)",ylab="Size
(cex)")}
Section 3
Plotting symbols: colour
• Colour of objects controlled by “col”, see
colours() for all ; palette()
function ()
{x<-rep(1:5,times=5)
y<-rep(1:5,each=5)
plot(x,y,type="p",pch=3,lwd=x,ce
x=y,
col=1:25,
xlim=c(0,5),ylim=c(0,5),bty="
n", ann=F)
title(main="R plotting symbols:
size and width",
xlab="Width (lwd)",
ylab="Size (cex)")}
Section 3
Types of lines, lty
function () {x<-1:10
plot(x,y=rep(1,10),ylim=c(1,6.5),
type="l",lty=1,axes=F,ann=F)
title(main="Line Types")
for(i in 2:6){lines(x,y=rep(i,10),lty=i)}
axis(2,at=1:6,tick=F,las=1)
linenames<paste("Type",1:6,":",c("solid","dashed
","dotted","dotdash","longdash","two
dash"))
text(2,(1:6)+0.25,labels=linenames)
box()
}
Section 3
Using lattice graphics
par(mfrow=c(2,2))
hist(D$wg,
main='Histogram',xlab='Weight Gain',
ylab ='Frequency', col=heat.colors(14))
boxplot(wg.7a$wg, wg.8a$wg, wg.9a$wg,
wg.10a$wg, wg.11a$wg, wg.12p$wg,
main='Weight Gain', ylab='Weight Gain
(lbs)',
xlab='Shift', names =
c('7am','8am','9am','10am','11am','12pm
'))
plot(D$metmin,D$wg,main='Met Minutes
vs. Weight Gain', xlab='Mets
(min)',ylab='Weight Gain (lbs)',pch=2)
plot(t1,D2$Intel,type="l",main='Closing
Stock Prices',xlab='Time',ylab='Price
$')
lines(t1,D2$DELL,lty=2)
Section 3
Grid Graphics
Object Oriented - graphical components can be reused and
recombined: library(grid); library(lattice)
•
•
•
•
•
•
•
A standard set of graphical primitives:
grid.rect(...)
grid.lines(...)
grid.polygon(...)
grid.circle(...)
grid.text(...)
library(help = grid) for details
0
0.25
0.5
text
0.75
1
R Package for Elegant Graphics–
Lattice, ggplot2
Ggplot2: R Graphics
Paul Murrell
Ggplot2: Elegant Graphics
for Data Analysis (Use R)
Hadley Wickham
Thanks you !
www.r-project.org
Descargar

An introduction to R