##########
Page 6

demo(image)

demo(colors)

##########
Page 7

library(help="base")

##########
Page 8

a<-available.packages()
length(unique(a[,1])) # we will explain this syntax later

##########
Page 10

X<-3 # X now has the value 3 (a number); it doesn't matter if there are spaces
before or after the <-
X<-"3" # X now has the character 3, which cannot be used in mathematical
operations
X

X<-log(3) # X now has the value the logarithm of 3 to the base e (the default
base)
X<-3+4 # X now has the value 7
X<-sum(1:53) # X is now the sum of the numbers (integers) from 1 to 53, i.e.
1431
X<-X*3 # X is now 1431 multiplied by 3, i.e. 4293
X<-sum(1:53)/3 # X is now 1431 divided by 3, i.e. 477

##########
Page 11

X<-1/0
X

2/0==1/0

X<-"My aunty's name is Jennifer" # X is now a character string, a single letter
is also a string
X

X<-"My aunty's name is Jennifer"
X

(117+117^2+log(54))/(1+log(54)) # logarithms are given in base e as default

u<-(117+117^2+log(54))/(1+log(54))

(117+117^2+log(54))/(1+log(54))->u
u

##########
Page 14

X<-c("p","q","r","s")

for(i in 20:30) {x<-c(x,i)} # creates a list of the integers from 20 to 30

for(i in 20:30) x<-c(x,i)
x

##########
Page 15

x<-1 # the initial (starting) value of x
y<-NULL
while(x < 10){ # calculates x2 for values of x from 1 to 9
y<-c(y,x^2 )
x<-x+1 # incrementing x by one each cycle
}
y

x<-0
repeat{
x<-x+1
if(x>10000) break # the break command exits from the loop
}

##########
Page 16

help("diversity",try.all.packages=TRUE)

names<-c("Adam","Alice","Brian","Barbara","Charles","Carla","Christina","David")

##########
Page 17

length(names)

sex<-c("M","F","M","F","M","F","F","M")

ages<-c(18,19,22,21,20,29,20,19)

heights<-c(167.5,160,170.5,154,179,166,158,160.5)

mean(ages)

mean(heights) # for the arithmetic mean of their heights

max(ages) # the oldest (maximum) age

min(ages) # the youngest (minimum) age

range(ages) # the range of ages

var(ages) # the variance of the ages

quantile(heights) # quantiles, but normally with many more data points

##########
Page 18

pos<-which(names=="Barbara")
pos

ages[pos]

ages[which(names=="Barbara")]

not20<-which(ages !=20)
names[not20]

##########
Page 19

sum(ages)

Male.friends<-length(which(sex=="M"))
Male.friends

Female.friends<-length(which(sex=="F"))
Female.friends

summary(as.factor(sex))

##########
Page 20

plot(NULL,axes=FALSE,xlim=c(0,2),ylim=c(min(heights),max(heights)),xlab="",ylab="",main="My friends' heights")

##########
Page 21

text(c(rep(1,8)),heights,names)

rep(1,8)

plot(NULL,xlim=c(0,2),ylim=c(min(heights),max(heights)),xlab="",ylab="Height (cm)",main="My friends' heights")

plot(NULL,xlim=c(0,2),axes=F,ylim=c(min(heights),max(heights)),xlab="",ylab="Height (cm)",main="My friends' heights")
axis(2,c(150,155,160,165,170,175,180),at=c(150,155,160,165,170,175,180))

##########
Page 22

axis(2,c(150,155,160,165,170,175,180),at=c(150,155,160,165,170,175,180),xpd=NA)
text(rep(1,8),heights,names)

par(mfrow=c(1,2)) # telling the plotting device that we want one row of plots with two columns

##########
Page 23

sex.as.number<-NULL

sex[5]

for(i in 1:8){
if(sex[i]=="M") sex.as.number[i]<-1
if(sex[i]=="F") sex.as.number[i]<-2
} # end of i loop
sex.as.number

##########
Page 24

ifelse(sex[i]=="M",sex.as.number[i]<-1,sex.as.number[i]<-2)

sex.as.number<-sex
sex.as.number<-gsub("M",1,sex.as.number)
sex.as.number<-gsub("F",2,sex.as.number)
sex.as.number

sex.as.number<-as.integer(sex.as.number)
sex.as.number

colours1<-c("blue","deeppink")
text(c(rep(1,8)),heights,names,col=colours1[sex.as.number])

##########
Page 25

for(i in 1:8){
arrows(0.65,heights[i],0,heights[i],code=2,col=colours1[sex.as.number[i]],lwd=2, length=0.07)} # end i loop

pdf("rplot.pdf") # you can call it what you like

graphics.off() # turns off the current open graphics device

dev.new() # new [graphics] device

##########
Page 26

colours() # or colors()

##########
Page 27

data<-list(c(1,2,2,4,1),c("A","A","B","C","A"))
data

##########
Page 28

table(data)

friends<-c("Jack","Piya","Lek","Ali","Peter","Robert","Sue","Emma","Mike")
eyes<-c("grey","brown","brown","brown","blue","brown","brown","blue","blue")
hair<-c("brown","black","black","black","blonde","brown","brown","blonde","blonde")
table(list(c(eyes),c(hair)))

Hylobates_lar<-152.4
Gorilla_gorilla<-c(171,191,48.4,107.9,195.2)
Pan_troglodytes<-c(548,608.8,391.1,258.3,142.9,274.5)
Pan_paniscus<-456.1
Pongo_pygmaeus<-61
Homo_sapiens<-63.5
names<-c(Hylobates_lar,Gorilla_gorilla,Pan_troglodytes,Pongo_pygmaeus,Homo_sapiens)

##########
Page 29

x<-as.table(Hylobates_lar,Gorilla_gorilla,Pan_troglodytes,Pongo_pygmaeus,Homo_sapiens)
x

table(Hylobates_lar,Gorilla_gorilla,Pan_troglodytes,Pongo_pygmaeus,Homo_sapiens)

x<-as.matrix(Hylobates_lar,Gorilla_gorilla,Pan_troglodytes,Pongo_pygmaeus,
Homo_sapiens)
x

samples<-lengths(list(Hylobates_lar,Gorilla_gorilla,Pan_troglodytes,
Pongo_pygmaeus, Homo_sapiens))
samples

L<-max(samples)
sperm_conc<-matrix(nrow=L,ncol=5) # by default, cells of the matrix are
filled with NA
sperm_conc

##########
Page 30

sperm_conc[1:length(Hylobates_lar),1]<-Hylobates_lar
sperm_conc[1:length(Gorilla_gorilla),2]<-Gorilla_gorilla
sperm_conc[1:length(Pan_troglodytes),3]<-Pan_troglodytes
sperm_conc[1:length(Pongo_pygmaeus),4]<-Pongo_pygmaeus
sperm_conc[1:length(Homo_sapiens),5]<-Homo_sapiens
sperm_conc

mean(sperm_conc[,2],na.rm=TRUE)

##########
Page 32

f<-"Adam,male,18,167.5
Alice,female,19,160
Brian,male,22,170.5
Barbara,female,21,154
Charles,male,20,179
Carla,female,29,166
Christina,female,20,158
David,male,19,160.5"
write(f,file="friends.csv") # for more on writing and reading files see Chapter 31.
Fr<-read.csv("friends.csv",header=FALSE)
Fr

##########
Page 33

Ff<-read.csv(textConnection(f),header=FALSE,col.names=c("name","sex",
"age","height"))

colnames(Fr)<-c("name","sex","age","height")

class(Fr$name) # class(Fr$V1) if we have not changed the names

class(Fr$sex)

class(Fr$age)

class(Fr$height)

attach(Fr)
mean(age) # age is the age of the third column in our dataframe F

detach(Fr)

mean(age)

##########
Page 34

with(Fr,mean(age))

age

sites<-c(rep("Benton",7),rep("Warrenton",5))
ttx_resist<-c(0.29,0.77,0.96,0.64,0.7,0.99,0.34,0.17,0.28,0.2,0.2,0.37)

TTXmatrix<-cbind(sites,ttx_resist)
TTXmatrix

##########
Page 35

class(TTXmatrix)

TTXdf<-as.data.frame(TTXmatrix)

class(TTXdf$ttx_resist)

TTXdf$ttx_resist<-as.numeric(TTXdf$ttx_resist)
TTXdf$ttx_resist

TTXdf<-as.data.frame(TTXmatrix) # starting again with our original dataframe
TTXdf$ttx_resist<-as.numeric(as.character(TTXdf$ttx_resist))
class(TTXdf$ttx_resist)

##########
Page 36

a<-matrix(1:6,nrow=3,ncol=2,dimnames=list(c("x","y","z"),c("p","q")))
a

class(a)

is.array(a)

is.table(a)

b<-as.table(a)
b

is.matrix(b)

is.array(b)

is.table(b)

##########
Page 38

mosses<-c(53,85,69,98,130,90)
liverworts<-c(47,67,85,113,190,239)

plot(mosses)

##########
Page 39

plot(NULL,axes=F,ylab="",xlab="",main="",xlim=c(0,6),ylim=c(0,10))
text(rep(0,10),c(0:9),c(0:9),cex=2)
points(rep(1,10),c(0:9),pch=c(0:9),cex=2)
text(rep(2,10),c(0:9),c(10:19),cex=2)
points(rep(3,10),c(0:9),pch=c(10:19),cex=2)
text(rep(4,6),c(0:5),c(20:25),cex=2)
points(rep(5,6),c(0:5),pch=c(20:25),cex=2)

plot(mosses,ylab="Number of records",xlab="Site number",pch=15,col="blue",
main ="Altitudinal distribution of bryophytes")

##########
Page 40

points(liverworts,pch=16,col="brown")

##########
Page 41

plot(mosses,ylab="Number of records",xlab="Site number",pch=15,col="blue",
main ="Altitudinal distribution of bryophytes",ylim=c(0,250))
points(liverworts,pch=16,col="brown")

plot(mosses,ylab="Number of records",xlab="Site number",pch=15,col="blue",main=
"Altitudinal distribution of bryophytes",ylim=c(min(mosses,liverworts),
max(mosses,liverworts)))
points(liverworts,pch=16,col="brown") # figure will be same as the previous

##########
Page 42

