# R code for:
# A protocol for data exploration to avoid common statistical problems
# Methods in Ecology and Evolution
#
# Alain F. Zuur (1,2), Elena N. Ieno (1,2), Chris S. Elphick (3)
#
# 1 Highland Statistics Ltd., 6 Laverock Road, Newburgh, AB41 6FN, UK
# 2 University of Aberdeen, Oceanlab, Main Street, Newburgh, AB41 6AA, UK
# 3 University of Connecticut, Department of Ecology and Evolutionary Biology and Center for Conservation Biology, 75 N. Eagleville Road, U-43, Storrs, CT 06269-3043, USA
#
#
# This file was produced by:
# Alain Zuur (highstat@highstat.com)
# www.highstat.com
#
# A detailed explanation of the R code used in the paper can be found in:
# A Beginner's Guide to R (2009).
# Zuur, AF, Ieno, EN, Meesters, EHWG. Springer
# http://www.springer.com/statistics/computational/book/978-0-387-93836-3
#
#
######################################################################
# DISCLAIMER
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
######################################################################
#
# The R code below was tested on R version 2.9.1
# October 2009
#We assume that the data has been stored in the following directory:
#setwd("D:\\applicat\\HighlandStatistics\\MyPapers\\DataExploration")
#Change this line to where ever you saved the data files.
#Figure 2
Sparrows <- read.table(file = "SparrowsElphick.txt", header = TRUE)
op <- par(no.readonly = TRUE)
par(mfrow= c (1,2), mar = c(5,4,2,1))
boxplot(Sparrows$wingcrd, ylab = "Wing length (mm)")
dotchart(Sparrows$wingcrd, xlab = "Wing length (mm)",
ylab = "Order of the data")
par(op)
#Figure 3
library(lattice)
Z <- cbind(Sparrows$wingcrd, Sparrows$tarsus, Sparrows$head,
Sparrows$culmen, Sparrows$nalospi, Sparrows$wt)
colnames(Z) <- c("wing length", "tarsus length", "head length",
"culmen length", "nalospi to bill tip", "weight")
dotplot(as.matrix(Z), groups = FALSE,
strip = strip.custom(bg = 'white',
par.strip.text = list(cex = 0.8)),
scales = list(x = list(relation = "free"),
y = list(relation = "free"),
draw = FALSE),
col = 1, cex = 0.5, pch = 16,
xlab = "Value of the variable",
ylab = "Order of the data from text file")
#
###############################################################
#Figure 4
#Godwit intake rates
Godwits <- read.table(file="Godwits.txt", header=TRUE)
#Data info:
#Sex
#sex 1=female 0 should go out
#sex 2=male
#Age
#1= adult
#2= juvenile
#0=UNKNOWN
#Location
#locationa=0
#locationb=1
#Period
#period=0 ;southern summr
#period=1; prepare for migration
#period=2; sourthern winter
library(lattice)
Godwits$fSEX <- factor(Godwits$SEX, levels = c(0, 1, 2),
labels = c("Not", "Female", "Male"))
Godwits$fPERIOD <- factor(Godwits$PERIOD, levels = c(0, 1, 2),
labels = c("Summer", "Pre-migration", "Winter"))
bwplot(mgconsumed ~ fPERIOD | fSEX, data = Godwits,
strip = strip.custom(bg = 'white'), subset = SEX!=0,
cex = .5, layout = c(2, 1),
xlab = "Migration period", ylab = "Intake rate",
par.settings = list(
box.rectangle = list(col = 1),
box.umbrella = list(col = 1),
plot.symbol = list(cex = .5, col = 1)),
scales = list(x = list(relation = "same"),
y = list(relation = "same")))
##############################################################
#Figure 5
Sparrows <- read.table(file="SparrowsElphick.txt", header=TRUE)
Sparrows$fMonth<-factor(Sparrows$Month,
levels = c(5, 6, 7, 8, 9, 10),
labels = c("May", "June", "July", "August",
"Sept.", "Oct."))
Sparrows$I1 <- Sparrows$fMonth =="June" |
Sparrows$fMonth =="July" |
Sparrows$fMonth =="August"
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 1))
hist(Sparrows$wt[Sparrows$I1],
xlab = "Weight (g)", breaks = 30,
main = "", ylab = "Frequency")
par(op)
op <- par(no.readonly = TRUE)
library(lattice)
histogram( ~ wt | fMonth, type = "count",
xlab = "Weight (g)",
ylab = "Frequency",
nint=30,layout=c(1,3),
strip.left = strip.custom(bg = 'white'),
strip = F,
col.line = "black", col = "white",
scales = list(x = list(relation = "same"),
y = list(relation = "same"),
draw = TRUE),
subset = fMonth =="June" | fMonth == "July" |fMonth == "August",
data = Sparrows)
par(op)
###############################################################
#Figure 6: Not an R graph
###############################################################
#Figure 7:
RiceField <- read.table(file="ElphickBirdData.txt", header = TRUE)
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 1))
#par(mar = c(4, 4, 3, 2))
plot(table(round(RiceField$AREA * RiceField$AQBIRDS)),
type = "h",
xlim = c(0, 100),
xlab = "Observed values", ylab = "Frequency")
par(op)
###############################################################
#Figure 8
RiceField <- read.table(file="ElphickBirdData.txt", header = TRUE)
#These are all the species
AllS <- c(
"TUSW", "GWFG", "WHGO", "CAGO", "MALL",
"GADW", "GWTE", "CITE", "UNTE", "AMWI", "NOPI",
"NOSH", "RIDU", "CANV", "BUFF", "WODU", "RUDU",
"EUWI", "UNDU", "PBGB", "SORA", "COOT", "COMO",
"AMBI", "BCNH", "GBHE", "SNEG", "GREG", "WFIB",
"SACR", "AMAV", "BNST", "BBPL", "KILL", "LBCU",
"GRYE", "LEYE", "LBDO", "SNIP", "DUNL", "WESA",
"LESA", "PEEP", "RUFF", "UNSH", "RBGU", "HEGU",
"CAGU", "GUSP")
#Determine species richness
Richness <- colSums(RiceField[,AllS] > 0, na.rm = TRUE)
#Remove all covariates
Birds <- RiceField[,AllS]
#To reduce the of variables in the figure, we only used the
#20 species that occured at more than 40 sites.
#As a result, N = 20. Else it becomes a mess.
Birds2 <- Birds[, Richness > 40]
N <- ncol(Birds2)
AllNames <- names(Birds2)
A <- matrix(nrow = N, ncol = N)
for (i in 1:N){
for (j in 1:N){
A[i,j] <- sum(RiceField[,AllS[i]]==0 & RiceField[,AllS[j]]==0, na.rm=TRUE)
}}
A1 <- A/2035
print(A1, digits = 2)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
## [1,] 1.00 0.98 0.98 0.99 0.69 0.96 0.91 0.98 0.99 0.89 0.82 0.83 0.99 1.00 0.99 0.99 1.00 0.99
## [2,] 0.98 0.98 0.97 0.98 0.69 0.95 0.90 0.97 0.98 0.88 0.81 0.82 0.98 0.98 0.98 0.98 0.98 0.98
## [3,] 0.98 0.97 0.98 0.98 0.69 0.95 0.90 0.97 0.98 0.88 0.82 0.82 0.98 0.98 0.98 0.98 0.98 0.98
## [4,] 0.99 0.98 0.98 1.00 0.69 0.96 0.91 0.98 0.99 0.89 0.82 0.83 0.99 1.00 0.99 1.00 1.00 0.99
## [5,] 0.69 0.69 0.69 0.69 0.69 0.69 0.67 0.69 0.69 0.68 0.66 0.65 0.69 0.69 0.69 0.69 0.69 0.69
## [6,] 0.96 0.95 0.95 0.96 0.69 0.96 0.89 0.95 0.96 0.87 0.81 0.82 0.96 0.96 0.96 0.96 0.96 0.96
## [7,] 0.91 0.90 0.90 0.91 0.67 0.89 0.91 0.90 0.91 0.85 0.80 0.81 0.91 0.91 0.91 0.91 0.91 0.91
## [8,] 0.98 0.97 0.97 0.98 0.69 0.95 0.90 0.98 0.98 0.88 0.81 0.83 0.98 0.98 0.98 0.98 0.98 0.98
## [9,] 0.99 0.98 0.98 0.99 0.69 0.96 0.91 0.98 1.00 0.89 0.82 0.83 1.00 1.00 1.00 1.00 1.00 0.99
## [10,] 0.89 0.88 0.88 0.89 0.68 0.87 0.85 0.88 0.89 0.89 0.80 0.80 0.89 0.89 0.89 0.89 0.89 0.89
## [11,] 0.82 0.81 0.82 0.82 0.66 0.81 0.80 0.81 0.82 0.80 0.82 0.77 0.82 0.82 0.82 0.82 0.82 0.82
## [12,] 0.83 0.82 0.82 0.83 0.65 0.82 0.81 0.83 0.83 0.80 0.77 0.83 0.83 0.83 0.83 0.83 0.83 0.83
## [13,] 0.99 0.98 0.98 0.99 0.69 0.96 0.91 0.98 1.00 0.89 0.82 0.83 1.00 1.00 1.00 1.00 1.00 0.99
## [14,] 1.00 0.98 0.98 1.00 0.69 0.96 0.91 0.98 1.00 0.89 0.82 0.83 1.00 1.00 1.00 1.00 1.00 1.00
## [15,] 0.99 0.98 0.98 0.99 0.69 0.96 0.91 0.98 1.00 0.89 0.82 0.83 1.00 1.00 1.00 1.00 1.00 0.99
## [16,] 0.99 0.98 0.98 1.00 0.69 0.96 0.91 0.98 1.00 0.89 0.82 0.83 1.00 1.00 1.00 1.00 1.00 1.00
## [17,] 1.00 0.98 0.98 1.00 0.69 0.96 0.91 0.98 1.00 0.89 0.82 0.83 1.00 1.00 1.00 1.00 1.00 1.00
## [18,] 0.99 0.98 0.98 0.99 0.69 0.96 0.91 0.98 0.99 0.89 0.82 0.83 0.99 1.00 0.99 1.00 1.00 1.00
## [19,] 0.97 0.96 0.96 0.97 0.69 0.94 0.89 0.96 0.97 0.88 0.82 0.82 0.97 0.97 0.97 0.97 0.97 0.97
## [20,] 0.98 0.97 0.97 0.98 0.69 0.95 0.90 0.97 0.98 0.88 0.81 0.83 0.98 0.98 0.98 0.98 0.99 0.98
## [,19] [,20]
## [1,] 0.97 0.98
## [2,] 0.96 0.97
## [3,] 0.96 0.97
## [4,] 0.97 0.98
## [5,] 0.69 0.69
## [6,] 0.94 0.95
## [7,] 0.89 0.90
## [8,] 0.96 0.97
## [9,] 0.97 0.98
## [10,] 0.88 0.88
## [11,] 0.82 0.81
## [12,] 0.82 0.83
## [13,] 0.97 0.98
## [14,] 0.97 0.98
## [15,] 0.97 0.98
## [16,] 0.97 0.98
## [17,] 0.97 0.99
## [18,] 0.97 0.98
## [19,] 0.97 0.96
## [20,] 0.96 0.99
rownames(A1) <- AllNames
colnames(A1) <- AllNames
library(lattice)
panel.corrgram.2 <- function(x, y, z, subscripts, at = pretty(z), scale = 0.8, ...)
{
require("grid", quietly = TRUE)
x <- as.numeric(x)[subscripts]
y <- as.numeric(y)[subscripts]
z <- as.numeric(z)[subscripts]
zcol <- level.colors(z, at = at, ...)
for (i in seq(along = z))
{
lims <- range(0, z[i])
tval <- 2 * base::pi *
seq(from = lims[1], to = lims[2], by = 0.01)
grid.circle(x = x[i], y = y[i], r = .5 * scale,
default.units = "native")
grid.polygon(x = x[i] + .5 * scale * c(0, sin(tval)),
y = y[i] + .5 * scale * c(0, cos(tval)),
default.units = "native",
gp = gpar(fill = zcol[i]))
}
}
levelplot(A1,xlab=NULL,ylab=NULL,
at=do.breaks(c(0.5,1.01),101),
panel=panel.corrgram.2,
scales=list(x=list(rot=90)),
colorkey=list(space="top"),
col.regions=colorRampPalette(c("red","white","blue")))
#Grey colours
levelplot(A1,xlab=NULL,ylab=NULL,
at=do.breaks(c(0.5,1.01),101),
panel=panel.corrgram.2,
scales=list(x=list(rot=90)),
colorkey=list(space="top"),
col.regions=colorRampPalette(c(grey(0.8),grey(0.5),grey(0.2))))
####################################################################
#Figure 9 & Table 1
Sparrows2 <- read.table(file = "VegSamplesV1.txt", header = TRUE)
#Different Sparrow object
names(Sparrows2)
## [1] "Year" "Site" "UniversalPlotName" "Banded"
## [5] "PtCountsum" "Avgmaxht" "Avgdens" "ht.thatch"
## [9] "S.patens" "Distichlis" "S.alternifloraShort" "S.alternifloraTall"
## [13] "Juncus" "Bare" "Other" "Phragmites"
## [17] "Shrub" "Tallsedge" "Water"
# [1] "Year" "Site" "UniversalPlotName"
# [4] "Banded" "PtCountsum" "Avgmaxht"
# [7] "Avgdens" "ht.thatch" "S.patens"
#[10] "Distichlis" "S.alternifloraShort" "S.alternifloraTall"
#[13] "Juncus" "Bare" "Other"
#[16] "Phragmites" "Shrub" "Tallsedge"
#[19] "Water"
#Load our own library files
source("HighstatLib.r")
#Seclect covraites
Z<-Sparrows2[,c("Avgmaxht", "Avgdens", "ht.thatch",
"S.patens", "Distichlis", "S.alternifloraShort",
"S.alternifloraTall", "Juncus", "Bare", "Other",
"Phragmites", "Shrub", "Tallsedge", "Water")]
corvif(Z) #Part of Table 1
## Correlations of the variables
##
## Avgmaxht Avgdens ht.thatch S.patens Distichlis S.alternifloraShort
## Avgmaxht 1.00000000 -0.22841524 0.466218046 -0.17785017 -0.02337749 -0.44042742
## Avgdens -0.22841524 1.00000000 0.062575238 0.67411203 0.13160006 -0.35875821
## ht.thatch 0.46621805 0.06257524 1.000000000 0.16053731 0.02232166 -0.40140786
## S.patens -0.17785017 0.67411203 0.160537307 1.00000000 -0.10290139 -0.41441467
## Distichlis -0.02337749 0.13160006 0.022321660 -0.10290139 1.00000000 -0.17327011
## S.alternifloraShort -0.44042742 -0.35875821 -0.401407857 -0.41441467 -0.17327011 1.00000000
## S.alternifloraTall 0.60074164 -0.51803067 0.254760430 -0.55521017 -0.21907684 -0.15659005
## Juncus 0.17293819 0.32490496 -0.003660471 0.03727048 0.04757495 -0.31812489
## Bare -0.23941028 -0.20344863 -0.155596279 -0.05769548 -0.24631400 -0.03848335
## Other -0.36157539 0.13420907 -0.194936964 0.06081853 -0.12816334 0.08430546
## Phragmites 0.38995648 -0.06910767 0.330766492 0.08310766 0.08813589 -0.32519571
## Shrub -0.07595686 -0.19996137 -0.141009710 -0.16031227 0.21464963 0.16155126
## Tallsedge 0.28552986 -0.07218760 -0.066039061 0.07970131 -0.09316033 -0.15043833
## Water -0.25038911 -0.12840907 -0.120468982 -0.17287729 -0.21243584 0.06676282
## S.alternifloraTall Juncus Bare Other Phragmites Shrub
## Avgmaxht 0.600741637 0.172938189 -0.23941028 -0.36157539 0.38995648 -0.07595686
## Avgdens -0.518030668 0.324904964 -0.20344863 0.13420907 -0.06910767 -0.19996137
## ht.thatch 0.254760430 -0.003660471 -0.15559628 -0.19493696 0.33076649 -0.14100971
## S.patens -0.555210170 0.037270482 -0.05769548 0.06081853 0.08310766 -0.16031227
## Distichlis -0.219076843 0.047574949 -0.24631400 -0.12816334 0.08813589 0.21464963
## S.alternifloraShort -0.156590051 -0.318124889 -0.03848335 0.08430546 -0.32519571 0.16155126
## S.alternifloraTall 1.000000000 -0.231319555 -0.01242215 -0.20152578 0.02700616 -0.17308212
## Juncus -0.231319555 1.000000000 -0.14345423 -0.11858986 0.08711016 0.04511243
## Bare -0.012422149 -0.143454233 1.00000000 0.18374555 -0.05895322 -0.09935402
## Other -0.201525783 -0.118589860 0.18374555 1.00000000 -0.12488041 0.02376765
## Phragmites 0.027006156 0.087110160 -0.05895322 -0.12488041 1.00000000 0.09161061
## Shrub -0.173082116 0.045112426 -0.09935402 0.02376765 0.09161061 1.00000000
## Tallsedge -0.050009373 -0.070738224 0.14131312 0.26940954 0.09853599 -0.04012283
## Water 0.002440565 -0.124503346 0.12718566 -0.07270756 -0.11213772 -0.09872519
## Tallsedge Water
## Avgmaxht 0.28552986 -0.250389111
## Avgdens -0.07218760 -0.128409073
## ht.thatch -0.06603906 -0.120468982
## S.patens 0.07970131 -0.172877290
## Distichlis -0.09316033 -0.212435845
## S.alternifloraShort -0.15043833 0.066762825
## S.alternifloraTall -0.05000937 0.002440565
## Juncus -0.07073822 -0.124503346
## Bare 0.14131312 0.127185661
## Other 0.26940954 -0.072707562
## Phragmites 0.09853599 -0.112137724
## Shrub -0.04012283 -0.098725195
## Tallsedge 1.00000000 -0.059059880
## Water -0.05905988 1.000000000
##
##
## Variance inflation factors
## Warning in summary.lm(object): essentially perfect fit: summary may be unreliable
## GVIF
## Avgmaxht 6.120018
## Avgdens 3.206401
## ht.thatch 1.671224
## S.patens 159.350658
## Distichlis 53.754540
## S.alternifloraShort 121.463782
## S.alternifloraTall 159.382860
## Juncus 44.995377
## Bare 12.058665
## Other 5.817015
## Phragmites 3.749027
## Shrub 2.781804
## Tallsedge 4.409398
## Water 17.067777
#Run linear regression
M1<-lm(Banded~Avgmaxht + Avgdens + ht.thatch + S.patens +
Distichlis + S.alternifloraShort + S.alternifloraTall +
Juncus + Bare + Other + Phragmites + Shrub + Tallsedge +
Water, data = Sparrows2)
summary(M1) #Part of Table 1
##
## Call:
## lm(formula = Banded ~ Avgmaxht + Avgdens + ht.thatch + S.patens +
## Distichlis + S.alternifloraShort + S.alternifloraTall + Juncus +
## Bare + Other + Phragmites + Shrub + Tallsedge + Water, data = Sparrows2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.142 -5.431 1.011 5.240 18.634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.612e+02 8.780e+01 -1.836 0.0730 .
## Avgmaxht 4.195e-01 3.548e-01 1.183 0.2432
## Avgdens 6.158e-02 1.719e-01 0.358 0.7219
## ht.thatch 8.427e-04 6.241e-01 0.001 0.9989
## S.patens 1.626e+00 8.562e-01 1.899 0.0640 .
## Distichlis 1.703e+00 8.558e-01 1.990 0.0527 .
## S.alternifloraShort 1.668e+00 8.462e-01 1.971 0.0549 .
## S.alternifloraTall 1.449e+00 8.521e-01 1.700 0.0960 .
## Juncus 2.009e+00 8.354e-01 2.405 0.0203 *
## Bare 1.620e+00 8.615e-01 1.880 0.0666 .
## Other 1.800e+00 9.807e-01 1.836 0.0730 .
## Phragmites 1.851e+00 1.003e+00 1.846 0.0715 .
## Shrub 6.300e-02 1.248e+00 0.050 0.9600
## Tallsedge 1.471e+00 1.172e+00 1.255 0.2160
## Water 2.092e+00 1.070e+00 1.955 0.0568 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.45 on 45 degrees of freedom
## Multiple R-squared: 0.393, Adjusted R-squared: 0.2041
## F-statistic: 2.081 on 14 and 45 DF, p-value: 0.03214
#Chop out covariates
Z<-Sparrows2[,c("ht.thatch", "S.patens", "Distichlis",
"S.alternifloraShort", "Juncus", "Bare", "Other",
"Phragmites", "Shrub", "Tallsedge", "Water")]
corvif(Z) #Part of Table 1
## Correlations of the variables
##
## ht.thatch S.patens Distichlis S.alternifloraShort Juncus Bare
## ht.thatch 1.000000000 0.16053731 0.02232166 -0.40140786 -0.003660471 -0.15559628
## S.patens 0.160537307 1.00000000 -0.10290139 -0.41441467 0.037270482 -0.05769548
## Distichlis 0.022321660 -0.10290139 1.00000000 -0.17327011 0.047574949 -0.24631400
## S.alternifloraShort -0.401407857 -0.41441467 -0.17327011 1.00000000 -0.318124889 -0.03848335
## Juncus -0.003660471 0.03727048 0.04757495 -0.31812489 1.000000000 -0.14345423
## Bare -0.155596279 -0.05769548 -0.24631400 -0.03848335 -0.143454233 1.00000000
## Other -0.194936964 0.06081853 -0.12816334 0.08430546 -0.118589860 0.18374555
## Phragmites 0.330766492 0.08310766 0.08813589 -0.32519571 0.087110160 -0.05895322
## Shrub -0.141009710 -0.16031227 0.21464963 0.16155126 0.045112426 -0.09935402
## Tallsedge -0.066039061 0.07970131 -0.09316033 -0.15043833 -0.070738224 0.14131312
## Water -0.120468982 -0.17287729 -0.21243584 0.06676282 -0.124503346 0.12718566
## Other Phragmites Shrub Tallsedge Water
## ht.thatch -0.19493696 0.33076649 -0.14100971 -0.06603906 -0.12046898
## S.patens 0.06081853 0.08310766 -0.16031227 0.07970131 -0.17287729
## Distichlis -0.12816334 0.08813589 0.21464963 -0.09316033 -0.21243584
## S.alternifloraShort 0.08430546 -0.32519571 0.16155126 -0.15043833 0.06676282
## Juncus -0.11858986 0.08711016 0.04511243 -0.07073822 -0.12450335
## Bare 0.18374555 -0.05895322 -0.09935402 0.14131312 0.12718566
## Other 1.00000000 -0.12488041 0.02376765 0.26940954 -0.07270756
## Phragmites -0.12488041 1.00000000 0.09161061 0.09853599 -0.11213772
## Shrub 0.02376765 0.09161061 1.00000000 -0.04012283 -0.09872519
## Tallsedge 0.26940954 0.09853599 -0.04012283 1.00000000 -0.05905988
## Water -0.07270756 -0.11213772 -0.09872519 -0.05905988 1.00000000
##
##
## Variance inflation factors
## Warning in summary.lm(object): essentially perfect fit: summary may be unreliable
## GVIF
## ht.thatch 1.483803
## S.patens 1.405993
## Distichlis 1.355857
## S.alternifloraShort 2.131070
## Juncus 1.300832
## Bare 1.231097
## Other 1.189983
## Phragmites 1.250227
## Shrub 1.156870
## Tallsedge 1.182988
## Water 1.186686
#Linear regression on subset
M2 <- lm(Banded ~ ht.thatch + S.patens +
Distichlis + S.alternifloraShort +
Juncus + Bare + Other + Phragmites + Shrub + Tallsedge +
Water, data = Sparrows2)
summary(M2) #Part of Table 1
##
## Call:
## lm(formula = Banded ~ ht.thatch + S.patens + Distichlis + S.alternifloraShort +
## Juncus + Bare + Other + Phragmites + Shrub + Tallsedge +
## Water, data = Sparrows2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.1448 -6.3075 0.4815 7.2374 20.0593
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.93887 8.38424 1.185 0.241685
## ht.thatch 0.13163 0.59680 0.221 0.826366
## S.patens 0.08011 0.08162 0.982 0.331243
## Distichlis 0.15929 0.13794 1.155 0.253875
## S.alternifloraShort 0.12045 0.11375 1.059 0.294972
## Juncus 0.59087 0.14415 4.099 0.000159 ***
## Bare 0.03854 0.27935 0.138 0.890862
## Other 0.03050 0.45018 0.068 0.946264
## Phragmites 0.65102 0.58768 1.108 0.273473
## Shrub -1.59457 0.81712 -1.951 0.056853 .
## Tallsedge 0.48892 0.61606 0.794 0.431324
## Water 0.11326 0.28638 0.395 0.694233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.61 on 48 degrees of freedom
## Multiple R-squared: 0.3331, Adjusted R-squared: 0.1802
## F-statistic: 2.179 on 11 and 48 DF, p-value: 0.03169
M2 <- lm(Banded ~ Juncus + Shrub, data = Sparrows2)
drop1(M2, test = "F")
## Single term deletions
##
## Model:
## Banded ~ Juncus + Shrub
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 5815.2 280.43
## Juncus 1 2015.82 7831.0 296.29 19.7589 4.11e-05 ***
## Shrub 1 341.07 6156.3 281.85 3.3431 0.07272 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(M1)
## (Intercept) Avgmaxht Avgdens ht.thatch S.patens
## -1.611578e+02 4.195265e-01 6.157904e-02 8.426682e-04 1.625527e+00
## Distichlis S.alternifloraShort S.alternifloraTall Juncus Bare
## 1.702737e+00 1.668078e+00 1.448936e+00 2.009191e+00 1.619609e+00
## Other Phragmites Shrub Tallsedge Water
## 1.800322e+00 1.850829e+00 6.300294e-02 1.470594e+00 2.091850e+00
step(M2)
## Start: AIC=280.43
## Banded ~ Juncus + Shrub
##
## Df Sum of Sq RSS AIC
## <none> 5815.2 280.43
## - Shrub 1 341.07 6156.3 281.85
## - Juncus 1 2015.82 7831.0 296.29
##
## Call:
## lm(formula = Banded ~ Juncus + Shrub, data = Sparrows2)
##
## Coefficients:
## (Intercept) Juncus Shrub
## 19.0220 0.5354 -1.3237
M3 <- lm(Banded ~ Juncus+Shrub, data = Sparrows2)
summary(M3) #Part of Table 1
##
## Call:
## lm(formula = Banded ~ Juncus + Shrub, data = Sparrows2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.022 -7.495 1.283 7.978 20.978
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.0220 1.5944 11.930 < 2e-16 ***
## Juncus 0.5354 0.1204 4.445 4.11e-05 ***
## Shrub -1.3237 0.7240 -1.828 0.0727 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.1 on 57 degrees of freedom
## Multiple R-squared: 0.2822, Adjusted R-squared: 0.2571
## F-statistic: 11.21 on 2 and 57 DF, p-value: 7.858e-05
#Figure 9
Z <- as.vector(as.matrix(Sparrows2[, c("Avgmaxht", "Avgdens",
"ht.thatch", "S.patens", "Distichlis",
"S.alternifloraShort", "S.alternifloraTall", "Juncus",
"Bare", "Other", "Phragmites", "Shrub", "Tallsedge", "Water")]))
#Setup the data in vector format for the xyplot
Y10 <- rep(Sparrows2$Banded, 14)
MyNames <- names(Sparrows2[,c("Avgmaxht", "Avgdens", "ht.thatch",
"S.patens", "Distichlis", "S.alternifloraShort",
"S.alternifloraTall", "Juncus", "Bare", "Other",
"Phragmites", "Shrub", "Tallsedge", "Water")])
ID10 <- rep(MyNames, each = length(Sparrows2$Banded))
library(lattice)
ID11 <- factor(ID10, labels = c("% Juncus gerardii",
"% Shrub", "Height of thatch", "% Spartina patens",
"% Distichlis", "% Bare ground", "% Other vegetation",
"% Phragmites australis", "% Tall sedge", "% Water",
"% Spartina alterniflora (short)",
"% Spartina alterniflora (tall)",
"Maximum vegetation height",
"Vegetation stem density"),
levels = c("Juncus", "Shrub", "Avgmaxht", "S.patens",
"Distichlis", "Bare", "Other", "Phragmites",
"Tallsedge", "Water", "S.alternifloraShort",
"S.alternifloraTall", "ht.thatch", "Avgdens"))
xyplot(Y10 ~ Z | ID11, col = 1,
strip = function(bg='white',...) strip.default(bg='white',...),
scales = list(alternating = T,
x = list(relation = "free"),
y = list(relation = "same")),
xlab = "Covariates",
par.strip.text = list(cex = 0.8),
ylab = "Banded",
panel=function(x, y, subscripts,...){
panel.grid(h =- 1, v = 2)
panel.points(x, y, col = 1, pch = 16)
if(ID10[subscripts][1] != "Tallsedge") {panel.loess(x,y,col=1,lwd=2)}
})
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0022278
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0022278
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0022278
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0022278
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0022278
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.0472
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0049
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0049
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0049
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0049
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.0049
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.012343
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.012343
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.012343
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.012343
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : radius 0.012343
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : all data on boundary of
## neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : pseudoinverse used at -0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : neighborhood radius 0.1111
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : reciprocal condition number 1
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, : zero-width neighborhood. make
## span bigger
##################################################################
#Figure 10
Sparrows <- read.table(file = "SparrowsElphick.txt", header = TRUE)
source(file = "HighstatLib.r")
MyNames <- c("wing chord", "tarsus length", "head length",
"culmen length", "nalospi to bill tip", "weightt")
pairs(Sparrows[,c(1, 3, 4, 5, 6, 7)],
lower.panel = panel.cor,
cex.labels=1.3,
labels=MyNames)
###################################################################
#Figure 11
Sparrows <- read.table(file = "SparrowsElphick.txt", header = TRUE)
#Take the data from species 1, Sex = 0 and Wing length >= 65
I1 <- Sparrows$SpeciesCode == 1 & Sparrows$Sex != "0" & Sparrows$wingcrd < 65
Wing1<- Sparrows$wingcrd[I1]
Wei1 <- Sparrows$wt[I1]
Mon1 <- factor(Sparrows$Month[I1])
Sex1<- factor(Sparrows$Sex[I1])
#Define Month and Sex as categorical variables
fMonth1 <- factor(Mon1,levels=c(5,6,7,8,9),
labels=c("May","Jun","Jul","Aug","Sep"))
fSex1 <- factor(Sex1, levels=c(4,5),labels=c("Male","Female"))
M1 <- lm(Wei1 ~ Wing1*fMonth1*fSex1)
summary(M1)
##
## Call:
## lm(formula = Wei1 ~ Wing1 * fMonth1 * fSex1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4629 -0.6407 -0.0633 0.5371 6.4472
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.5335 16.5746 2.446 0.0147 *
## Wing1 -0.3363 0.2817 -1.194 0.2329
## fMonth1Jun -34.5083 16.7505 -2.060 0.0397 *
## fMonth1Jul -37.5463 16.9226 -2.219 0.0268 *
## fMonth1Aug -45.8166 17.8946 -2.560 0.0106 *
## fMonth1Sep -14.1815 40.0069 -0.354 0.7231
## fSex1Female -25.9335 38.8650 -0.667 0.5048
## Wing1:fMonth1Jun 0.5808 0.2847 2.040 0.0416 *
## Wing1:fMonth1Jul 0.6324 0.2877 2.198 0.0282 *
## Wing1:fMonth1Aug 0.7757 0.3050 2.543 0.0112 *
## Wing1:fMonth1Sep 0.1963 0.7159 0.274 0.7840
## Wing1:fSex1Female 0.4529 0.6931 0.653 0.5136
## fMonth1Jun:fSex1Female 23.7940 39.2265 0.607 0.5443
## fMonth1Jul:fSex1Female 42.7893 39.3326 1.088 0.2769
## fMonth1Aug:fSex1Female 42.3729 40.2315 1.053 0.2925
## fMonth1Sep:fSex1Female NA NA NA NA
## Wing1:fMonth1Jun:fSex1Female -0.4210 0.6996 -0.602 0.5475
## Wing1:fMonth1Jul:fSex1Female -0.7792 0.7015 -1.111 0.2670
## Wing1:fMonth1Aug:fSex1Female -0.7752 0.7175 -1.080 0.2802
## Wing1:fMonth1Sep:fSex1Female NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.097 on 887 degrees of freedom
## Multiple R-squared: 0.4038, Adjusted R-squared: 0.3923
## F-statistic: 35.33 on 17 and 887 DF, p-value: < 2.2e-16
anova(M1)
## Analysis of Variance Table
##
## Response: Wei1
## Df Sum Sq Mean Sq F value Pr(>F)
## Wing1 1 531.39 531.39 441.6123 < 2.2e-16 ***
## fMonth1 4 46.50 11.63 9.6621 1.191e-07 ***
## fSex1 1 64.08 64.08 53.2578 6.495e-13 ***
## Wing1:fMonth1 4 39.18 9.80 8.1405 1.900e-06 ***
## Wing1:fSex1 1 5.08 5.08 4.2228 0.04018 *
## fMonth1:fSex1 3 26.96 8.99 7.4698 6.090e-05 ***
## Wing1:fMonth1:fSex1 3 9.58 3.19 2.6548 0.04741 *
## Residuals 887 1067.31 1.20
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Make the coplot
coplot(Wei1 ~ Wing1 | fMonth1 * fSex1, ylab = "Weight (g)",
xlab = "Wing length (mm)",
panel = function(x, y, ...) {
tmp <- lm(y ~ x, na.action = na.omit)
abline(tmp)
points(x, y) })
##################################################################
#Figure 12
Waders <- read.table(file = "wader.txt", header = TRUE)
#Define the time axis
Time <- seq(1,25)
par(mfrow = c(2, 2), mar = c(5, 4, 3, 2))
plot(Time, Waders$C.fuscicolis, type = "l", xlab = "Time (2 weeks)",
ylab = "C. fuscicollis abundance")
acf(Waders$C.fuscicolis, main = "C. fuscicollis ACF")
plot(Time, Waders$L.dominicanus, type = "l", xlab = "Time (2 weeks)",
ylab = "L. dominicanus abundance")
acf(Waders$L.dominicanus, main = "L. dominicanus ACF")
#################################################################