plot(c(1:6,1:6),c(mosses,liverworts),pch=c(rep(15,6),rep(16,6)),col=
c(rep("blue",6), rep("brown",6)),ylab="Number of records",xlab="Site
number",main ="Altitudinal distribution of bryophytes") # figure will be
same as previous

lines(c(1:6),mosses,col="blue",lty=1) # lty= specifies the line type: 1 for
normal, 2–6 for different types of dashes

lines(c(1:6),liverworts,col="brown",lty=2)

altitude<-c(400,600,800,1000,1200,1300)

plot(c(altitude,altitude),c(mosses,liverworts),pch=c(rep(15,6),rep(16,6)),
col=c(rep("blue",6),rep("brown",6)),ylab="Number of records",xlab=
"Altitude a.s.l. (m)",main ="Altitudinal distribution of bryophytes")
lines(altitude,mosses,col="blue",lty=1)
lines(altitude,liverworts,col="brown",lty=2)

##########
Page 43

A<-"Liverworts"
B<-"Mosses"
legend("topright",pch=c(15,16),col=c("blue","brown"),c(A,B))

legend("topleft",pch=c(15,16),col=c("blue","brown"),c(A,B))

##########
Page 44

legend("topleft",pch=c(15,16),col=c("blue","brown"),c(A,B),cex=0.8,
inset=0.05) # inset is the proportion of the plot area to inset the legend by

##########
Page 45

mosses<-c(53,85,69,98,130)

mosses<-c(53,85,69,98,130,90)

length(mosses)

length(liverworts)

a<-c(1:5); b<-c(1:5); c<-c(1:3)
rbind(a,b,c)

##########
Page 46

mosses<-c(53,85,69,NA,130,90)
liverworts<-c(47,67,85,113,190,239)

plot(altitude,mosses,pch=15,col="blue",ylab="Number of bryophytes",xlab="Altitude a.s.l. (m)",main ="Altitudinal distribution of bryophytes",ylim=c(min(na.omit (mosses,liverworts)),max(mosses,liverworts, na.rm=TRUE))) # see Information Box 6.4
points(altitude,liverworts,pch=16,col="brown")
lines(altitude,mosses,col="blue",lty=1)
lines(altitude,liverworts,col="brown",lty=2)

a<-c(1,2,3,4,5,6,NA,7,8)
mean(a)

mean(a,na.rm=TRUE)

mean(na.omit(a))

is.na(a)

##########
Page 47

mosses<-c(53,85,69,98,130,90); liverworts<-c(47,67,85,113,190,239)
bryos<-c(mosses,liverworts) # concatenating the two into a single vector
alt<-c(altitude,altitude) # concatenating the matched altitudes

use_colours<-c(rep("blue",6),rep("brown",6))
use_symbols<-c(rep(15,6),rep(16,6))

##########
Page 48

moss_species<-c(12,19,20,24,27,19)
liverwort_species<-c(17,19,23,31,42,46)

par(mar=c(5.1,4.1,4.1,4.1))
plot(alt,bryos,pch=use_symbols,col=use_colours,ylab="Number of records",xlab="Altitude a.s.l. (m)",main ="Altitudinal distribution of bryophytes",ylim=c(min(bryos),max(bryos)))
lines(altitude,mosses,col="blue",lty=1)
lines(altitude,liverworts,col="brown",lty=2)
legend("topleft",pch=c(15,16),col=c("blue","brown"),c(A,B),cex=0.8,inset=0.05)

par(new=TRUE) # tells R to plot the new graph on top of the existing one
plot(NULL,axes=F,ylab="",xlab="",main="",xlim=c(400,1300),ylim=c(min(c(moss_species,liverwort_species)),max(c(moss_species,liverwort_species))))
points(altitude,moss_species,pch=3,col="blue")
lines(altitude,moss_species,lty=3,col="blue") # using a different dotted line style
points(altitude,liverwort_species,pch=4,col="red")
lines(altitude,liverwort_species,lty=4,col="red") # using a different dotted line style
axis(4,c(10,20,30,40,50),at=c(10,20,30,40,50))
axis(4,las=0) # las=0 means print text vertically, las=1 means horizontally
mtext(side=4,line=2.5,"Numbers of species") # using mtext (=margin text) to print axis name on right (‘side=4’) and 2.5 lines from the axis

##########
Page 49

d<-read.csv("Wright_y_Muller_Landau_2006.csv",header=TRUE)
head(d)

##########
Page 50

colnames(d)

range(d[,1])

mean(d[,1]); median(d[,1])

par(mfrow=c(2,2))

levels(d[,3])

as.numeric(d[,3])

levels(as.numeric(d[,3]))

plot(d[,1],d[,2],pch=as.numeric(d[,3]),col=as.numeric(d[,3]),xlab=colnames(d)[1],ylab=colnames(d)[2])

plot(log10(d[,1]),d[,2],pch=as.numeric(d[,3])+14,col=as.numeric(d[,3])+3,type="p",xlab="Log(rural popn density)",ylab=colnames(d)[2])

##########
Page 51

plot(d[,1],d[,2],log="x",pch=as.numeric(d[,3])+14,col=as.numeric(d[,3])+3,axes=FALSE,xlim=c(0.1,1000),ylim=c(0,100),xlab="Rural popn density",ylab="% potential forest remaining")
axis(1,c("0.1","1.0","10","100","1000"),at=c(10^-1,10^0,10^1,10^2,10^3),cex.
axis=0.8)
axis(2,at=c(0:10)*10,cex.axis=0.8)

legend("bottomleft",pch=as.numeric(levels(as.factor(as.numeric(d[,3]))))+14,col=as.numeric(levels(as.factor(as.numeric(d[,3]))))+3,levels((d[,3])),bty="n")

as.factor(as.numeric(d[,3]))

levels(as.factor(as.numeric(d[,3])))

as.numeric(levels(as.factor(as.numeric(d[,3]))))+14

install.packages("sfsmisc")
library(sfsmisc)
eaxis(1,at=c(10^-1,10^0,10^1,10^2,10^3),cex.axis=0.8)
axis(2,at=c(0:10)*10,cex.axis=0.8)

##########
Page 53

data<-rbind(mosses,liverworts) # combining the original vectors into a matrix
barplot(data,beside=TRUE) # using the ‘beside’ argument to put bars side
by side

ants<-c("Oecophylla smaragdina","Pheidole spp.","Monomorium destructor","Paratechina longicornis","Odontoponera denticulata","Tapinoma melanocephalum","Monomorium pharoensis")
unburnt<-c(2957,2020,556,445,418,401,308)
burnt<-c(696,1425,445,1450,1256,269,476)
data<-cbind(unburnt,burnt)
barplot(data,ylab="Numbers of individuals",sub="Dipterocarp forest condition")

##########
Page 54

par(mfrow=c(2,2))
bardensity<-c(0,10,15,20,40,30,100)
barangle<-c(45,135,0,90,45,130,90)
barplot(data,ylab="Numbers of individuals",xlab="Dipterocarp forest condition",density=bardensity,angle=barangle,lwd=2,col="black")
barplot(data,ylab="Numbers of individuals",xlab="Dipterocarp forest condition",sub="effects on ant abundance",col=c("black","red","blue","green","yellow","grey50","white")) # adding a subtitle to x axis using the ‘sub=’ argument

##########
Page 55

par(mfrow=c(1,1))
barplot(data,ylim=c(0,9500),ylab="Numbers of individuals",
xlab="Dipterocarp forest condition",col=c("black","red","blue","green", "yellow","grey50","white"))
legend ("top",ants,fill=c("black","red","blue","green","yellow","grey50","white"),cex=0.7)

legend("top",ants,fill=c("black","red","blue","green","yellow","grey50","white"),cex=0.7, text.font=3)

barplot(data,ylab="Numbers of individuals",xlab="Dipterocarp forest condition",xlim=c(0,4.2),col=c("black","red","blue","green","yellow","grey50", "white"))
legend("right",ants,fill=c("black","red","blue","green","yellow","grey50","white"),cex=0.7,text.font=3)

##########
Page 56

legend(0,8100,pch=c(15,16),col=c("blue","brown"),legend=expression(italic("A. scientific name"),("other ants")),xpd=NA) # legend positionspecified here by x, y coordinates

##########
Page 57

r1<-c(23,5,9,0)
r2<-c(7,0,3,0)
r3<-c(0,0,12,21)
r4<-c(165,276,170,122)
r5<-c(0,0,7,12)

##########
Page 58

bird.names<-c("GH","WH","PH","BH")
hornbill.forage<-rbind(r1,r2,r3,r4,r5)
hornbill.forage

colnames(hornbill.forage)<-bird.names
hornbill.forage

rownames(hornbill.forage)<-c("cracking","probing","hawking","plucking",
"snatching")
hornbill.forage

hornbill.forage<-t(hornbill.forage)
hornbill.forage

barplot(hornbill.forage, beside=TRUE)

##########
Page 59

barplot(hornbill.forage,beside=TRUE,legend.text=bird.names,args.legend=list(x="topleft",bty="n")) # the legend position argument is called ‘x’

##########
Page 60

image<-c("poacher","domestic_dog","ranger","villager","tourist","barking_deer","wild_pig","macaque","porcupine","sambar_deer","gaur","elephant","serow","pangolin","Indian_civet","black_bear","sun_bear","mongoose","dhole","clouded_leopard","golden_cat","binturong")
number_of_images<-c(25,17,16,7,5,60,60,37,30,28,13,9,3,2,37,21,16,11,8,8,5,2)

par(mfrow=c(2,2))
barplot(number_of_images,col="black",names=image,horiz=TRUE)

barplot(number_of_images,col="black",names=image,horiz=TRUE,las=2)

##########
Page 61

par(mar=c(5.1,9.1,4.1,2.1))
barplot(number_of_images,col="black",names=image,horiz=TRUE,las=2)

number_of_images<-rev(number_of_images)
image<-rev(image)
barplot(number_of_images,col="black",names=image,horiz=TRUE,las=2,xlab="Number of camera trap photos",ylab="‘Animal’ photographed")

text(xpd=NA,-32,12,"‘Animal’ photographed",srt=90) # the argument srt is used to give a y axis label with 90 degree text rotation

barplot(number_of_images,col="black",names=image,horiz=TRUE,las=2,xlab="Number of camera trap photos",ylab="")
title(ylab="‘Animal' photographed",line=7,xpd=NA) # within the title call we use ‘line’ to move the positions of the x labels

##########
Page 63

error.bars<-function(xv,z,nn){xv<-barplot(yv,ylim=c(0,(max(yv)+max(z))),names=nn,ylab=deparse(substitute(yv)))
g=(max(xv)-min(xv))/50 # the 50 determines the horizontal lengths of the bar at top and bottom of each error bar
for(i in 1:length(xv)){
lines(c(xv[i],xv[i]),c(yv[i]+z[i],yv[i]-z[i])) # vertical part of the line
lines(c(xv[i]-g,xv[i]+g),c(yv[i]+z[i],yv[i]+z[i]))
lines(c(xv[i]-g,xv[i]+g),c(yv[i]-z[i],yv[i]-z[i]))
} # end i loop
} # end function

flies<-c("L. dux","C. megacephala","C. rufifacies","C. nigripes","L. cuprina","M. domestica")
males<-c(6032,4371,5367,4801,3678,3488)
females<-c(6073,5635,5207,4778,3620,3427)
standarddev_males<-c(385,344,385,323,155,110)*1.96 # data in the original paper are 1 S.D.
standarddev_females<-c(207,514,244,277,209,231)*1.96

data<-matrix(c(males,females),2,6,byrow=TRUE) # the 2 and 6 are the number of rows and columns in the matrix respectively
errors<-matrix(c(standarddev_males,standarddev_females),2,6,byrow=TRUE)
bardata<-barplot(data,beside=TRUE,ylim=c(0,max(data)+max(errors)), ylab="Mean number ommatidea (+/- 95% confidence)",names=flies,cex.names=0.7,xlab="Fly species",col=c("lightblue","pink"))
arrows(bardata,data+errors,bardata,data-errors,angle=90,code=3,length=0.1)
legend("topright",fill=c("lightblue","pink"),c("males","females"))

##########
Page 64

bardata<-barplot(data,beside=TRUE,ylim=c(0,max(data)+max(errors)),ylab="Mean number ommatidea (\u00B1 95% confidence)",names=flies,cex.names=0.7, xlab="Fly species",col=c("lightblue","pink"))

##########
Page 65

TWIG<-c("Andaman_Islands","Cambodia","Laos","Myanmar","Nicobar_Islands","South_China_Sea","Thailand","Vietnam","Borneo","Christmas_Island", "Cocos_Keeling_Island","Jawa","Lesser_Sunda_Islands","Malaya","Maluku","Philippines","Sulawesi","Sumatra","Bismarck_Archipelago","New_Guinea", "Solomon_Islands")
number_of_genera<-c(453,686,767,1304,277,39,1605,1401,1396,172,62,1403,812,1457, 898,1406,1000,1300,509,1424,665)
number_of_genera_per_area<-c(943,598,622,812,779,291,1097,1061,841,980,438,1307,802,1354,1026,1062,850,877,612,844,861)

##########
Page 67

par(mar=c(2,4.1,4.1,2.1))
par(mfrow=c(1,2))
pie(number_of_genera,labels=as.character(TWIG),cex=0.5)
pie(number_of_genera_per_area,labels=as.character(TWIG),cex=0.5)

install.packages("circlize")
library(circlize)

num_genera<-c(2,27,70,117,113,84,75,81,52,47,26,18,7,4)
log_body_mass_base2<-c(1:14)
n<-length(num_genera)
m<-cumsum(num_genera) # cumsum gives the cumulative sum
x<-360*m/sum(num_genera)
x<-c(0,x) # adjacent pairs from cumsum are the beginning and end angles except for the first one so we append a zero
plot(NULL,xlim=c(–8,15),ylim=c(–12,12),axes=F,ann=F, type="n",xlab="",ylab="") # set up a blank plot area with size based on max mass (log base 2)
for(i in 1:n){
draw.sector(start.degree=x[i],end.degree=x[i+1],rou1=0,rou2=log_body_mass_base2[i],center=c(0,0),clock.wise=FALSE,col=gray(1-(0.1+(0.9*num_genera[i]/max(num_genera)))),border="black",lwd=1,lty=0)}

##########
Page 68

for(i in 1:n){
draw.sector(start.degree=x[i],end.degree=x[i+1],rou1=0,rou2=log_body_mass_base2[i],center=c(0,0),clock.wise=FALSE,col=heat.colors(100, 0.1+(0.9*num_genera[i]/max(num_genera))),border="black",lwd=1,lty =0)}

for(i in 1:n){
draw.sector(start.degree = x[i],end.degree=x[i+1],rou1= 0,rou2= log_body_mass_base2[i],center=c(0,0),clock.wise=FALSE,col=rainbow(256,0.9) [round(256*(0.1+(0.9*num_genera[i]/max(num_
genera))))],border="black",lwd=1,lty=0)}

aggregate(count ~ spray,data=InsectSprays,sum)

tapply(InsectSprays$count, InsectSprays$spray, sum)

##########
Page 69

Tilapia<-5
Mystus_singaringan<-c(1,2,2,2)
Osteochilus_vittatus<-c(26,30,11,2,28,6,19,12,4,4,18,13,29,21,9,8,21,19,28)
Henicorhynchus_siamensis<-c(46,22,33,28,50,42,31,4,13,8,11,8,18)
Barbodes_gonionotus<-c(8,11,66)
Anabas_testudineus<-c(12,26)
Trichogaster_microlepis<-c(9,5)
Trichogaster_pectoralis<-c(21,12,6,3)

counts<-c(Tilapia,Mystus_singaringan,Osteochilus_vittatus,
Henicorhynchus_siamensis,Barbodes_gonionotus,Anabas_testudineus,
Trichogaster_microlepis,Trichogaster_pectoralis)

sample_sizes<-c(130.0,28.95,27.2,17.92,26.12,81.65,78.31,59.46,62.84,
58.28,45.03,35.02,46.46,43.54,24.64,31.52,57.05,58.87,59.7,39.02,
28.64,65.03,43.28,36.2,64.57,67.43,70.42,89.3,75.44,51.38,78.92,30.69,
21.69,21.84,23.45,18.21,49.42,110,105.63,56.25,51.55,91.42,32.78,
39.83,45.76,62.03,34.32,14.77) # in grams of tissue analysed

species<-c("Tilapia",rep("Mystus singaringan",4),rep("Osteochilus vittatus",19), rep("Henicorhynchus siamensis",13),rep("Barbodes gonionotus",3),rep("Anabas testudineus",2),rep("Trichogaster microlepis",2),rep("Trichogaster pectoralis",4))

##########
Page 70

counts_per_gram<-counts/sample_sizes

boxplot(counts_per_gram ~ factor(species),ylab="Metacercaria per gram")

##########
page 71

boxplot(counts_per_gram ~ factor(species),ylab="Metacercaria per
gram",las=2)

par(mar=c(13,4,4,2))
boxplot(counts_per_gram ~ factor(species),ylab="Metacercaria per gram",las=2,cex=0.8)

##########
page 72

rownames(as.matrix(sort(tapply(counts_per_gram,species,mean))))

Fish<-rownames(as.matrix(sort(tapply(counts_per_gram,species,mean))))

boxplot(counts_per_gram ~ factor(species,Fish),ylab="Metacercaria per gram",las=2,cex=0.8)

##########
Page 73

boxplot(counts_per_gram ~ factor(species,Fish),notch=TRUE,ylab="Metacercaria per gram",las=2,cex=0.8)

##########
Page 74

par(mar=c(5.1, 4.1, 4.1, 2.1))
plot(TukeyHSD(aov(counts_per_gram ~ factor(species))))

##########
Page 75

levels(as.factor(species))

LETTERS[1:length(levels(as.factor(species)))]

##########
Page 76

species<-factor(species,labels=LETTERS[1:length(levels(as.factor(species)))])

x<-factor(species,labels=LETTERS[1:length(levels(as.factor(species)))])
plot(TukeyHSD(aov(counts_per_gram ~ factor(x))))

TukeyHSD(aov(counts_per_gram ~ factor(x)))[1]

##########
Page 78

Narum<-c(0.0001,0.0010,0.0062,0.0101,0.0214,0.0227,0.0273,0.0292,0.0311,0.0323, 0.0441,0.0490,0.0573,0.1262,0.5794)

p.adjust(Narum,"bonferroni"); length(which(p.adjust(Narum,"bonferroni")<0.05))

p.adjust(Narum,"hochberg"); length(which(p.adjust(Narum,"hochberg")<0.05))

p.adjust(Narum,"BH"); length(which(p.adjust(Narum,"BH")<0.05))

##########
Page 79

install.packages("ggplot2")
library(ggplot2)

##########
Page 80

head(Loblolly)

attach(Loblolly)
plot(age,height,col=as.factor(Seed))

##########
Page 81

ggplot(Loblolly,aes(x=age,y=height,color=as.factor(Seed))) + geom_point()

##########
Page 82

detach(Loblolly)

remove.packages("ggplot2")

install.packages("ggpubr")
library("ggpubr")

x<-c(0.14990444,0.14352947,0.18572451,0.22692724,0.32947406,0.39478415,0.43702197,0.44205266,0.45779166,0.4772302,0.4990459,0.5148496,0.54387033, 0.5666838,0.61728454,0.6657855,0.72059345) # ants in fish diet

y<-c(0.074910976,0.6254242,0.3017297,0.24208963,0.14281006,1.0448253,0.62397754,1.1371646,1.2069718,0.8333741,1.0302496,
0.95308137,0.70407337,0.62443894,1.1078949,0.39062706,1.2676938) # ornamentation of food dummy

g<-as.data.frame(cbind(x,y))
ggscatter(g,x="x",y="y",xlab="ants in diet",ylab="similarity",add="reg.line",conf.int=TRUE,cor.coef=TRUE,cor.method=c("spearman"))

cor.test(g[,1],g[,2],method=c("spearman"))

##########
Page 85

install.packages("BiocManager")
library(limma)

khao<-readLines("KhaoYai_woodpeckers.txt")

##########
Page 86

head(khao)

ky<-"Blythipicus pyrrhotis, Bay woodpecker
Chrysocolaptes guttacristatus, Greater flameback
Chrysophlegma flavinucha, Greater yellownape
Dendrocopos analis, Freckle-breasted woodpecker
Dendrocopos atratus, Stripe-breasted woodpecker
Dinopium javanense, Common flameback
Dryobates cathpharius, Crimson-breasted woodpecker
Dryocopus javensis, White-bellied woodpecker
Jynx torquilla, Eurasian wryneck
Picumnus innominatus, Speckled piculet
Picus chlorolophus, Lesser yellownape
Picus erythropygius, Black-headed woodpecker
Picus guerini, Black-naped woodpecker
Sasia ochracea, White-browed piculet
Yungipicus canicapillus, Grey-capped pygmy woodpecker"
ky

##########
Page 87

khao<-unlist(strsplit(ky,"\n"))
khao

allpeckers<-unique(c(khao,doi,kaeng))

length(allpeckers)

allpeckers<-union(khao,union(kaeng,doi))

in_all_three_parks<-intersect(doi,intersect(kaeng,khao))
in_all_three_parks

##########
Page 88

data<-matrix(NA,nrow=length(allpeckers),ncol=3)
colnames(data)<-c("kaeng","khao","doi") # assigning column names
data<-as.data.frame(data) # as.data.frame coerces the matrix into the dataframe structure
for(i in 1: length(allpeckers)){ # the next three lines use ‘ifelse’ to determine whether each woodpecker species occurs in each of the three national parks
ifelse (allpeckers[i] %in% kaeng,data$kaeng[i]<-1,data$kaeng[i]<-0)
ifelse (allpeckers[i] %in% khao,data$khao[i]<-1,data$khao[i]<-0)
ifelse (allpeckers[i] %in% doi,data$doi[i]<-1,data$doi[i]<-0)}
data

##########
Page 89

apply(data,1,sum) # the 1 tells apply to work on each row

apply(data,2,sum) # the 2 tells apply to work on each column

m2<-vennCounts(data)

vennDiagram(m2,names=c("Kaeng Krachan","Khao Yai","Doi Inthanon"),main= "Woodpecker species",circle.col=c("coral","paleturquoise3","bisque3"),lwd=3) # we livened this up a bit by colouring the circles and making lines thicker
mtext ("Aves: Picidae")

##########
Page 90

X<-structure(vd1,class="VennCounts")

plot(NULL,xlim=c(0,2),ylim=c(0,2),axes=F,xlab="",ylab="",main="Woodpecker
overlap between three Thai N.P.s")

cx1<-0.7
cx2<-1.3
cx3<-1
cy1<-0.8
cy2<-0.8
cy3<-0.8+(((cx2-cx1)*3^0.5)/2) # the height of an equilateral triangle is √3 / 2 × its side

##########
Page 91

theta<-seq(0,2*pi,0.01)
radius<-0.5
x<-radius*cos(theta) # obtain the circle’s x and y coordinates using trig functions
y<-radius*sin(theta)
lines(x+cx1,y+cy1,lwd=3,col="blue")
lines(x+cx3,y+cy3,lwd=3,col="darkgreen")
lines(x+cx2,y+cy2,lwd=3,col="red")
text(0.7,0.22,"Kaeng Krachan N.P.",col="blue")
text(0.45,1.8,"Khao Yai N.P.",col="darkgreen")
text(cy3,0.22,"Doi Inthanon N.P.",col="red")

symbols(0.7,0.8,circles=10,fg="blue",cex=20,lwd=4)

##########
Page 92

all<-intersect(intersect(kaeng,khao),doi)
text(1,1.04,length(all))

kk_ky_only<-length(intersect(kaeng,khao))-length(all)
kk_di_only<-length(intersect(kaeng,doi))-length(all)
di_ky_only<-length(intersect(doi,khao))-length(all)
text(1,0.65,kk_di_only)
text(0.75,1.1,kk_ky_only)
text(1.25,1.1,di_ky_only)

kk<-length(kaeng)-length(all)-kk_ky_only-kk_di_only
ky<-length(khao)-length(all)-di_ky_only-kk_ky_only

##########
Page 93

di<-length(doi)-length(all)-di_ky_only-kk_di_only
text(0.6,0.65,kk)
text(1.4,0.65,di)
text(1,1.4,ky)

##########
Page 96

number_of_genera_per_area<-c(943,598,622,812,779,291,1097,1061,841,980,438,1307,802,1354,1026,1062,850,877,612,844,861)human_footprint<-c(34.90,21.95,20.74,19.91,22.58,0.62,24.60,28.12,20.27,33.66,3.48,41.05,31.53,26.24,23.88,33.79,30.68,28.99,26.00,16.72,29.53)
plot(human_footprint,number_of_genera_per_area,xlab="Mean human footprint value", ylab="Number of genera standardized by area",pch=15)

##########
Page 105

chisq.test(c(80,105,107,96,90,122)) # the default is that all scores are equally likely, H0

TAO<-c(63,14,44)
PAOD<-c(375,235,440)

freq_TAO<-121/1050 # = 11.5% of PAOD admissions.

expected_TAO<-freq_TAO*PAOD
expected_TAO

##########
Page 106

chisq.test(TAO,p=expected_TAO,rescale.p=TRUE) # if you don’t include ‘rescale.p=TRUE’ the probabilities must add to be precisely 1

chisq.test(TAO,p=PAOD,rescale.p=TRUE)

chisq.test(TAO,p=PAOD/sum(PAOD)) # the probabilities add up precisely to 1

##########
Page 107

males<-c(782383,23580702)
females<-c(391841,23828880)
clonorchis<-rbind(males,females) # rbind sticks the rows together
clonorchis

clonorchis<-matrix(c(782383,391841,23580702,23828880),nrow=2)
clonorchis

chisq.test(clonorchis)

##########
Page 108

p<-seq(0,1,0.001)
plot(NULL,xlim=c(0,1),ylim=c(0,1),xlab="p",ylab="Frequency")
lines(p,p^2,col="red",lwd=2)
lines(p,(1-p)^2,col="blue",lwd=2)
lines(p,2*p*(1-p),col="purple",lwd=2)
text(0.2,0.95,expression('q'^2*' Propn(aa)'),col="blue") # use expression and place text to be superscripted between ^and*
text(0.8,0.95,"q^2 Prop(AA)",col="red")
text(0.5,0.68,"2pq Prop(Aa)",col="purple")
arrows(0.2,0.93,0.13,0.7569,code=2,col="blue",length=0.1) # 0.7569 is
our x value^2
arrows(0.8,0.93,0.87,0.7569,code=2,col="red",length=0.1)
arrows(0.5,0.66,0.5,0.5,code=2,col="purple",length=0.1)

##########
Page 110

obs<-c(1469,138,5)
Tot<-2*(1469+138+5)
fA<-(2*1469+138)/Tot

fa<-(2*5+138)/Tot
fa+fA # to check that the allele frequencies do add up to 1

expected<-c(fA^2,2*fA*fa,fa^2) # these are expected phenotype proportions

chisq.test(obs,p=expected)

install.packages("DescTools")
library(DescTools)

##########
Page 111

GTest(clonorchis)

##########
Page 112

frog_males<-c(395,390,410,390,395,455,400,370,460,490)
frog_females<-c(335,360,315,335,405,265,315,390,405,380)
t.test(frog_males,frog_females)

boxplot(cbind(frog_males,frog_females),col="greenyellow")

##########
Page 113

Chinese<-c(0.875,0.781,0.75,0.677,0.655,0.5,0.4,0.333,0.167)
Thai<-c(0.818,0.435,0,0.771,0.829,0.5,0.677,0.429,0.6)
herbalists<-rbind(Chinese,Thai)
t.test(Chinese,Thai,paired=TRUE)

##########
Page 114

var.test(Chinese,Thai)

without_ants_beetle_damage <-c(0.44,0.49,0.59,0.70,0.73,0.94,1.03,1.09,1.42,2.10)
with_ants_beetle_damage<-c(0.19,0.39,0.03,0.04,0.19,0.36,0.42,0.20,0.21,0.09,0.05,0.16,0.07,0.05,0.06,0.06,0.01,0.14)
boxplot(without_ants_beetle_damage,with_ants_beetle_damage,ylab="Percent leaf damage",names=c("without","with"),col="olivedrab")
text(1.5,-0.45,"Presence of tree ant nest(s)",cex=1.1,xpd=NA) # xpd=NA is included because we want to write a text label outside of the plotting area

##########
Page 115

par(mfrow=c(2,2))
hist(without_ants_beetle_damage)
hist(with_ants_beetle_damage)

##########
Page 116

shapiro.test (with_ants_beetle_damage)

ks.test(with_ants_beetle_damage,"pnorm")

install.packages("nortest")
library(nortest)
ad.test(with_ants_beetle_damage)

##########
Page 117

wilcox.test(without_ants_beetle_damage,with_ants_beetle_damage)

##########
Page 118

binom.test(4,28)

##########
Page 121

CB<-read.csv("CarribTot.csv")
par(mfrow=c(2,2))
plot(CB)
plot(log(CB,10),xlab="log(Area)",ylab="log(Total species)")
mod_CB<-lm(log(CB$Tot_butterflies,10) ~ log(CB$Area_square_km,10))
summary(mod_CB)

##########
Page 122

abline(mod_CB,col="blue",lwd=2)

##########
Page 123

a<-simulate(mod_CB,10)
plot(CB,pch=15)
for(i in 1:10) points(cbind(CB[1],10^a[i]),pch=3,col="red") # because 'modelt' was based on log10 we have to back convert by taking 10 to the power of the model values

nests<-c(0,0,0,0,0,0,0,0,0,0,1,1,2,2,2,2,2,4,5,5,7,7,9,11,15,20,24,30)
total_damage_percent<-c(3.72,3.11,2.15,1.40,1.12,0.94,0.55,0.59,0.73,0.83,0.19,0.44,0.19,0.64,0.66,0.68,1.03,0.19,0.20,0.13,0.10,0.59,0.09,0.05,0.06,0.06,0.04,0.14)
leafbeetle_damage_percent<-c(0.44,0.49,0.59,0.70,0.73,0.94,1.03,1.09,1.42,2.10,0.19,0.39,0.03,0.04,0.19,0.36,0.42,0.20,0.21,0.09,0.05,0.16,0.07,0.05,0.06,0.06,0.01,0.14)
plot(nests,total_damage_percent,xlab="Number of weaver ant nests",ylab="Percentage leaf damage",pch=1)
points(nests,leafbeetle_damage_percent,pch=4)

##########
Page 125

modelA<-lm(leafbeetle_damage_percent ~ nests)

##########
Page 126

par(mfrow=c(2,2))
plot(modelA)

##########
Page 127

modelB<-lm(leafbeetle_damage_percent ~ log(nests+1))
plot(modelB)

##########
Page 128

plot(log(nests+1),total_damage_percent^0.333,xlab="ln(weaver ant nests+1)",ylab= "Percentage leaf damage ^0.33",pch=1,col="red")
points(log(nests+1),leafbeetle_damage_percent^0.333,pch=4,col="green")
modelC<-lm(total_damage_percent^(1/3) ~ log(nests + 1))
summary(modelC)

modelT<-lm(leafbeetle_damage_percent^(1/3) ~ log(nests+1))
summary(modelT)

##########
Page 129

abline(lm(total_damage_percent^(1/3) ~ log(nests+1)),lwd=2,col="red")
abline(modelT,lwd=2,col="green")

modelAnts<-lm(g$y ~ g$x)
install.packages("car")
library(car)
qqPlot(modelAnts)

##########
Page 131

Ochlero<-read.table("culicid_DT.txt",header=TRUE)
head(Ochlero)

with(Ochlero,plot(Temperature_C,ln.Development_time.,col=localities,
pch=as.numeric(localities)+7)) # using with to avoid typing "Ochlero" every time

##########
Page 133

Sar<-subset(Ochlero,localities=="Sarmiento")
BAires<-subset(Ochlero,localities=="BuenosAires")

regrSar<-lm(log(log(Sar$ln.Development_time.)) ~ Sar$Temperature_C)
regrBA<-lm(log(log(BAires$ln.Development_time.)) ~ BAires$Temperature_C)

summary(regrSar)[[4]][2]; summary(regrBA)[[4]][2]

summary(regrSar)[[4]][4]; summary(regrBA)[[4]][4]

Z<-abs(summary(regrSar)[[4]][2]-summary(regrBA)[[4]][2])/sqrt(summary(regrSar)[[4]][4]^2+summary(regrBA)[[4]][4]^2)
Z

qt(0.975,83)

##########
Page 134

age<-c(6.96,6.96,7.94,5.89,7.88,10.90,13.83,14.89,14.84,16.87,17.94,20.81,22.93, 24.86,25.83)
length<-c(12.78,13.54,14.12,28.15,26.93,18.70,27.46,28.22,38.0,29.28,28.72,35.39, 36.74,48.15,38.28)
plot(age,length,xlim=c(0,25),ylim=c(0,50),ylab="Length (cm)",xlab="Age(months)")

##########
Page 135

modelG<-nls(length ~ C*(1-exp(-a*age)),start=c(C=100,a=0.5))
summary(modelG)

xv<-seq(0,25,0.1)

yv<-predict(modelG,list(age=xv)); lines(xv,yv,lwd=2,col="blue")

##########
Page 138

pairs(ChickWeight, panel=panel.smooth,lwd=3) # panel=panel.smooth fits a line to the data

##########
Page 139

class(T)

class(T)="logical"

##########
Page 140

knor<-read.table("krill.txt",header=TRUE)
head(knor)

par(mfrow=c(2,2))
plot(knor$wwt_grams,knor$ulO2_min,ylab="Oxygen consumption (ul/min)",xlab="Wet weight(grams)")
plot(log10(knor$wwt_grams),log10(knor$ulO2_min),ylab="log(Oxygen consumption (ul/min))",xlab="log(Wet weight(grams))")

par(mfrow=c(2,2))

##########
Page 141

model1<-lm(log(knor$ulO2_min) ~ log(knor$wwt_grams))
plot(model1)

summary(model1)

##########
Page 142

model3I<-lm(log(knor$ulO2_min) ~ knor$wwt_grams+I(knor$wwt_grams^2)+I(knor$wwt_grams^3))

modelpoly3<-lm(log(knor$ulO2_min) ~ poly(log(knor$wwt_grams),3))

summary(model3I)

##########
Page 143

summary(modelpoly3) # the model with orthogonal variables

##########
Page 144

model2<-lm(log(knor$ulO2_min) ~ poly(log(knor$wwt_grams),2))
anova(model2,modelpoly3)

##########
Page 145

AIC(model2,model3I)

##########
Page 149

d<-read.csv("fledge_number.csv")
par(mfrow=c(2,2))
plot(d,xlab="Clutch initiation date",ylab="Number of fledglings",pch=17,col="blue")

freq<-table(d[,2])
freq

barplot(freq,xlab="Number of fledglings",ylab="Frequency")

mean(d[,2])

var(d[,2])

##########
Page 150

install.packages("vcd")
library(vcd)
gf<-goodfit(freq,type="poisson",method="ML")
plot(gf,type="standing",xlab="Number of fledged chicks",ylab="Frequency",scale="raw")
plot(gf,type="hanging",xlab="Number of fledged chicks",ylab="Frequency",scale="raw")

model_card<-glm(d$fledged ~ d$date,family="poisson")
summary(model_card)

##########
Page 151

PV<-read.csv("Passiflora_visits.csv",header=TRUE)
head(PV)

plot(PV,pch=16,col="red")

##########
Page 152

mean(PV$Total_visitors); var(PV$Total_visitors)

modelQP<-glm(PV$Total_visitors ~ PV$Flower_diam_mm,family="quasipoisson")
summary(modelQP)

##########
Page 153

modellm<-lm(PV$Total_visitors ~ PV$Flower_diam_mm)
summary(modellm)

##########
Page 155

data<-InsectSprays
oneway.test(InsectSprays$count ~ InsectSprays$spray)

##########
Page 156

modelSprays<-lm(InsectSprays$count ~ InsectSprays$spray)
summary(modelSprays)

sprayMeans<-tapply(InsectSprays$count,InsectSprays$spray,mean)
sprayMeans

aovSprays<- aov(InsectSprays$count ~ InsectSprays$spray)
TukeyHSD(aovSprays)

##########
Page 158

tj<-read.table("Tunjai_regeneration.txt",header=TRUE)
head(tj)

colnames(tj)[9]<-"PercentEstab"

arcsine_establishment<-asin(sqrt(tj$PercentEstab/100)) # division by 100 converts percentages to proportions and we will bind the proportions onto the right-hand end of our dataframe

##########
Page 160

tj<-cbind(tj,arcsine_establishment) # the vector name will automatically be assigned as the column name

levels(tj$Size_categ) # seed size (intermediate, large and small)

levels(tj$Coat) # seed coat thickness (medium, thick and thin)

levels(tj$Shape) # shape (flat, oval and round)

boxplot(tj$arcsine_establishment ~ tj$Size_categ)

##########
Page 161

is.ordered(tj$Size_categ) # commands that start with ‘is.’ are asking whether or not the following vector has the given property, in this case, is it ordered?

tj$Size_categ<-ordered(tj$Size_categ,levels=c("S","I","L"))
tj$Coat<-ordered(tj$Coat,levels=c("Tn","M","Tk"))

is.ordered(tj$Size_categ)

tj$Size_categ<-factor(tj$Size_categ,levels=c("S","I","L"),ordered=TRUE)
tj$Size_categ

is.ordered(tj$Size_categ)

boxplot(tj$arcsine_establishment ~ tj$Size_categ)

##########
Page 162

modelA<-aov(tj$arcsine_establishment ~ tj$Size_categ)

par(mfrow=c(2,2)) # two rows and two columns of graphs to see all four plots at the same time rather than having to step through them one by one

plot(modelA)

##########
Page 163

summary(modelA)

modelB<-aov(tj$PercentEstab ~ tj$Size_categ)
summary(modelB)

##########
Page 164

tj$logit_establishment<-with(tj, log((PercentEstab/100)/(1-(PercentEstab/100))))

modelC<-aov(tj$logit_establishment ~ tj$Size_categ)
plot(modelC,which=2) # the which argument allows us to select which plots to view

##########
Page 165

summary(modelC)

##########
Page 166

Goby_growth<-read.csv("Chapter_14_Malone_gobies.csv")
with(Goby_growth,plot(Initial_length.mm,Growth.mm,pch=c(rep(1,22),rep(19,24)),cex=3,col=c(rep("black",22),rep("grey50",24))))
legend("topright",c("Visual implant","Acrylic paint"),pch=c(1,19),bty="n",col=c("black","grey50"),cex=1.2)

##########
Page 167

par(mfrow=c(2,2))
goby_model<-lm(Goby_growth$Growth.mm~Goby_growth$Treatment*Goby_growth$Initial_length.mm)
plot(goby_model)

##########
Page 168

summary(goby_model)

##########
page 169

sir<-read.csv("CHAP14_Sirindhornia_count.csv")
modelPoiss<-glm(sir$pollinaria ~ sir$flowers*sir$species-1,family="poisson")
summary(modelPoiss)

##########
Page 172

psyllids<-c(0:9)
no_Bracon<-c(5,139,153,109,70,28,20,2,3,1)
Bracon<-c(63,57,61,36,20,9,3,0,0,0)

PropBracon<-Bracon/(Bracon+no_Bracon)
plot(psyllids,PropBracon,pch=15,xlab="Number of psyllids in gall",ylab="Proportion of galls containing a Bracon")

##########
Page 173

data<-cbind(Bracon,no_Bracon)
model_galls<-glm(data~psyllids,family=binomial) # we can omit "family="
summary(model_galls)

##########
Page 174

xv<-seq(0,9,0.01)
yv<-predict(model_galls,list(psyllids =xv),type="response")
lines(xv,yv,lwd=2,col=4)

nPsyllid<-NULL
Bracon_status<-NULL

##########
Page 175

for(i in 1:10){
np<-i-1 # because the first row in data corresponds to zero psyllids
if(data[i,1]>=1){
for(j in 1:data[i,1]){
nPsyllid<-c(nPsyllid,np)
Bracon_status<-c(Bracon_status,0)}}
if(data[i,2]>=1){
for(j in 1:data[i,2]){
nPsyllid<-c(nPsyllid,np)
Bracon_status<-c(Bracon_status,1)
}} # end if statement and j loop
} # end i loop
xpandedData<-cbind(nPsyllid,Bracon_status)

plot(jitter(xpandedData[,1],1,0.15),jitter(xpandedData[,2],1,0.03),
ylab="Bracon in gall",xlab="Total number of psyllids in gall",main="",
pch=20,cex=0.3,axes=F)
axis(1,c(0:9),at=c(0:9)); axis(2,c("absent","present"),at=c(0,1))
Total number

##########
Page 177

mantid<-c(2,8,10,12,3,5,2,0,15,6,2,2,1,3,7,16,5,17,9,8,3,8,4,6,9,4,10,10,17,14,32,6,12)
flower<-c(4,2,4,7,3,5,4,2,8,10,1,12,0,1,6,11,4,9,9,10,2,4,8,4,4,7,4,8,4,5,28,5,5)
control<-c(0,0,0,0,2,2,0,0,0,4,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1,0,1,0,1,0,2,0,0)
m<-as.data.frame(cbind(as.numeric(mantid),as.numeric(flower),as.numeric(control)))
colnames(m)<-c("Total.Mantid","Total.Flower","Total.Control")
head(m)

m<-read.table("mantid.txt",header=TRUE)[,2:4]

colMeans(m[,1:3])

rowSums(m)

##########
Page 178

apply(m,1,sum) # the 1 means that the function sum will be applied to the rows

apply(m,2,sum) # the 2 applies it to the columns

par(mfrow=c(2,3)) # Creating a plot area so we can put three histograms side-by- side and not have them too tall
par(mar=c(2,2,2,2)) # setting the margin around each plot to 2, reducing the border around the plots so they look neater
brks<-seq(1:34)-1.5 # creates a sequence from –0.5 to 32.5, which spans the full range of data with breaks always between other integer counts
hist(m$Total.Mantid,breaks=brks,main="Visits to mantid")
abline(v=mean(m$Total.Mantid),col="blue",lwd=2,lty=2)
hist(m$Total.Flower,breaks=brks,main="Visits to flowers")
abline(v=mean(m$Total.Flower),col="blue",lwd=2,lty=2)
hist(m$Total.Control,breaks=brks,main="Visits to controls")
abline(v=mean(m$Total.Control),col="blue",lwd=2,lty=2)

##########
Page 179

reps<-nrow(m)
type_list<-c(rep("Mantid",reps),rep("Flower",reps))
visitors<-c(m$Total.Mantid,m$Total.Flower)
# new.data<-as.data.frame(cbind(type_list,visitors))
model1<-glm(visitors ~ type_list,family=poisson)
par(mfrow=c(2,2))
plot(model1)

##########
Page 181

model2<-glm(visitors ~ type_list,family=quasipoisson)# plot(model2)
summary(model2)

##########
Page 182

wilcox.test(mantid,flower,paired=T)

library(MASS)
data(snails)
spA<-snails[1:48,]; spB<-snails[48:96,]
par(mfrow=c(2,2))
boxplot(spA$Deaths~spA$Temp,main="Species A",xlab="Temperature")
boxplot(spA$Deaths~spA$Rel.Hum,main="Species A",xlab="Relative Humidity")
boxplot(spB$Deaths~spB$Temp,main="Species B",xlab="Temperature")
boxplot(spB$Deaths~spB$Rel.Hum,main="Species B",xlab="Relative Humidity")

##########
Page 183

snails$deadORalive<-cbind(snails$Deaths,20-snails$Deaths)

##########
Page 184

head(snails$deadORalive)

m1<-glm(deadORalive ~ Species/(Temp+Rel.Hum),family=binomial,
data=snails)

m2<-glm(deadORalive~Species/(Temp+Rel.Hum),family=quasibinomial,data=snails)
summary(m1)

summary(m2)

##########
page 187

Random<-NULL # this will become our vector of 10,000 random numbers
seed<-9997 # less than 10,000 and each seed will lead to a different random number sequence
A<-16877 # this number is 7^5

##########
Page 188

M<-39119 # a large prime number; 2147483647 = 2^31-1 is often used
for(i in 1:10000){
seed<-M*A*seed %% M/M # %% means the modulus, i.e. the remainder after division; M*A*seed is likely too large for the function round to be used; %/% gives the integer part of a division
s<-as.character(seed)
Random<-c(Random,as.numeric(substr(s,nchar(s),nchar(s)))) # gets the last digit of seed
} # end i loop
hist(Random,breaks=c(0:10)-0.5) # plots the frequencies of values in bins from minus 0.5 to 9.5; 0 goes in the 1st bin, 1 in the 2nd, etc.

##########
Page 189

cerana<-c(1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,0,1,1,0,1,1,0,0,0,
1,0,1,0,0,0,0,1, 1,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,
1,1,1,1,0,1,0,0,1,1,0,0,0,1,0,1, 1,0,0,1,1,1,0,1,1,1,1,0,1,1,1)
dorsata<-c(0,1,1,1,0,0,1,1,0,0,0,0,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,
1,0,0,1,1,1,1,0, 1,0,1,1,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,
1,0,0,0,0,0,0,0,1,0,0,0,0,1,0, 0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0)
florea<-c(0,0,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,1,0,1,1,0,0,0,1,
1,1,0,0,0,1,0,0,0,0,1,0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,0,0,0,1,1,1,1,1,1,1,
0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,1, 0,1,0,0,1,1,0,0,1,1,1,0,0,1,0)
mellifera<-c(0,1,0,1,0,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,
1,0,0,0,0,1,1,1,1,0,0,1, 1,0,0,1,0,1,1,1,0,1,0,0,1,1,1,0,0,1,0,0,1,0,1,1,1,1,1,1,
0,0,1,1,1,0,0,0,1,1,1,1,1, 1,1,1,1,1,0,0,0,0,1,0,0,0,0,1,1)
data<-cbind(cerana,dorsata,florea,mellifera)

raw_score<-sum(rowSums(data[,1:4])^2)
raw_score

n<-nrow(data) # getting the number of plant species
randomizations<-NULL
for(i in 1:100000){
d1<-sample(data[,1],n,replace=FALSE) # randomizing the 1s and 0s for each bee species
d2<-sample(data[,2],n,replace=FALSE)
d3<-sample(data[,3],n,replace=FALSE)
d4<-sample(data[,4],n,replace=FALSE) # technically we only really need to shuffle the values for three of the four bee species because once three have been randomized, the fourth will be random with respect to each of these anyway – but it doesn’t matter if we do
shuffled_score<-sum((d1+d2+d3+d4)^2)
randomizations<-c(randomizations,shuffled_score)
} # end i loop

##########
Page 190

a<-1:5
b<-LETTERS[1:5]
c<-c("cat","dog")
paste(c(a,b,c))

paste(c(a,b,c),collapse="") # nothing to separate the elements

paste(c(a,b,c),collapse=" ") # a space to separate the elements, etc.

paste("The value of p was ",round(1/240,4))

x<-c("a", "b", "c", "d")
y<-c("w", "x", "y", "z")
paste(x,y)

paste(x,y,sep="*")

paste(x,y,sep="*",collapse=" ")

hist(randomizations)
abline(v=raw_score,lwd=2,col=2)
text(528,12000,"observed",srt=270,col="red")
equals<-length(which(randomizations==raw_score))
below<-length(which(randomizations<raw_score))+0.5*equals
above<-length(which(randomizations>raw_score))+0.5*equals
text(480,16000,paste("less than = ",below)) # x positions will have to be adjusted because of output from randomization
text(560,16000,paste("greater than = ",above),xpd=NA)
text(560,12000,paste("P = ",round(above/(above+below),4)),col="blue")

##########
Page 192

treatments<-sample(rep(c("A","B","C","D"),16))

noquote(matrix(treatments,nrow=8))

m<-as.data.frame(matrix(rep(c("A","B","C","D"),16),nrow=8))
m

##########
Page 193

m2<-noquote(sapply(m[,1:8],sample))
m2

##########
Page 195

Lap1<-c(0,0,0,0,0,0,0,0.025,0.030); Lap2<-c(0,0,0,0,0,0,0,0,0.030); Lap3<-c(0,0,0,0,0,0,0,0.350,0); Lap4<-c(0,0.214,0.250,0,0.140,0,0.024,0.475,0.156); Lap5<-c(0.167,0.119,0.286,0.278,0.160,0.136,0.119,0,0.345); Lap6<-c(0.833,0.667, 0.214,0.611,0.580,0.705,0.666,0,0.345); Lap7<-c(0,0,0.250,0.111,0.120,0.091,0.167,0.150,0.094); Lap8<-c(0,0,0,0,0,0.068,0.024,0,0)
Mpi1<-c(0,0.024,0,0,0.040,0.023,0.050,0,0); Mpi2<-c(0.167,0,0.107,0,0.060,0.068,0.025,0,0); Mpi3<-c(0.833,0.262,0.107,0.056,0.360,0.045,0.125,0.263,0.500); Mpi4<-c(0,0.547,0.607,0.721,0.460,0.546,0.650,0.737,0.406); Mpi5<-c(0,0.143, 0.143,0.167,0.080,0.273,0.150,0,0.094); Mpi6<-c(0,0.024,0.036,0.056,0,0.045,0,0,0)
Pgi1<-c(0,0,0.036,0,0.040,0,0.190,0.265,0.056); Pgi2<-c(0,0,0,0,0,0.023,0,0,0); Pgi3<-c(0.667,0.095,0.214,0.111,0.140,0.931,0.525,0.059,0.167); Pgi4<-c(0.333,0.429,0.214,0.278,0.240,0.023,0.190,0.529,0.527); Pgi5<-c(0,0.309,0.429,0.444, 0.480,0.023,0.095,0.147,0.250); Pgi6<-c(0,0.167,0.107,0.167,0.100,0,0,0,0)

Allozyme_data<-matrix(NA,nrow=9,ncol=20) # create a blank matrix of appropriate
dimensions

Allozyme_data[,1]<-Lap1; Allozyme_data[,2]<-Lap2; Allozyme_data[,3]<-Lap3; Allozyme_data[,4]<-Lap4; Allozyme_data[,5]<-Lap5; Allozyme_data[,6]<-Lap6; Allozyme_data[,7]<-Lap7; Allozyme_data[,8]<-Lap8; Allozyme_data[,9]<-Mpi1;Allozyme_data[,10]<-Mpi2; Allozyme_data[,11]<-Mpi3; Allozyme_data[,12]<-Mpi4; Allozyme_data[,13]<-Mpi5; Allozyme_data[,14]<-Mpi6; Allozyme_data[,15]<-Pgi1; Allozyme_data[,16]<-Pgi2; Allozyme_data[,17]<-Pgi3; Allozyme_data[,18]<-Pgi4; Allozyme_data[,19]<-Pgi5; Allozyme_data[,20]<-Pgi6

rownames<-c("Trat","Chon_Buri","Chumphon_5","Chumphon_6","Surat_Thani","Satun", "Trang","Sydney_Aust","Magnetic_Island")
colnames(Allozyme_data)<-c("Lap1","Lap2","Lap3","Lap4","Lap5","Lap6", "Lap7","Lap8","Mpi1","Mpi2","Mpi3","Mpi4","Mpi5","Mpi6","Pgi1","Pgi2","Pgi3","Pgi4","Pgi5","Pgi6")

pca1<-prcomp(t(Allozyme_data),scale.=TRUE)
head(pca1$rotation,3) # just to check

##########
Page 196

colours1<-c("orange","orange","orange","orange","orange","orange","orange","blue","blue")

par(mfrow=c(2,2))
plot(pca1$rotation[,1],pca1$rotation[,2],col=colours1,pch=15)
plot(pca1$rotation[,2],pca1$rotation[,3],col=colours1,pch=15)
plot(pca1$rotation[,1],pca1$rotation[,3],col=colours1,pch=15)
plot(pca1$rotation[,2],pca1$rotation[,4],col=colours1,pch=15)

##########
Page 197

localities<-c("Trat","Chon_Buri","Chumphon_5","Chumphon_6","Surat_Thani", "Satun","Trang","Sydney_Aust","Magnetic_Island")
plot(pca1$rotation[,1],pca1$rotation[,2],col=colours1,pch=15)
xoffset<-0.01*(max(pca1$rotation[,1])-min(pca1$rotation[,1]))
yoffset<-0.02*(max(pca1$rotation[,2])-min(pca1$rotation[,2]))
text(pca1$rotation[,1]+xoffset,pca1$rotation[,2]+yoffset,localities,cex=0.7,xpd=NA)

data(iris)
par(mfrow=c(2,2))
PCA1<-prcomp(iris[,1:4],center=TRUE,scale.=TRUE)
biplot(PCA1)
PCA_log<-prcomp(log(iris[,1:4]),center=TRUE,scale.=TRUE)
biplot(PCA_log)

##########
Page 198

library(ggplot2)
pca2_plot<-autoplot(PCA1,data=iris,colour='Species')
pca2_plot

##########
Page 201

PNG_butterflies<-c(1,1,56,53,1,1,2,8,4,1,15,7,2,22,6,9,1,16,28,49,2,5,562,30,63,7,1922,321,231,59,4,1,1,6,2,2,5,2,3,4,3,163,3,61,2,58,6,5,15,22,85,50,8,6,97,3,140,37,4,3,17,2,1,41,1,1,45,3,4,45,14,5,17,11,1,2,1,17,48,24,100,2,4,2,52,17,1,1,1,12,4,9,9,21,11,4,7,37,51,7,9,66,1,4,8,51,39,2,19,15,1,3,16,75,10,141,26,21,2,3,689,15,9,9,3,2,6,3,149,1,5,79,1,1,2,18,154,10,113,3,3,3,1,1,142,2,9,1,9,16)
par(mfrow=c(2,2)) # setting the plot area to be 2 × 2 graphs
hist(PNG_butterflies,col="grey50")

length(which(PNG_butterflies==1))

s<-seq(0:11)-1 # annoyingly seq assumes that the first value is 1 so we have to take one off to start at 0
s2<-2^s
s3<-s2-0.5
hist(PNG_butterflies,breaks=s3,col="grey50")

log_butterflies<-log(PNG_butterflies,2)
hist(log_butterflies,col="grey50")

install.packages("vegan")
library(vegan) # use install.packages("vegan") if you have not already installed it
fisherfit(PNG_butterflies)

plot(fisherfit(PNG_butterflies))

##########
Page 202

data<-read.table(file="PNG_transects.txt",header=TRUE)
head(data)

##########
Page 203

PA<-replace(PA<-data,PA>1,1)

species_recorded<-NULL
hold_temp<-NULL
for(i in 1:ncol(PA)){ # stepping through the columns
temp<-grep(1,PA[,i]) # finding which rows in the column have a specimen, i.e. a 1
hold_temp<-c(hold_temp,temp) # adding current column ‘temp’
species_recorded<-c(species_recorded,length(unique(hold_temp)))
# removing duplicates to get total number of species recorded up to and including current column
}
plot(species_recorded,xlab="Transect number",ylab="Total number of species")

##########
Page 205

N<-100
M<-matrix(NA,nrow=N,ncol=ncol(PA)) # to hold the results of N (100) random sequences
for(i in 1:N){
to_add_seq<-sample(1:50,50,replace=FALSE)
species_recorded<-NULL
hold_temp<-NULL
for(j in to_add_seq){ # automatically steps through all values in the list to_add_seq
temp<-grep(1,PA[,j])
hold_temp<-c(hold_temp,temp)
species_recorded<-c(species_recorded,
length(unique(hold_temp))) }
M[i,]<-species_recorded} # end i loop
head(M)

##########
Page 206

plot(NULL,xlim=c(1,50),ylim=c(1,93),xlab="Transect number",ylab="Total number of species")
for(i in 1:100) points(c(1:50),M[i,],pch=15,cex=0.2)

points(c(1:50),colMeans(M),col="red",pch=16)

##########
Page 207

lines(smooth.spline(c(1:50),colMeans(M),spar=0.4,all.knots=F),lwd=2,lty=1,col="blue")

##########
Page 208

propPNGbutts<-PNG_butterflies/sum(PNG_butterflies)

Simpson<-1-sum(propPNGbutts^2)
Simpson

ShannonH<--sum(propPNGbutts*log(propPNGbutts,2))
ShannonH

##########
Page 209

species_observed<-NULL; Simpson_transect<-NULL; Shannon_transect<-NULL
for(i in 1:length(data)){
temp<-data[which(data[,i]>0),i]
species_observed<-c(species_observed,length(temp))
Simpson_transect<-c(Simpson_transect,1-sum((temp/sum(temp))^2))
Shannon_transect<-c(Shannon_transect,-sum(sum((temp/sum(temp)) * log(temp/sum(temp),2))))}
par(mfrow=c(2,2))
plot(species_observed,Simpson_transect,xlab="Number of butterflies spp in transect",ylab="Simpson's Index",pch=18,col="cyan3")
plot(species_observed,Shannon_transect,xlab="Number of butterflies spp in transect",ylab="Shannon's Index",pch=18,col="chocolate")

plot(Shannon_transect,Simpson_transect,xlab="Shannon",ylab="Simpson",col="red")

##########
Page 210

install.packages("vegan")# if not already installed
library(vegan)
SW<-diversity(propPNGbutts,index="shannon",base=exp(1))
SW

ShannonH<--sum(propPNGbutts*log(propPNGbutts)) # log(propPNGbutts,2.718281828459)); also note the minus sign before sum
ShannonH

##########
Page 211

cocc<-read.table("British_ladybirds.txt",header=TRUE)
head(cocc)

##########
Page 212

Number_of_species<-nrow(cocc)
Number_of_species

Number_of_genera<-length(levels(cocc$Genus))
Number_of_genera

species_per_genus<-as.data.frame(rowSums(table(cocc)))
head(species_per_genus)

a<-sort(sample(c(2:Number_of_species),Number_of_genera-1)) # the default function sample has replace=FALSE, which is what we want; the function sort arranges elements into ascending order
c<-a-0.5
c

hold<-floor(c[1]) # setting the number of species in the first genus to the smallest random number generated
for(i in 1:length(c)-1) hold<-c(hold,c[i+1]-c[i])
hold<-c(hold,Number_of_species-floor(c[i+1])) # here we are adding the last species to the end of the last genus

##########
Page 213

length(hold)==Number_of_genera

sum(hold)==Number_of_species

par(mfrow=c(2,2)) # we will be plotting the results of using two different approaches
plot(NULL,ylim=c(0,100),xlim=c(-0,Number_of_species+1),axes=FALSE,ylab="",xlab="")
rug(0:53,ticksize=0.03,side=1,lwd=0.5) # the first 53 small tick marks
rug(c(5,10,15,20,25,30,35,40,45,50),ticksize=0.05,side=1,lwd=2,line=0.5,lend="square") # second bigger and fatter tick marks every five species
twocolours<-c(rep(c("red","yellow"),1+Number_of_genera/2)) # set up a list of alternating colours
remember<-NULL # 'remember' will store all our results for analysis at the end
for(j in 1:100){ # start 100 random replicates loop
a<-sort(sample(c(2:Number_of_species),Number_of_genera-1))
c<-a-0.5
hold<-floor(c[1])
for(i in 1:length(c)-1){ # a loop nested within the j loop
hold<-c(hold,c[i+1]-c[i])
} # end i loop
hold<-c(hold,Number_of_species-floor(c[i+1]))# Now to plot each replicate
m<-0
for(i in 1:Number_of_genera){
lines(c(m,cumsum(hold)[i]),c(j,j),lwd=3.5,col=twocolours[i],lend="square")
m<-cumsum(hold)[i]
} # end i loop
remember<-rbind(remember,hold)
} # end of the 100 randomizations j loop

##########
Page 214

Styracaceae<-c(3,2,1,2,3,1,1,2,2,5,8,87) # numbers of spp. in same order as genus list
Styracaceae<-sort(Styracaceae)
num_genera<-length(Styracaceae)
nsp<-sum(Styracaceae)
list3<-c(1:(nsp-1)) # note the list stops one short of the total number of species
mem<-NULL

##########
Page 215

for(i in 1:100){
breaks<-sample(list3,(num_genera-1)) # creating break points that are integers between 1 and num_species minus 1
breaks<-sort(breaks)
# if a random break happens to be number 1, the vector ‘d’ has a zero added to its left-hand end so that one minus zero will give the number of species in the first random genus
d<-c(min(list3)-1,breaks) 
# if a random break point happens to be the number of species minus one, the last length of the broken stick must just include the last species. This is allowed for in ‘e’ by adding max(list3)+1 (i.e. nsp).
e<-c(breaks,max(list3)+1)
r<-sort(e-d)
mem<-rbind(mem,r)} # end of i loop

plot(1:num_genera,Styracaceae,xlab="Rank genus based on no. spp.",ylab="Number of species",ylim=c(0,nsp),col="blue",pch=7)
for(i in 1:100) points(1:num_genera,mem[i,],col="red",pch=4)

##########
Page 216

random_is_bigger<-NULL
for(i in 1:num_genera){
random_is_bigger<-c(random_is_bigger,length(which
(mem[,i]> Styracaceae[i]))/100)} # end i loop
random_is_bigger

temp<-t(apply(mem,1,function(x){ifelse(x>Styracaceae,1,0)}))
prop_out_of_100<-colSums(temp)/100
prop_out_of_100c

install.packages("devtools")
devtools::install_github("hafen/hbgd")
devtools::install_github("stefvanbuuren/brokenstick")

##########
Page 218-219

Input=("
day grass_cover oil_pad
1 16 28
2 16 28
3 16 28
4 16 28
5 16 28
6 16 28
7 16 28
8 16 27
9 16 26
10 16 25
11 16 24
12 16 24
13 16 24
14 16 24
15 16 24
16 16 23
17 16 21
18 16 20
19 15 19
20 15 17
21 15 17
22 14 14
23 14 12
24 14 11
25 14 8
26 14 8
27 12 8
28 9 NA
")
killdeer<-read.table(textConnection(Input),header=TRUE)

par(mfrow=c(2,2))
plot(killdeer$day,killdeer$grass_cover,ylim=c(0,28),pch=16,col="darkgreen", xlab="Survival (days)",ylab="Number of nests")
points(killdeer$day,killdeer$oil_pad,pch=16,col="brown") # now plotting as proportion surviving
plot(killdeer$day,killdeer$grass_cover/max(killdeer$grass_cover),ylim =c(0,1),pch=16,col="darkgreen",xlab="Survival (days)",ylab="Number of nests")

points(killdeer$day,killdeer$oil_pad/max(killdeer$oil_pad,na.rm=TRUE),pch=18, col="brown")
par(mfrow=c(1,1))

##########
Page 220

install.packages("survival")

library(survival)

grass_survive<-killdeer[nrow(killdeer),2] # the number of nests remaining at
the end

##########
Page 221

data<-c(rep(killdeer[nrow(killdeer),1],grass_survive),rep(0,grass_survive),
rep(1,grass_survive))
grass_survival<-matrix(data,ncol=3,byrow=FALSE) # creating matrix with just surviving nests
z<-length(levels(as.factor(killdeer[,2])))-1 # minus 1 because we have already filled the matrix with the first level
temp<-grass_survive
for(i in 1:z){
d<-max(which(killdeer[,2]>temp))
diff<-killdeer[d,2]-temp
for(j in 1:diff) grass_survival<-rbind(grass_survival,c(d,1,1))
temp<-killdeer[d,2]}
colnames(grass_survival)<-c("time","status","habitat")

oil<-Filter(Negate(anyNA),killdeer$oil_pad)
oil

# unlist(lapply(killdeer$oil_pad,function(x) x[!is.na(x)]))
# killdeer$oil_pad[sapply(killdeer$oil_pad,Negate(anyNA))]

oil_survive<-oil[length(oil)]
data<-c(rep(killdeer[length(oil),1],oil_survive),rep(0,oil_survive),rep(2,oil_survive))
oil_survival<-matrix(data,ncol=3,byrow=FALSE)
temp<-oil_survive
for(i in 1:13){
d<-max(which(killdeer[,3]>temp))
diff<-killdeer[d,3]-temp
for(j in 1:diff) oil_survival<-rbind(oil_survival,c(d,1,2))
temp<-killdeer[d,3]} # end i loop
killdeerSurv<-as.data.frame(rbind(grass_survival,oil_survival))
killdeerSurv # these are the data in the format required by the function Surv

##########
Page 222

par(mfrow=c(1,1))
kdfit<-survfit(Surv(killdeerSurv$time,killdeerSurv$status) ~ killdeerSurv$habitat)
plot(kdfit,xlab="Days",ylab="Nest survivorship",lwd=2,lty=c(1,2),col=c("darkgreen", "brown"))

##########
Page 223

summary(kdfit)

##########
Page 224

kdgrass<-Surv(killdeerSurv$time[1:16],killdeerSurv$status[1:16])~killdeerSurv$habitat[1:16]
kdoil<-Surv(killdeerSurv$time[17:44], killdeerSurv$status[17:44])~killdeerSurv$habitat[17:44]
g95<-survfit(kdgrass)
o95<-survfit(kdoil)
plot(g95,col="darkgreen",xlab="Days",ylab="Nest survivorship",main="Kaplan-Meier estimate with 95% confidence bounds",lwd=2)
lines(o95,col="brown",lwd=2)

##########
page 225

kdc<-coxph(Surv(killdeerSurv$time,killdeerSurv$status) ~ killdeerSurv$habitat)
summary(kdc)

##########
Page 226

ggsurvplot(kdfit, data=killdeerSurv, pval=TRUE, pval.method=TRUE, conf.int=TRUE,palette=c("springgreen","orangered1"))

ggsurvplot(kdfit, data=killdeerSurv, fun="event", linetype="strata", pval=TRUE, conf.int=TRUE,pval.method=TRUE, palette="Brewer")

##########
Page 227

date()

Sys.Date()

format(Sys.Date(),"%d %b %Y")

format(Sys.Date(),"%d %m %Y")

format(Sys.Date(),"%d %b %y")

format(Sys.Date(),"%D")

##########
Page 228

Sys.time()

format(Sys.time(),"%d %b %Y")

date<-as.POSIXlt(Sys.time(),format="%Y/%m/%d")
date

Sys.timezone(location=TRUE)

substring(Sys.time(),3,4)

date<-as.POSIXlt(Sys.Date(),format="%Y/%m/%d")
date

julian(date)

tmp<-as.Date("31Aug19",format="%d%b%y")
format(tmp,"%j")

CorsBT<-read.csv("Corsican_blue_tit.csv")
CorsBT


##########
Page 229

dat<-as.Date("27March1986",format="%d%B%Y")
dat

dat<-as.Date("27-March-1986",format="%d-%B-%Y")
dat

as.Date(CorsBT[2,2],format="%d %B %Y")

CorsBT[,1:2]<-apply(CorsBT[,1:2],2,function(x){x<-as.Date(x,format="%d %B %Y")})
CorsBT

CorsBT$Building=as.Date(CorsBT$Building,origin='1970-01-01')
CorsBT$Laying=as.Date(CorsBT$Laying,origin='1970-01-01')

##########
Page 230

CorsBT

class(CorsBT$Building)

start_each_year<-c("1986-01-01", "1987-01-01", "1987-01-01", "1988-01-01", "1988-01-01", "1988-01-01", "1986-01-01", "1987-01-01", "1987-01-01", "1988-01-01", "1988-01-01", "1988-01-01")
building<-NULL # a new variable with lower case "b"
laying<-NULL # ditto
for(i in 1:nrow(CorsBT)){
building<-c(building,CorsBT$Building[i]-as.Date(start_each_year[i]))
laying<-c(laying,CorsBT$Laying[i]-as.Date(start_each_year[i]))
}
building

laying

plot(NULL,xlab="Pair number",ylab="Nest building to laying period", xlim=c(1,length(building)),ylim=c(min(c(building,laying)),c(max(c(building,laying)))))
linecolours<-c(rep("darkgreen",6),rep("orange",6))
for(i in 1:length(building)) lines (c(i,i),c(building[i],laying[i]),lwd=4,col=linecolours[i])
text(3,136,"Mainland",col="darkgreen",cex=2)
text(10,86,"Corsica",col="orange",cex=2)

##########
Page 231

substring(start_each_year,1,4)

years<-format(as.Date(start_each_year,format="%Y-%m-%d"),"%Y")
years

for(i in 1:length(building)) text(i,laying[i]+3,years[i],cex=0.5)

CorsBT$Laying-CorsBT$Building

##########
Page 232

hope<-read.csv("HOPE_Date_Order_Table_1.csv",header=TRUE)
head(hope,12)

##########
Page 233

burial<-as.Date(hope$Burial.date,"%m/%d/%y")
head(burial,14)

burial[1:11]<-sub("20","19",burial[1:11])
head(burial,12)

data<-data.frame(burial,hope$Sex,hope$Age)
head(data)

year<-format(as.Date(data$burial,format="%d/%m/%Y"),"%Y")
head(year)

data$Year<-as.numeric(year)

##########
Page 234

obs1<-readLines("European_corn_borer_2003.txt")
head(obs1)

obs2<-scan("European_corn_borer_2003.txt")

head(obs2,40)

##########
Page 235

par(mfrow=c(1,2))
hist(obs2,breaks=max(obs2)-min(obs2)+1) # breaks is set to the total number of days covered
plot(density(obs2))
par(mfrow=c(1,1))

##########
Page 236

MD<-c(31,28,31,30,31,30,31,31,30,31,30,31)
sum(MD) # just to check

FD<-c(1,1+cumsum(MD)[-length(MD)])
FD

months_in_range<-which(FD>=min(obs2) & FD<=max(obs2))
months_in_range

plot(density(obs2))

MN<-month.name[1:12]
MN

MN[months_in_range] # or more directly month.name[months_in_range]

FD[months_in_range]

rug(FD[months_in_range],lwd=2,lend="square")

##########
Page 237

FirstOfMonth<-paste("01",MN[months_in_range]) # paste automatically puts a space between the items being pasted unless we use collapse=""
FirstOfMonth

text(FD[months_in_range],rep(0.0005,length(months_in_range)),FirstOfMonth,cex=0.9, col="blue")

text(mean(range(obs2)),-0.0032,"Julian date in 2003",xpd=NA,cex=1.2) # mean(range(x)) gives the midpoint of the x axis

##########
Page 238

day_in_year<-c(0.9023162,1.9475052,1.9111508,17.947441,17.61407,19.799334,33.986637, 33.615097,39.078434,49.624847,50.02184,58.024895,65.629875,65.61569,76.569275,81.67198,82.0228,95.510284,97.67119,98.09835,114.81702,113.7471,113.76201,
133.79655,129.08539,145.89676,145.54774,152.80917,161.92287,171.76254,177.9617,191.43936,194.37025,209.6631,209.6791,228.60011,225.7012,247.55856,251.5441,254.8538,270.53964,270.5327,286.1364,286.54538,289.83145,302.56824,302.5624,308.39548,318.93423,318.5678,327.66803,334.57318,334.956,346.61084,350.22012,350.6073)
log_number<-c(1.434,0.454,-0.294,-0.123,0.513,0.506,0.109,-0.041,-0.056,-0.415, 0.258,0.034,-0.887,-1.179,-0.655,-0.596,-0.872,-0.678,-1.187,0.108,-0.671,-0.199,0.107,0.099,0.601,1.73,2.044,1.55,1.69,1.782,1.916,1.909,2.253,2.118,2.447,2.012,2.327,2.35,1.907,2.55,3.01,2.86,1.63,2.55,2.707,2.445,2.325,2.422,1.906,1.86,1.726,1.397,1.779,1.741,1.0525,1.524)

plot(day_in_year,log_number)
lines(lowess(day_in_year,log_number,f=.2,iter=5,delta=1),col="red",lwd=2)

m<-lowess(day_in_year,log_number)
lines(m,col="green",lwd=2)

model<-loess(log_number ~ day_in_year)
xv<-seq(0,366,0.1) # create a sequence 0:366 in increments of 0.1
yv<-predict(model,xv) # use the loess model to predict y values for the xv date values
lines(xv,yv,col="blue",lwd=2)

##########
Page 240

install.packages("maptools")
library(maptools)
data(wrld_simpl)
plot(wrld_simpl)
countries<-wrld_simpl$NAME # contains a list of all the country names in the
file wrld_simpl

countries[18]

##########
Page 241

Australia<-400/3700 # proportion of Peru’s total in Australia
France<-260/3700
Japan<-300/3700
Nigeria<-1319/3700
Peru<-3700/3700
Thailand<-1287/3700
USA<-750/3700

nc<-length(countries)
col.map<-numeric(nc) # this just creates a vector of as many zeros as there are countries (nc), i.e. it’s the same as "col.map<-c(rep(0,nc))"

col.map[which(countries=="Australia")]<-rgb(Australia,1-Australia,0)
col.map[which(countries=="France")]<-rgb(France,1-France,0)
col.map[which(countries=="Japan")]<-rgb(Japan,1-Japan,0)
col.map[which(countries=="Nigeria")]<-rgb(Nigeria,1-Nigeria,0)
col.map[which(countries=="Peru")]<-rgb(Peru,1-Peru,0)
col.map[which(countries=="Thailand")]<-rgb(Thailand,1-Thailand,0)
col.map[which(countries=="United States")]<-rgb(USA,(1-USA),0)
plot(wrld_simpl,col=col.map)

##########
Page 242

install.packages("rworldmap")
library(rworldmap)

data(countryExData) # using the data function in base R to make data available by name
country2Region(regionType="Stern") # to report which countries make up regions, not shown
sternEnvHealth<-country2Region(inFile=countryExData, nameDataColumn="ENVHEALTH",joinCode="ISO3", nameJoinColumn="ISO3V10",regionType="Stern",FUN="mean")


sPDF<-joinCountryData2Map(countryExData,joinCode="ISO3",
nameJoinColumn= "ISO3V10") # ISOs are international standardized country and region codes, see https://en.wikipedia.org/wiki/ISO_3166-1

##########
Page 243

mapCountryData(sPDF,nameColumnToPlot="BIODIVERSITY",mapTitle="World Biodiversity",oceanCol="lightblue",missingCountryCol="white",borderCol="black")

mapCountryData(sPDF,nameColumnToPlot="BIODIVERSITY",mapTitle="Asia",oceanCol="lightblue",missingCountryCol="white",mapRegion="Asia",borderCol="black")

##########
Page 244

mapCountryData(sPDF,nameColumnToPlot="BIODIVERSITY",mapTitle="Australia", oceanCol="lightblue",missingCountryCol="white",xlim=c(110,160),ylim=c(-48,-5), borderCol="black") # for southern hemisphere ylim goes from south to north

mapBubbles(dF=getMap())

##########
Page 245

mapBubbles(dF=getMap(),oceanCol="lightblue",landCol="wheat")

head(getMap()$NAME)

##########
Page 246

mammals<-read.table(text="
Country 'Mammal species' 'threatened spp' Lat Lon
Canada 426 40 60 -95
Brazil 648 80 -34 -64
Mexico 923 96 24 -102
Peru 647 63 -9 -75
Colombia 442 58 3.9 -73
'United States' 440 40 39 -99
Argentina 374 38 -35 -65
Ecuador 372 47 -1 -78
Bolivia 362 21 -17 -64
Chile 143 19 -37 -72
",header=TRUE) # read.table treats sets of spaces like a tab, and does not need textConnection
mapBubbles(dF=mammals,nameZSize="Mammal.species",nameZColour="Country", oceanCol="lightblue",landCol="wheat",nameX="Lon",nameY="Lat",main="Number of Mammal species",symbolSize=0.5)

##########
Page 247

map<-read.csv("Thailand_Border.csv",header=TRUE)

head(map)

##########
Page 248

map[,1]<-(map[,1]-min(map[,1])*(1/(max(map[,1])-min(map[,1]))))
map[,2]<-(map[,2]-min(map[,2])*(1/(max(map[,2])-min(map[,2]))))

plot(NULL,xlim=c(0,1),ylim=c(0,1),axes=F,xlab="",ylab="",main="")
polygon(map)

##########
Page 249

map[,1]<-map[,1]*8.34+97.33 # remember multiplications are always calculated before additions
map[,2]<-map[,2]*14.838+5.612

plot(NULL,xlim=c(94,108.83),ylim=c(5.5,20.5),axes=F,xlab="",ylab="",main="Thailand") # to keep proportions correct we use the same overall range for longitude
polygon(map,col="palegreen1")

##########
Page 250

placidus<-readLines("A_placidus.txt")
head(placidus)

words<-strsplit(placidus," ")
words

##########
page 251

words<-unlist(words)
head(words, 10)

degrees.All<-grep("°",words)

words[degrees.All]

dec.mins<-words[degrees.All+1]
dec.mins

placidus.coords<-cbind(words[degrees.All],dec.mins)
placidus.coords

##########
page 252

placidus.coords[,2]<-gsub("'","",placidus.coords[,2])

placidus.coords[,1]<-gsub("°","",placidus.coords[,1])

placidus.coords[,2]<-as.numeric(placidus.coords[,2])*(100/60) # 60 minutes in a degree, so now they are all decimals of a degree

placidus.coords[,2]<-gsub("\\.","",placidus.coords[,2])

new.degrees<-paste(placidus.coords[,1],".",placidus.coords[,2],sep="")

lat<-(c(1:10)*2)-1 # there are 10 rows in the table
long<-(c(1:10)*2)
points(as.numeric(new.degrees[long]),as.numeric(new.degrees[lat]),pch=10)
mtext ("Aleiodes placidus",font=3) # this is an easy way to add a subtitle to a graph

##########
Page 253

install.packages("dismo")
library(dismo)
bhut<-gbif("Bhutanitis","lidderdalii")

colnames(bhut) # there are a lot of columns, some with long names, which is why the output is ugly

##########
Page 255

localities<-subset(bhut,select=c("lat","lon"))
head(localities)

b<-complete.cases(subset(bhut,select=c("lat","lon")))
localities[b,]

##########
Page 256

newmap<-getMap(resolution="low")
plot(newmap)
points(localities$lon,localities$lat,col="red",cex=2) # use pch to change point type

##########
Page 257-258

tr1<-"(((((BBTH004_16_CCDB_24026_A4_91Tachinidae[90Drosophila]:1.14675336480776657311,(BBTH006_16_CCDB_24026_A6_91Belvosia:0.59297714057027572920,BBTH055_16_CCDB_24026_E7_92Tachinidae[91Belvosia]:0.52736765819178677006):
0.32998852276617751667):0.31714205003432383023,((BBTH084_16_CCDB_24026_G12_97Tachinidae:0.00006954085997979414,BBTH085_16_CCDB_24026_H1_97Tachinidae:0.00006954085997979414):0.00006954085997979414,(BBTH080_16_CCDB_24026_G8_97
Tachinidae:0.00006954085997979414,BBTH091_16_CCDB_24026_H7_97Tachinidae:0.00006954085997979414)"

tr2<-gsub("16_CCDB_24026_","",tr1)
tr2

tr3<-gsub("BBTH.+?24026","BBTH24026",tr1)
tr3

##########
Page 259

tr4<-gsub("\\[.+?\\]","",tr3) # this means remove everything between the
square brackets, [ and ]
tr4

tr5<-strsplit(tr4,"")

tr5

##########
Page 260

tr5<-unlist(tr5)
head(tr5,40) # the 40 tells R how many items in the list tr5 to show

x<-"attgacccggaggaggaatta"
first_G<-regexpr("g",x)
first_G

attr(,"match.length")

attr(,"index.type")

attr(,"useBytes")

all_GGA<-gregexpr("gga",x)
unlist(all_GGA)

all_GGA[[1]][1:length(all_GGA[[1]])]

##########
Page 261

colon<-which(tr5==":")
colon

close<-which(tr5==")")
comma<-which(tr5==",")
close

comma

tr5[colon[1]:(comma[1]-1)] # we have to use () to take the 1 only off the value of comma[1]

##########
Page 262

tr4<-gsub("\\[.+?\\]","",tr3)
repeat{
ss<-unlist(strsplit(tr4,""))
colon<-which(ss==":") # [1] 31 79 128 152 176 228 277 301 351 400
L<-length(unlist(colon)) # 10
if(L==0) break # condition to break out of repeat loop
close<-which(ss==")") # [1] 151 175 300 423
comma<-which(ss==",") # [1] 54 102 199 251 324 374
a<-which(comma>colon[1]) # [1] 1 2 3 4 5 6
b<-which(close>colon[1]) # [1] 1 2 3 4
c<-min(comma[a[1]],close[b[1]]) # 54
n<-substr(tr4,colon[1],c-1) # [1] ":1.14675336480776657311"
tr4<-gsub(n,"",tr4) # removing substring n from tr4
} # end repeat loop
tr4

tr4<-gsub("\\[.+?\\]","",tr3)
tr5<-unlist(strsplit(tr4,""))
colons<-grep(":",tr5)
ends<-colons+22 # 22 including the colon
for(i in length(colons):1) tr4<-gsub(substr(tr4,colons[10],ends[10]),"",tr4)
tr4

##########
Page 262-263

tr4<-gsub("\\[.+?\\]","",tr3)
tr5<-unlist(strsplit(tr4,""))
colons<-grep(":",tr5)
ends<-colons+22 # 22 including the colon
positions<-NULL
for(i in 1:length(colons)){
temp<-colons[i]:ends[i]
positions<-c(positions,temp)}
tr5<-tr5[-positions]

tr4<-paste(tr5,collapse="")
tr4

names<-gsub("\\(","",tr4) # removing (
names<-unlist(strsplit(names,")")) # splitting string at )
names<-unlist(strsplit(names,",")) # splitting string at commas
names

names<-names[which(names!="")]
names

##########
Page 263-264

rearrange<-NULL
for(i in 1:length(names)){
temp<-unlist(strsplit(names[i],"_"))
# temp = "BBTH24026" "A4" "91Tachinidae"
rearrange<-c(rearrange,paste(temp[2],temp[1],temp[3],sep="_"))
} # end i loop
rearrange

sort(rearrange)

for(i in 1:length(rearrange)){
if(substr(rearrange[i],3,3)=="_") {
m<-substr(rearrange[i],1,1)
n<-substr(rearrange[i],2,3)
new<-paste(m,"0",n,sep="")
rearrange[i]<-sub(substr(rearrange[i],1,3),new,rearrange[i])} # end if
} # end i loop
rearrange

sort(rearrange)

##########
Page 265

data<-readLines(textConnection("Aleiodes_bicolor_MRS1008 TATTTTATATTTTTTATTT
Aleiodes_praetor_UK_MRS67_ XXXXXXXXXGTTTTATAT
Aleiodes_rugulosus_CollHH1599_Norway xxxxxxxxxATTTTGTATTTTTT
Aleiodes_seriatus_MRS252_France xxxxxxxxxxxxxxxxxxxxxxxxx
Aleiodes_seriatus_MRS254_France xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
Aleiodes_seriatus_MRS263_France xxxxxxxxxATTTTATACTTTTTATTTGG
Aleiodes_seriatus_MRS264_France xxxxxxxxxATTTTATACTTTTTATTTGG
Aleiodes_seriatus_MRS136_France GATATTGGAATTTTATATTT
MRS239_Aleiodes_seriatus_Russia xxxxxxxxxGTTTTATACTTCTTATTT
Aleiodes_seriatus_MRS222_Germany xxxxxxxxxATTTtaTaCTTTTTATT
Aleiodes_sibiricus_MRS313_Sweden xxxxxxxxxxTTTTGTATTTTTTATT
Aleiodes_signatus_MRS378_Sweden xxxxxxxxxxxTTTATATTTTTTATT
Aleiodes_signatus_MRS712_Sweden GATATTGGTATTTTATATTTTTTA
Aleiodes_unipunctator_CollHH1603_Norway xxxxxxxATTTTATATTTTTTATG"))

grep("MRS",data)

regexpr("MRS",data)

##########
Page 266

unlist(gregexpr(pattern ="MRS[0-9]+",data))

unlist(regmatches(x=data,gregexpr("MRS[0-9]+",text=data)))

startwithMRS<-function(b){
mrs<-unlist(regmatches(x=b,gregexpr("MRS[0-9]+",text=b)))
z<-gsub(mrs,"",b)
z<-paste(mrs,"_",z,sep="")
return(z)} # end startwithMRS function
MRSatstart<-unlist(lapply(data, function(x) startwithMRS(x)))

mrs<-unlist(regmatches(x=data[3],gregexpr("MRS[0-9]+",text=data[3])))
mrs

mrs==character(0)

##########
Page 267

startwithMRS<-function(b){
mrs<-unlist(regmatches(x=b,gregexpr("MRS[0-9]+",text=b)))
if(length(mrs)>0) { # the MRS string occurs on the line
z<-gsub(mrs,"",b) # remove the original MRS number occurrence
z<-paste(mrs,"_",z,sep="") # paste the MRS number at the beginning
} # end if
else z<-b # the MRS string does not occur on the line
z # because the function has produced just the one value, 'z', we can omit 'return(z)'
} # end startwithMRS function
MRSatstart<-unlist(lapply(data,function(x) startwithMRS(x)))
MRSatstart


##########
Page 269

text1<-"We want to write some R code recognizes animal or plant family, subfamily, tribe or subtribe names in a piece of text. For animals these always end in -idae, -inae, -ini and -ina respectively. The Noctuidae, Erebidae, Sphingidae and Saturniidae are families of moths, the Lymantriinae is a subfamily of Erebidae, the Lymantriini, Leucomini, Orgyiini and Nygmiini are tribes of the Lymantriinae, etc."

text2<-unlist(strsplit(text1," "))

lymantrids<-grep("^Lyman",text2)
lymantrids

families<-text2[grep("idae$",text2)]
families

temp<-gsub("[[:punct:]]$","",as.character(text2))
temp

##########
Page 270

families<-temp[grep("idae$",temp)]
families

punctuation<-grep("[[:punct:]]$",as.character(text2))
for(i in 1:length(punctuation)){
word<-text2[punctuation [i]]
word<-paste(substr(word,1,nchar(word)-1),substr(word,nchar(word),
nchar(word)),sep=" ")
text2[punctuation [i]]<-word} # end i loop
text2<-unlist(strsplit(text2," "))
text2

##########
Page 271

newpunct<-grep("^[[:punct:]]$",as.character(text2)) # the ^ and $ make sure it is both the first and last character
newtext<-text2[1] # we assume the first "word" is not a punctuation mark
for(i in 2: length(text2)){
ifelse(i %in% newpunct,newtext<-paste(newtext,text2[i],sep=""),newtext
<-paste(newtext," ",text2[i],sep=""))} # end if
newtext

Fly<-"Tachinid"
Fly<-toupper(Fly)
Fly

Fly<-tolower(Fly)
Fly

casefold(Fly,upper=TRUE)

grepl("[[:upper:]]","Fly")

grepl("[[:upper:]]",tolower(Fly))

##########
Page 272

parasitoids<-c("tachinid","braconid","ichneumonid","eulophid")
parasitoids<-gsub("(\\w)(\\w*)","\\U\\1\\L\\2",parasitoids,perl=TRUE)
parasitoids

parasitoids<-c("tachinid","braconid","ichneumonid","eulophid")
for(i in 1:length(parasitoids)){
parasitoids[i]<-gsub(substr(parasitoids[i],1,1),toupper(substr(parasitoids
[i],1,1)), parasitoids[i])} # end i loop
parasitoids

parasitoids<-c("tachinid","braconidae","ichneumonid","eulophidae")

family_names<-grep(".ae$",parasitoids)
family_names

substr(parasitoids[family_names],1,1)<toupper(substr(parasitoids
[family_names],1,1))
parasitoids

family_names<-grep(".ae$",parasitoids,value=TRUE)
family_names

##########
Page 273

punct<-'[]\\?!\"\'#$%&(){}+*/:;,._`|~\\[<=>@\\^-]'

punct2<-sub("%","",punct) # removing the % symbol
punct2

tx<-"Numbers, and percentages, of the fish were: 9 (14%), 12 (19%), and 2
(2%), respectively"
gsub(punct2,"",tx)

##########
Page 275

install.packages("ape")
library(ape)

##########
Page 276

packageVersion("ape")

insecta<-"(dragonflies,(cockroaches,(bugs,(wasps&bees,((beetles,lacewings),
(flies, butterflies&moths))))));" # note the semicolon at the end of the line and that the taxon names do not need to be in inverted commas
insectTree<-read.tree(text=insecta) # or we can read from a file, read.tree
({file name})

summary(insectTree)

##########
Page 277

str(insectTree)

attributes(insectTree)

insectTree$edge

install.packages("phytools")
library(phytools)
plotTree(insectTree,offset=1) # this phytools function allows us to add some other information
tiplabels() # rather like adding points or lines to an open plot in base R
nodelabels()

##########
Page 278

plot(insectTree)

##########
page 279

??plot

insectTree$edge.length
<-c(370,20,350,40,310,50,260,80,50,130,130,80,100,100)
write.tree(insectTree,"Insect_tree.txt")
readLines("Insect.tree")

##########
Page 280

random_tree25<-rtree(n=25)
random_tree25 # you probably won't get the same answer because the tree is random

plot(random_tree25)

random_tree25 [2]

random_tree25$tip.label

##########
Page 281

par(mfrow=c(2,2))
plot(random_tree25,edge.width=2,label.offset=0.1,type ="cladogram")
plot(random_tree25,edge.width=2,label.offset=0.1,type ="phylogram")
plot(random_tree25,edge.width=2,label.offset=0.1,type ="fan")
plot(random_tree25,edge.width=2,label.offset=0.1,type ="radial")

random_tree25$tip.label[3]

random_tree25$tip.label[3]<-"apple"

Nodonata<-drop.tip(insectTree,"dragonflies")

wrong<-root(Nodonata,7,resolve.root=TRUE)
plot(wrong)

##########
Page 284

p<-unlist(strsplit("ATCTCGGCGCGCATCGCGTACGCTACTAGC",""))
p

map=c("A"="T","T"="A","G"="C","C"="G")
map[p]

##########
Page 285

unname(map[p])

sapply(p,switch, "A"="T","T"="A","G"="C","C"="G")

RCOMP<-function(x){
r<-sapply(lapply(strsplit(x,NULL),rev),paste,collapse="") # this line splits the sequence to individual bases, reverses the order and recombines
y<-paste(map[unlist(strsplit(r,NULL))],collapse="") # this line replaces the bases with their complements
return(y)} # end RCOMP function

install.packages("seqinr")
library(seqinr)

comp(p)

rev(comp(p)) # gives the reverse complement, i.e. the complement read in
the reverse direction

##########
Page 286

utf8ToInt("abcdefghijklmnopqrstuvwxyz")

intToUtf8(utf8ToInt("abcdefghijklmnopqrstuvwxyz"))

c2s(comp(s2c("AAAATTTTGGGGCCCC"))) # comp accepts upper or lower case but returns lower case

c2s(rev(comp(s2c("aaaattttggggcccc"))))

allbases<-s2c("atgcmynnsk")
comp(allbases) # NAs are produced for non-ATGC bases

comp(allbases,ambiguous=TRUE) # no more NAs

install.packages("BiocManager")
library(Biostrings)
dna=DNAStringSet(c("ATCTCGGCGCGCATCGCGTACGCTACTAGC","ACCGCTA"))
complement(dna)

##########
Page 288

rnasq<-"UAUAUUUCAAAUAUAAAGUAAUUUAGAAAGUAAGAUCAAAGAUUAAU"
types<-"PPPPPUUUUPPPPPPPUPPUUUPPPPUUUUUUUUUUPPPPPPUUUUU"
U<-which(unlist(strsplit(types,""))=="U")
P<-which(unlist(strsplit(types,""))=="P")

Pgroups<-split(P,cumsum(c(1,diff(P) != 1)))
Pgroups

##########
Page 289

diff(P)

c(1,diff(P))

c(1,diff(P) !=1)

cumsum(c(1,diff(P) != 1))

randomDNA<-paste(sample(c("A","C","G","T"),10000,rep=TRUE,prob=c(0.4,0.1,0.1,0.4)), collapse="")

##########
Page 290

library(ape) # use install.packages("ape") if you have not done so in Chapter 23
Monilobracon1<-read.GenBank("AY529648",as.character=TRUE)
Monilobracon1

mode(Monilobracon1)

##########
Page 291

str(Monilobracon1)

Monilobracon1$AY529648[1:80]

s<-paste(Monilobracon1$AY529648,collapse="")
s

attributes(Monilobracon1)[1]

attributes(Monilobracon1)[2]

all_Monilos<-c("MH234998","MH260676","AY296646","AY529647","AY532320", "AY529649","AJ296046")
Mons<-read.GenBank(all_Monilos,as.character=TRUE)

data<-cbind(attr(Mons,"species"),names(Mons))
data

##########
Page 292

write.dna(Mons,file="Monilobracon.fas",format="fasta",nbcol=6,colsep=" ",colw=10)

SGC4<-getGeneticCode("SGC4")

barcode<-paste(Mons$MH260676,collapse="")
barcode

dna1<-DNAString(barcode) # we have to use the function DNAString to tell Biostrings what type of object our sequence is

dna2<-DNAString(substr(barcode,2,nchar(barcode)))
dna3<-DNAString(substr(barcode,3,nchar(barcode)))

translate(dna3,genetic.code=SGC4,no.init.codon=TRUE)

##########
Page 294

43500000000000001/1500000000000000

all.equal(8.00000012,round(8.00000012))

all.equal(8.00000011,round(8.00000011))

is.whole<-function(a) a%%1==0 # %% means modulo, i.e. remainder after division

is.whole<-function(a) round(a)==a
is.whole(7.0000000000000004)

is.whole(7.0000000000000005)

Xenarcha<-"000000010000001011101111000000000?01?100001000?0?0100"
Colastes<-"000000010000011001101011000000000?01?110101000?0?0100"
Gnamptodon<-"000000010000001001001001001000111?01?0?0001?00?0?0101"
wasps<-rbind(Xenarcha,Colastes,Gnamptodon)
wasps

##########
Page 295

newwasps<-wasps
for(i in 1:length(wasps)){
newdata<-NULL
for(j in 1:nchar(wasps[i])){
ifelse(is.whole(j/5),newdata<-paste(newdata,substr(wasps[i],j,j),
"",collapse="",sep=""),newdata<-paste(newdata,substr(wasps[i],
j,j), collapse="",sep="")) # end ifelse
} # end for j loop
if(i==1) ONE<-TRUE
if(i==2) TWO<-TRUE
newwasps[i]<-newdata} # end i loop
newwasps

tuatara<-readLines("Tuatara.txt")
tuatara

##########
Page 296

names<-NULL
for(i in 1:length(tuatara))
names<-c(names,unlist(strsplit(tuatara[i]," "))[1])
names

names<-unlist(lapply(strsplit(tuatara,split=" "),"[",1))
names

tuatara[6]<-gsub("*\\(.*?\\)","?", tuatara[6])

##########
Page 297

data=read.csv("ActiniaXYdata.csv",header=TRUE)

##########
Page 298

data$X<-round(data$X,1);data$Y<-round(data$Y,1) # multiple R commands
on one line separated by a semicolon – or we can just use data<-
round(data,1)
pdd<-c(4.9,1.3,5.3,1.6,3.9,3.8,1.8,3.5,3,2.9,3,3.9,3,3.3,2.1,4.5,3.2,4.2,6.1,3.8,4,
4,3.8,1.5, 3.9,4.3,4.1,3.9,1.8,3.2,4.2,3.2,1.9,4,2,3.8,1.5,5.9,3.9,4.1,3,3.6,
3.2,3.8,3.9,2.5,5.5,3,1.5,3.9,2.5,4,1.2,4,3.9,3.7,1.8,4.4,3.2,3.9,3.7,3.7,
4.4,3.9,4.6,4.2,5,4,4.6,5,3.7,3,4.1,1.5,4,3.2,3,1.7,4,3,3.1,3.8,4.2,5,2,4,
4.1,3.9,3.1,3,1.6,3.7,3.5,3,3.8,3.9,4,4.9)
plot(data,pch=19,xlim=c(0,50),ylim=c(0,50),col="darkred",cex=pdd/1.6,
xlab="X coordinate (cm)",ylab="Y coordinate (cm)") # cex is adjusted to
scale diameter of anemones correctly

rect(0,0,50,50,lty=4,lwd=0.6)

##########
Page 299

hypotenuse<-function(x1,y1,x2,y2) sqrt((x1-x2)^2+(y1-y2)^2)

hypotenuse(0,0,3,4)

N<-nrow(data)
distances<-matrix(nrow=N,ncol=N) # the cells of the matrix are automatically
filled with NAs if nothing is specified
for(i in 1:N){
for(j in 1:N){
distances[i,j]<-hypotenuse(data$X[i],data$Y[i],data$X[j],data$Y[j])
}}

NNDs<-NULL
for(i in 1:N) NNDs<-c(NNDs,min(distances[i,which(distances[i,]>0)]))

hist(NNDs,col="pink")

##########
Page 300

RE<-1/(2*sqrt(0.0392))
RE

R<-mean(NNDs)/RE
R # index of dispersion

##########
Page 301

t.test(NNDs,mu=RE)

SR<-0.26136/sqrt(N* 0.0392)
SR

z<-(mean(NNDs)-RE)/SR
z

plot(sort(NNDs)^2,-log(1-((1:N)/(N+1))),ylab="Weiss Index",xlab="NND^2",pch=18,main="Weiss plot of anemone spacing")
lines(sort(NNDs)^2,-log(1-((1:N)/(N+1))),col="red")

##########
Page 304

ricker<-function(Nt,r){ # for simplicity K is taken to be 1
NtPlus1<-Nt*exp(r*(1-Nt))
return(NtPlus1) # you can omit "return" and just write NtPlus1
}
minimum_r<-0.001 # set the starting value of r
maximum_r<-5
r_sequence<-seq(minimum_r,maximum_r,0.001) # creates a vector of r-values that will be stepped through:
plot(NULL,xlim=range(r_sequence),ylim=c(0,5),xlab="r",ylab="Nt",
main="Ricker logistic population model") # create a blank plotting field
for(r in r_sequence){ # note that here we have presented a sequence of real numbers so we cannot use a colon to separate extremes
Nt<-runif(1) # picks a random value for Nt between 0 and 1 to start with
for(i in 1:200) Nt<-ricker(Nt,r) # runs the simulation for 200 generations ignoring the result, you can think of this as a burn-in phase
for(i in 1:200){Nt<-ricker(Nt,r) # now calculates and plots the population size for each of the 200 generations after the burn in
points(r,Nt,pch=".",cex=1.5)} # plotting a full stop because there are going to be so many points
} # end of r loop

##########
Page 305

ricker<-function(Nt,r) Nt*exp(r*(1-Nt))

##########
Page 306

i<-1
while(i < 500000000){i<-i+1}; beep(8)

cat("press [escape key] to stop"); i<-1; while(i<2) beep(9)

system("say Program finished!")

Nt<-1000
Pt<-20
alpha<-0.002
r<-3
c<-1
K<-2000

##########
Page 307

hosts<-Nt # starting these two vectors with the initial host and parasitoid populations
parasitoids<-Pt
number_of_generations<-1000
for(i in 1: number_of_generations){
# THE MAIN BODY TEXT WILL GO HERE
}

Beddington<-function(fNt,fPt,fr,falpha,fK,fc){
NtPlus1<-fNt*exp((fr*(1-(fNt/fK)))-(falpha*fPt))
PtPlus1<-fc*fNt*(1-(exp(-falpha*fPt)))
answer<-c(NtPlus1,PtPlus1)
answer
}

Nt<-1000
Pt<-20
alpha<-0.002
r<-3
c<-1
K<-2000 # we put the user-defined constants first to make editing them easier
Beddington<-function(fNt,fPt,fr,falpha,fK,fc){
NtPlus1<-fNt*exp((fr*(1-(fNt/fK)))-(falpha*fPt))
PtPlus1<-fc*fNt*(1-(exp(-falpha*fPt)))
answer<-c(NtPlus1,PtPlus1)
answer # this passes the value of "answer" to where the function was called
} # end Beddington function

##########
Page 308

hosts<-Nt # starting these two vectors with the initial host and parasitoid populations
parasitoids<-Pt
number_of_generations<-1000
for(i in 1: number_of_generations){
answer<-Beddington(Nt,Pt,r,alpha,K,c)
hosts<-c(hosts,answer[1])
parasitoids<-c(parasitoids,answer[2])
Nt<-answer[1]
Pt<-answer[2]} # end of i loop
plot(1:length(hosts),hosts,ylim=c(0,max(hosts)),col="green",pch=20,xlab="Generations",ylab="Population size",main="Parasitoid-Host
Population Dynamics")
points(1:length(hosts),parasitoids,col="red",pch=20)
# there is no room for a legend so we have to insert the legend manually with a bit of trial and error for positions
points (number_of_generations/4.5,max(hosts)*1.07,col="green",pch=20,xpd=NA)
text(number_of_generations/3.5,max(hosts)*1.07,"Hosts",xpd=NA)
points(number_of_generations/2,max(hosts)*1.07,col="red",pch=20,xpd=NA)
text (number_of_generations/1.65,max(hosts)*1.07,"Parasitoids",xpd=NA)

##########
Page 309

plot(hosts,parasitoids,pch=20,xlab="Host population size",ylab="Parasitoid population size",main="Parasitoid-Host Phase Plot",cex=0.5)

plot(hosts,parasitoids,pch=20,xlab="Host population size",ylab="Parasitoid population size",main="Parasitoid-Host Phase Plot",cex=0.5)
lines(hosts,parasitoids)

##########
Page 311

Z<-30
H<-matrix (NA,nrow=Z,ncol=Z) # number of hosts
P<-matrix (NA,nrow=Z,ncol=Z) # number of parasitoids

G<-100

Beddington<-function(fNt,fPt,fr,falpha,fK,fc){
NtPlus1<-fNt*exp((fr*(1-(fNt/fK)))-(falpha*fPt))
PtPlus1<-fc*fNt*(1-(exp(-falpha*fPt)))
answer<-c(NtPlus1,PtPlus1)
answer}

##########
Page 312

random_asign_to_neighbours<-function(to_Add_f,Number_moving,neighbours){
breaks<-sort(sample(1:Number_moving,7,replace=TRUE))
cell<-rep(NA,8) # because we are going to assign values to numbered positions in the vector "cell" we need to define it first
for(i in 2:7){
cell[i]<-breaks[i]-breaks[i-1]} # end i loop
cell[1]<-breaks[1]
cell[8]<-Number_moving-breaks[7]
for(x in 1:8){
to_Add_f[neighbours[x,1],neighbours[x,2]]<-to_Add_f
[neighbours[x,1],neighbours [x,2]]+cell[x]}
to_Add_f # same as return(to_Add_f), the last variable will be returned
} # end of random_asign_to_neighbours function

plot_generation<-function(X,Z){ # X receives the host or parasitoid matrices
mcol<-max(X) # to colour each plotted symbol in no particular way
for(j in 1:Z){
for(k in 1:Z){
points(j,k,col=round(10*X[j,k]/mcol),pch=15)
}} # end j and k loops
} # end plot_generation function

H[1:Z,1:Z]<-100
P[1:Z,1:Z]<-10 # we could use P<-matrix(10,nrow=Z,ncol=Z)

alpha<-0.0025
r<-2.1 # intrinsic rate of host population increase
c<-1
K<-1000 # carrying capacity

##########
Page 313

proportion_to_migrate<-0.1 # 10% of individuals in each cell will migrate to
neighbouring cells
neighbours<-matrix(NA,nrow=8,ncol=2) # hold the relative coordinates of
neighbouring cells

save_runs<-c(10,30,60,100) # these will be used to make snapshots after these numbers of generations; we plot 10 and 100 at the end
hold_H<-matrix(NA,nrow=Z*Z,ncol=length(save_runs)) # these will store the four host population snapshots
hold_P<-matrix(NA,nrow=Z*Z,ncol=length(save_runs)) # ditto for parasitoids count<-0 # we will use this to set the column in the hold matrices to fill
par(mfrow=c(3,2)) # set graphics parameter for plotting spatial arrays during simulation

for(g in 1:G) { # stepping through generations
for(ba in 1:Z){ # stepping through rows
for(bb in 1:Z){ # stepping through columns
answer<-Beddington(H[ba,bb],P[ba,bb],r,alpha,K,c)
H[ba,bb]<-answer[1]
P[ba,bb]<-answer[2]
} # end bb loop
} # end ba loop

tomove_Host<-round(H* proportion_to_migrate)
tomove_Parasitoid<-round(P* proportion_to_migrate)
H<-round(H-tomove_Host) # subtracts the values of cells in tomove_Host from corresponding cells in H, leaving just the residents
P<-round(P-tomove_Parasitoid) # ditto for parasitoids

##########
Page 314

to_add_Host<-matrix(0,nrow=Z,ncol=Z)
to_add_Parasitoid<-matrix(0,nrow=Z,ncol=Z)
for(ba in 1:Z){ # stepping through rows
for(bb in 1:Z){ # stepping through columns
neighbours[1,1:2]<-c(ba-1,bb-1)
neighbours[2,1:2]<-c(ba-1,bb)
neighbours[3,1:2]<-c(ba-1,bb+1)
neighbours[4,1:2]<-c(ba,bb-1)
neighbours[5,1:2]<-c(ba,bb+1)
neighbours[6,1:2]<-c(ba+1,bb-1)
neighbours[7,1:2]<-c(ba+1,bb)
neighbours[8,1:2]<-c(ba+1,bb+1)
neighbours[which(neighbours==0)]<--999 # minus 999 is an impossible holding value
neighbours[which(neighbours==Z+1)]<-1
neighbours[which(neighbours==-999)]<-Z
to_add_Host<-random_asign_to_neighbours(to_add_
Host,H[ba,bb], neighbours)
to_add_Parasitoid<-random_asign_to_neighbours(to_add_
Parasitoid, P[ba,bb],neighbours)
} # end bb loop
} # end ba loop

H<-H + to_add_Host
P<-P + to_add_Parasitoid
plot(NULL,xlim=c(1,Z),ylim=c(1,Z),axes=F,xlab="",ylab="",main=paste
("generation =",g)) # the default sep argument in paste is " "
mtext("Host")
plot_generation (H,Z)
plot(NULL,xlim=c(1,Z),ylim=c(1,Z),axes=F,xlab="",ylab="",main=paste
("generation =",g))
plot_generation (P,Z)
mtext("Parasitoid")

##########
Page 315

if(g %in% save_runs){
count<-count+1
hold_H[,count]<-unlist(expand.grid(H)) # you must use unlist because expand.grid returns a list of length 1
hold_P[,count]<-unlist(expand.grid(P))
} # end if g
} # end G loop

getcoloursrainbow<-function(x){
purplef<-220
f<-unlist(as.list(rainbow(256))) # rainbow is a built-in colour palette comprising 256 shades
mn<-0
mx<-max(x,na.rm=T)
cor<-(x-mn)/(mx-mn)
ans<-f[(1+ceiling(purplef*(cor)))] # ceiling returns the integer immediately above a real number, cf. floor and round
ans} # end function getcoloursrainbow

##########
Page 315-316

par(mfrow=c(2,2))
point_colours<-getcoloursrainbow(hold_H[,1])
point_colours2<-matrix(point_colours,nrow=Z)
point_colours3<-getcoloursrainbow(hold_P[,1])
point_colours4<-matrix(point_colours3,nrow=Z)
plot(NULL,xlim=c(1,Z),ylim=c(1,Z),axes=F,xlab="",ylab="",main="Host")
for(j in 1:Z){
for(k in 1:Z){
points(j,k,col=point_colours2[j,k],pch=15)
}} # end j and k loops
plot(NULL,xlim=c(1,Z),ylim=c(1,Z),axes=F,xlab="",ylab="",main="Parasitoid")
for(j in 1:Z){
for(k in 1:Z){
points(j,k,col=point_colours4[j,k],pch=15)
}} # end j and k loops
point_colours<-getcoloursrainbow(hold_H[,4])
point_colours2<-matrix(point_colours,nrow=Z)
point_colours3<-getcoloursrainbow(hold_P[,4])
point_colours4<-matrix(point_colours3,nrow=Z)
plot(NULL,xlim=c(1,Z),ylim=c(1,Z),axes=F,xlab="",ylab="",main="Host")
for(j in 1:Z){
for(k in 1:Z){
points(j,k,col=point_colours2[j,k],pch=15)
}} # end j and k loops
plot(NULL,xlim=c(1,Z),ylim=c(1,Z),axes=F,xlab="",ylab="",main="Parasitoid")
for(j in 1:Z){
for(k in 1:Z){
points(j,k,col=point_colours4[j,k],pch=15)
}} # end j and k loops

##########
Page 316

centreX<--7 # minus 7; you do not need to put a space between the two -s
centreY<-45
LX<-8
lowY<-centreY-2
highY<-centreY+2
increment<-0.08
line_colour<-getcoloursrainbow(1:220)
for(i in 1:220){
X<-centreX-LX+(i*increment)
lines(c(X,X),c(lowY,highY),col=line_colour[i],xpd=NA)}
text(xpd=NA,-24,centreY,"Population = 0") # xpd=NA can go anywhere
text(xpd=NA,15,centreY,"Maximum population")

##########
Page 319

par(mfrow=c(3,2)) # to create a two columns of three rows plotting area
par(mar=c(2,2,2,2))
alleles<-c(rep("red",200),rep("blue",200),rep("green",200),rep("yellow",200),
rep("turquoise",200))
alleles<-sample(alleles,1000,replace=F)
plotlist<-c(1,5,10,25,50,100)

plotpoints<-function(all,gen){
plot(NULL,axes=F,xlim=c(0,100),ylim=c(0,100),ylab="",xlab="",main="")
for(j in 1:length(all)){ # stepping through vector of genotypes to plot
x<-sample(1:1000000,1)/10000 # randomly selecting and scaling it to axis length of 100 to minimize the probability of points precisely overlapping
y<-sample(1:1000000,1)/10000
points(x,y,pch=15,col=alleles[j])} # end j loop
text(xpd=NA,50,110,paste("generation = ",gen,sep=""),cex=0.8) # adding the generation number as a label centred above each plot
} # end plotpoints function

for(i in 1:100){
if(i %in% plotlist) plotpoints(alleles,i) # we don't need to use {} as only
one operation to be performed
breeders<-sample(alleles,100,replace=F)
alleles<-rep(breeders,10)} # end i loop

##########
Page 322

a<-rnorm(100000)

##########
Page 323

start_time<-proc.time() # start computer clock
for(i in 1:10000000) a[i]<-a[i] + 1 # loop through the vector
proc.time()-start_time # the computer clock returns a new value to proc.time

a<-rnorm(10000000)
start_time<-proc.time() # record computer clock time
a<-a+1
proc.time()-start_time

a<-c(1,2,3,4)
b<-c(0,6,3,12)
c<-c(0,2,3.1,1000)
d<-rbind(a,b,c)
d

apply(d,1,mean) # the 1 means by row, there are 3 rows so 3 means

apply(d,2,mean) # the 2 means apply the function to columns, there are 4
means

apply(d,2,function(x) length(x[x<10]))

apply(d,1,function(x) x>5) # the 1 in this case means by row so the answers are given as rows × columns

##########
Page 324

apply(d,2,function(x) x>5) # the 2 in this case means by column so the answers are given as columns × rows

sapply(c(1,2,3,4,5),function(x) x^2) # you could replace c(1,2,3,4,5) by 1:5

lapply(1:3,function(x) x^2) # gives the same answers as sapply but in the form of a list

unlist(lapply(1:3,function(x) x^2))

Tilapia<-5
Mystus_singaringan<-c(1,2,2,2)
Osteochilus_vittatus<-c(26,30,11,2,28,6,19,12,4,4,18,13,29,21,9,8,21,19,28)
Henicorhynchus_siamensis<-c(46,22,33,28,50,42,31,4,13,8,11,8,18)

names<-c("Tilapia","Mystus","Osteochilus","Henicorhynchus")
wize<-lengths(list(Tilapia,Mystus_singaringan,Osteochilus_vittatus,
Henicorhynchus_siamensis))
fishnames<-c(rep(names,wize))

##########
Page 325

metacercaria<-c(Tilapia,Mystus_singaringan,Osteochilus_vittatus,
Henicorhynchus_siamensis)
df<-data.frame(fishnames,metacercaria)
head(df,7)

tapply(df$metacercaria,df$fishnames,mean)

tapply(df$metacercaria,df$fishnames,range)

##########
Page 327

parasitoid<-c("Bracon1","Bracon1","Bracon1","Bracon1","Bracon1","agathidine1","agathidine1","agathidine1","agathidine1","Cotesia1","Cotesia1","Cotesia1","Cotesia1","Cotesia2","Cotesia3","Cotesia3","Cotesia3","Cotesia4","Aleiodes1","Aleiodes2","Macrocentrus1")
host<-c("Nymphalidae2","Erebidae1","Erebidae1","Erebidae1","Erebidae1","Erebidae2", "Erebidae2","Noctuid1","Lycaenidae1","Hesperidae1","Hesperidae1","Hesperidae1", "Hesperidae2","Noctuidae2","Geometridae1","Geometridae2","Erebidae1","Erebidae1","Erebidae1","Erebidae3","Erebidae4")

l<-levels(as.factor(parasitoid))
m<-levels(as.factor(host))

plot(NULL,axes=F,xlab="",ylab="",xlim=c(0,length(parasitoid)),ylim=c(0,5))

##########
Page 327-328

P<-NULL
H<-NULL
for(i in 1:length(l)) P<-c(P,length(which(parasitoid==l[i])))
for(i in 1:length(m)) H<-c(H,length(which(host==m[i])))
S<-cumsum(P)
S1<-c(0,S)
T<-cumsum(H)
T1<-c(0,T)
for(i in 1:length(l)) rect(S1[i],3.9,S1[i+1],5.7,col="red",lwd=6,border="white")
for(j in 1:length(m)) rect(T1[j],0.7,T1[j+1],2,col="green",lwd=6,border="white")
midT<-T1[1:length(T)]+(T-T1[1:length(T)])/2
midS<-S1[1:length(S)]+(S-S1[1:length(S)])/2
for(i in 1:length(l)) text(midS[i],4.1,l[i],srt=90,cex=0.8,xpd=NA,adj=0,font=3)
for(i in 1:length(m)) text(midT[i],1.9,m[i],srt=270,cex=0.8,xpd=NA,adj=0,font=3)
for(k in 1:length(l)){ # stepping through each different parasitoid
pp<-which (parasitoid==l[k])
np<-length(pp)
hostsofparasitoid<-host[pp]
links<-length(levels(as.factor(hostsofparasitoid)))
for(q in 1:links){
noRearings<-length(which(hostsofparasitoid==levels(as.factor
(hostsofparasitoid))[q]))
temph<-which(m==levels(as.factor(hostsofparasitoid))[q])
lines(c(midS[k],midT[temph]),c(3.8,2.2),col="grey50",lwd=noRearings*2,lend=3)
} # end q loop
} # end k loop
points(midT,rep(2.2,length(T)),pch=15,col="blue") # putting the squares at the end of the connecting lines
points(midS,rep(3.8,length(S)),pch=15,col="blue")

##########
Page 329

install.packages("cheddar")
library(cheddar)
options(stringsAsFactors=F) # this line is necessary because cheddar does not automatically recognize character strings as factors

properties<-data.frame(title="ThaiWasps")
properties

nodes<-data.frame(node=unique(c(host,parasitoid)),type=c(rep("host",length(unique (host))),rep("parasitoid",length(unique(parasitoid)))),colour=c(rep("green", length(unique(host))),rep("red",length(unique(parasitoid)))))
head(nodes)

all.records<-as.data.frame(cbind(host,parasitoid)) # this is just our data as a dataframe
unique.links<-as.data.frame(unique(all.records)) # unique finds unique rows in a dataframe
head(unique.links)

##########
Page 330

hp_pairs<-paste(all.records$host,all.records$parasitoid)
hp_pairs

weight<-as.vector(table(hp_pairs)) # using as.vector just gives us the numbers
weight

unique.links$weight<-weight
colnames(unique.links)<-c("resource","consumer","weight")

thaiFW<-list(nodes, properties,unique.links)
names(thaiFW)<-c("nodes","properties","trophic.links")
attr(thaiFW,"class")<-c("Community","list")

par(mfrow=c(1,3))

PlotCircularWeb(thaiFW,show.nodes.as="labels",node.labels="node",col=NPS(thaiFW)$colour,main="Thai host-parasitoid food web",xpd=NA)

PlotCircularWeb(thaiFW,show.nodes.as="labels",node.labels="node", col=NPS(thaiFW)$colour,link.lwd=TLPS(thaiFW)$weight,main="Thai host-parasitoid food web",xpd=NA)

##########
Page 331

PlotWebByLevel(thaiFW,pch=21,bg=NPS(thaiFW)$colour,link.lwd=TLPS(thaiFW)$weight, main="Thai host-parasitoid food web",xpd=NA)

##########
Page 332

install.packages("jpeg")
library(jpeg)

install.packages("tiff")
library(tiff)

imgjpg<-readJPEG("Oriental_Pied_Hornbill.jpg",native=TRUE)

##########
Page 333

par("mar")

par(mfrow=c(1,1))
plot(c(0,3*22.56),c(0,4*12.7),type="n",xlab="",ylab="")

for(i in 1:3){
for(j in 1:4){
rasterImage(imgjpg,(i-1)*22.56,(j-1)*12.7,i*22.56,j*12.7)
}} # end i and j loops

##########
Page 335

head(rnorm(100,20,4),20)

hist(rnorm(100,20,4))

hist(rnorm(1000000,20,4))

##########
Page 336

hist(rnorm(1000000,20,4))
xv<-seq(-5,45,0.1) # from –5 to 45 in steps of 0.1
yv<-dnorm(xv,mean=20,sd=4)*2000000
lines(xv,yv,col="red",lwd=2)

##########
Page 337

xv<-seq(-20,20,0.1) # from –20 to 20 in steps of 0.1
yv<-pnorm(xv,mean=0,sd=5)
plot(xv,yv,col="red",lwd=2,type="l",main="Cumulative values of the normal distribution",xlab="X",ylab="Cumulative P") # type= "l" tells R to plot lines
mtext("mean =1; S.D. =5")

##########
Page 337-338

lower95<-pnorm(-1.96*5,mean=0,sd=5)
upper95<-pnorm(1.96*5,mean=0,sd=5)
lines(c(-1.96*5,-1.96*5),c(0,lower95),col="blue",lwd=3,lty=2,lend=2)
lines(c(1.96*5,1.96*5),c(0,upper95),col="blue",lwd=3,lty=2,lend=2)
xrange<-1.96*5 - -1.96*5 # 19.6, same as xrange<-1.96*5--1.96*5
x_increment<-xrange/1000
for(i in 7:993){ # not from 1 to 1000 because we don't want to obliterate the blue dashed lines
x<--1.96*5+(i* x_increment)
lines(c(x,x),c(0,pnorm(x,mean=0,sd=5)),col="grey80",lwd=0.5)
} # end i loop
rect(-1.5,0.152,7.5,0.275,col="white",bty="n")
text(3,0.25,"95% of",cex=0.8)
text(3,0.21,"observations",cex=0.8)
text(3,0.17,"are in grey zone",cex=0.8)
lines(xv,yv,col="red",lwd=2) # redrawing the red line that has been slightly covered by grey

##########
Page 339

sample_size<-10000
degf<-10
x<-rt(sample_size,degf)
hist(x,col="grey30",ylab="Probability and cumulative probability",xlab="t",
prob=TRUE,ylim=c(0,1),main="Histogram of t with 10 D.F.")
curve(dt(x df=degf),lwd=2,col="green",add=TRUE)
curve(pt(x,df=degf),lwd=2,col="blue",add=TRUE)
lines(c(qt(0.975,degf),qt(0.975,degf)),c(0,0.975),col="red",lty=2,lwd=2)
lines(c(-100,qt(0.975,degf)),c(0.975,0.975),col="red",lty=2,lwd=2)
text(2.65,0.5,paste("P=0.05,critical (2-tailed) t = ",round(qt(0.975, degf),4),
sep=""),col="red",srt=270,cex=0.8)
lines(c(qt(0.025,degf),qt(0.025,degf)),c(0,0.025),col="red",lty=2,lwd=2)
lines(c(-100,qt(0.025,degf)),c(0.025,0.025),col="red",lty=2,lwd=2)

##########
Page 340

round(qt(0.95,1:10),4)
round(qt(0.99,1:10),4)
round(qt(0.975,1:10),4)
round(qt(0.995,1:10),4)

##########
Page 341

plot(seq(0,10,0.01),dlnorm(seq(0,10,0.01)),type="l",lwd=2,xlab="X",ylab="
Lognormal (x)")

xv<-seq(-10,10,0.1)
a<-1
b<--0.5 # changing beta to negative makes the function start at 1 and decrease to zero
P<-a/(a+ exp(a+b*xv))
plot(xv,P,ylab="Probability",xlab="X",type="l",lwd=2)

##########
Page 342-343

inc<-0.085
p<-c("red","blue","green","purple","brown","pink","seagreen2","violet","orange","turquoise") # a vector of colours for each mean value
plot(NULL,xlim=c(0,20),ylim=c(0,0.4),xlab="X",ylab="Poisson (x)",main="Poisson distributions with means 1:10")
for(i in 1:10){
for(j in 0:20){
lines(c(j+(i-1)*inc,j+(i-1)*inc),c(0,dpois(j,i)),lwd=3,col=p[i])
}}
par(new=TRUE)
for(i in 1:10) lines((0:20)+(i-1)*inc,dpois(0:20,i),col=p[i])
legendLabels=c(paste("mean =",1:10)) # we want a space after = so collapse not used
legend("topright",inset=0.05,legend=legendLabels,col=p,lty=1,cex=0.8)

##########
Page 344

par(mfrow=c(2,2))
xv<-seq(0.001,5,0.01)
plot(xv,dgamma(xv,0.5,0.1),type="l") # type= "l" tells R to plot lines not points
plot(xv,dgamma(xv,0.5,5),type="l")
plot(xv,dgamma(xv,1.5,0.5),type="l")
plot(xv,dgamma(xv,1.5,5),type="l")

par(mfrow=c(1,2))
sample_size<-10000
x<-rchisq(sample_size,5)
hist(x,col="grey30",ylab="Frequency")

##########
Page 345

hist(x,prob=TRUE,col="grey30",ylab="Probability")

##########
Page 346

par(mfrow=c(1,1))
hist(x,prob=T,col="grey30",ylab="Probability")
curve(dchisq(x,df=2),lwd=2,col="green",add=TRUE)
curve(dchisq(x,df=5),lwd=2,col="blue",add=TRUE)
curve(dchisq(x,df=10),lwd=2,col="red",add=TRUE)
leg_text<-c("d.f. = 2","d.f. = 5","d.f.=10")
legend("topright",lwd=2,col=c("green","blue","red"),leg_text)

##########
Page 348

txt<-readLines("/Users/donald/Desktop/books/testgrab.txt")

newdata<-read.table(file="/Users/donald/Desktop/DougChestersCo1/Data_
Release_1/ ChestersDataSumR.txt",header=T) # this will obviously only work on DQ’s computer

write(file="/Users/donald/Desktop/DougChestersCo1/newdata.txt")

##########
Page 349

write(x,file="newdata.txt") # if you just want the contents of ‘x’ saved – suitable for text

write(1:27,file="test.write.txt")
readLines("test.write.txt")

scan("test.write.txt",quiet=TRUE) # 'quiet' suppresses reporting how many items were read

write(file="newdata.txt",append=TRUE)

##########
Page 350

zamias<-read.delim("Panamanian_cycads_colon.txt",sep=";",header=F)
zamias

unlist(unname(zamias))

##########
Page 351

test<-readLines("PrimersTable.xlsx")

head(test,3) # this shows the top three lines in the file test

test<-read.table("PrimersTable.txt",header=TRUE)
test

##########
Page 352

install.packages("openxlsx")

library(openxlsx)

pelagic1<-read.xlsx("https://doi.org/10.1371/journal.pone.0072683.s003")

pelagic_respiration<-read.xlsx("https://doi.org/10.1371/journal.pone.0072683.s003",sheet=2) # note drop the .XLSX

ordinary dataframes to tibbles.
install.packages("readxl")
library(readxl)
tun<-read_excel("Tunjai_regeneration.xlsx")
head(tun)

##########
Page 353

tun$Species

##########
Page 354

install.packages("pdftools")
library(pdftools)
pdf_file<-("Wildlife-at-Doi-Inthanon-National-Park.pdf")
pdftext<-pdf_text(pdf_file)

next output to appear on the screen, e.g.
png("My sine graph",bg="wheat",width=600,height=600)
plot(seq(1,10,by=0.01),sin(seq(1,10,by=0.01)),pch=10,col="darkgreen",type="l",lwd=3)
dev.off()

##########
Page 359

layout(matrix(c(1,2,3,4),2,2)) # is equivalent to par(mfrow=c(2,2)) with the plot windows numbered 1 to 4 by row
layout(matrix(c(1,1,2,3),2,2,byrow=TRUE)) # gives a single wide plot at the top (called number 1), and numbers 2 and 3 on the lower row

corners<-matrix(c(0,.6,.6,1,.65,1,0,1,0,.6,0,.65),byrow=TRUE,nrow=3)
split.screen(corners)

screen(1)
plot(0:10,10:0,xlim=c(0,10),ylim=c(0,10),col="red",type="l")
screen(2)
plot(0:10,10:0,xlim=c(0,10),ylim=c(0,10),col="blue",pch=15)
screen(3)
plot(sample(1:100),sample(1:100),xlim=c(0,100),ylim=c(0,100),col=c("blue","green"),pch=0)

##########
Page 360

getwd()

setwd("C:/Users/Rachel/Downloads/wetransfer-126312/DonaldData") # for a PC-user
setwd("/Users/donaldquicke/Downloads") # for a MacOSX-user

##########
Page 361

my_files<-list.files(path="C:/Users/Rachel_R_data")

my_downloads<-list.files("/Users/donaldquicke/Downloads/")

list.files(path=getwd()) # or just list.files(getwd())

##########
Page 362

cat(c(0,6,8,4,3),"\n","egg, caterpillar, pupa, adult","\n",4*21,"\n",letters)

Sys.getlocale()

system('defaults write org.R-project.R force.LANG en_US.UTF-8')

##########
Page 367

pi

2*pi/360

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########

##########