Showing posts with label sub(). Show all posts
Showing posts with label sub(). Show all posts

Sunday, September 20, 2015

Import dbf to R, Manipulate Strings with grep & sub Function

Here's a set of historical species presence records of a certain geographical region (data-link). I wanted to manipulate / simplify strings (species names) and get an overview of the data. 
...The tasks were to split genera and epitheta, to exclude species with specific strings included and to get rid of unwanted text (author names). For graphical presentation of the species record history I did a plot with segments indicating the first and last year of a species record:



 



# for dbf import you'll need:
require(foreign)

# change to your directory: TLMFB_IL <- read.dbf("D:\\Downloads\\TLMFB_IL.dbf")

str(TLMFB_IL)

# get rid of hybrids:
dat_sub1 <- TLMFB_IL[-grep(" x ", TLMFB_IL$NAME),]

# check how many were dismissed:
length(TLMFB_IL$NAME) - length(dat_sub1$NAME)

# get genera and epitheta:
Gen <- sub(" .*", "", dat_sub1$NAME)
Epi <- sub(" .*", "", substring(dat_sub1$NAME, nchar(Gen)+2))

# get rid of species with unsure determination:
dat_sub2 <- dat_sub1[Epi != "spec." &
Epi != "sp." &
Epi != "cf." &
Epi != "" &
Epi != "?", ]

length(TLMFB_IL$NAME) - length(dat_sub2$NAME)

# check:
dat_sub2$NAME
TLMFB_IL$NAME

table(Epi != "spec." &
Epi != "sp." &
Epi != "cf." &
Epi != "" &
Epi != "?")

length(grep(" x ", TLMFB_IL$NAME))
length(dat_sub2$NAME)

# get genera and epitheta:
Gen <- sub(" .*", "", dat_sub2$NAME)
Epi <- sub(" .*", "", substring(dat_sub2$NAME, nchar(Gen)+2))

# get rid of authors:
sp <- paste(Gen, Epi)
length(sp)

# check arbitrary sample of 100 rows, window should be large

# enough to show columns next to each other:
id <- sample(1:3175, 100)
data.frame(sp = sp, orig = dat_sub2$NAME)[id,]

# add species names without authors to dataframe
dat_sub2$Sp <- sp

# there are some erronous values that should be discarded:
dat_sub3 <- dat_sub2[dat_sub2$date_long < 2010,]
str(dat_sub3)

# get max and min year at which species were recorded
y_min <- aggregate(. ~ Sp, min, data = dat_sub3[,c("Sp","date_long")])
y_max <- aggregate(. ~ Sp, max, data = dat_sub3[,c("Sp","date_long")])

# plot each species first and last record in line plot, data:
head(pldat <- data.frame(Sp = y_min[,1], y_min = y_min[,2],
y_max = y_max[,2], span = y_max[,2] - y_min[,2]))

# plot:
# example("segments")

# new ordering for plot:
pldat <- pldat[order(pldat$y_min, pldat$span),]
plot(x = c(min(pldat$y_min), max(pldat$y_max)), y = c(1, nrow(pldat)),
type = 'n', xlab = "Year", ylab ="", axes = F, frame.plot=FALSE)

segments(x0 = pldat$y_min, x1 = pldat$y_max, y0 = 1:nrow(pldat),
y1 = 1:nrow(pldat), col = "gray60")

axis(1, pretty(1800:2000), cex.axis = 0.75)
mtext(paste("Species Records Innsbruck\n", " (Sp - N = ", nrow(pldat), ")",
sep = ""), side = 2, line = -2)


# especially inbetween the 60's and 80's many species were
# re-recorded and newly added by only a few authors:

sixt_eight <- dat_sub3[dat_sub3$date_long > 1960 & dat_sub3$date_long < 1980, c("AUTOR", "Sp")]

data.frame(table(as.character(sixt_eight$AUTOR)))

To cite package ‘foreign’ in publications use:
R-core members, Saikat DebRoy , Roger Bivand and
others: see COPYRIGHTS file in the sources. (2011). foreign: Read Data Stored by Minitab, S,
SAS, SPSS, Stata, Systat, dBase, .... R package version 0.8-44.
http://CRAN.R-project.org/package=foreign
Read more »

Use of Classification Trees to Investigate Traits of Invasive Species

Which traits make an alien species invasive?
Due to what traits an alien species becomes established in a foreign flora?


This kind of questions could be analysed by the use of recursive partitioning and classification trees..
(the below example also includes some useful data manipulation techniques)...

data-download*

*data is courtesy of BiolFlor:
Klotz, S., Kühn, I. & Durka, W. [Hrsg.] (2002): BIOLFLOR - Eine Datenbank zu biologisch-ökologischen Merkmalen der Gefäßpflanzen in Deutschland. - Schriftenreihe für Vegetationskunde 38. Bonn: Bundesamt für Naturschutz.

library(coin)
library(zoo)
library(modeltools)
library(sandwich)
library(strucchange)
library(vcd)
library(colorspace)
library(Formula)
library(partykit)

## (!)note: package party should not be loaded at the same time as partykit

## read data, mind to change the path:
dat <- read.csv("E:\\R\\Data\\Alien_Plant_Traits.csv",
sep = ";")

## delete unneeded variable:
dat$Repr_Literatur <- NULL
dat$Diasporen..bzw..Germinulentyp <- NULL
dat$Heteromorphietyp <- NULL
dat$Ausprägung <- NULL
dat$Status <- NULL


## source data uses "," as decimal seperator, as R uses
## "." we will change this appropriately:
dat$Länge_MW_mm <- as.numeric(sub(",", ".", dat$Länge_MW_mm))

## explore dataset:
head(dat[,1:10])
str(dat)

# Name : species names
# Status : floristic status; "Neophyt" = alien species, which is the whole set
# Einführung : mode of introduction
# Einbürg_grad_Max : level of naturalization
# Lebensform : lifeform
# Ausprägung : manifestation of lifeform, "immmer" = always, "manchmal" = sometimes
# Blattanatomie : leaf anatomy
# Blühdauer_Mon : flowering duration
# N_Blühphasen : number of flowering phases
# Länge_MW_mm : length of diaspore / germinule, "MW_mm" = mean value in milimeters
# Länge_Bezug : "Germinule" = length given is that of germinule,
# "Diaspore" = that of diaspore
# Länge_Bezug_Einschr : restricted knowledge, maybe erroneus data, "_FALSCH" = false, "_WAHR" = true
# Reproduktionstyp : reproduction type, "Samen" = seed, "vegetative" = clonal proliferation,
# "Samen/Sporen" = both

## ad diaspores / germinules:
## there are 2 cases - (1) a diaspore (= unit of propagation)
## is at the same time the germinule (= unit of seedling germination).
## this is the case when single seeds are propagated.
## (2) the diaspores disintegrates into many germinules.
## this is the cases when fruits or parts of them
## which contain many seeds are propagated.

## length measurments ("Länge_MW_mm") were taken for either
## germinules or diaspores or both (see "Länge_Bezug") and sometimes
## data may be erronous as indicated by the variable "Länge_Bezug_Einschr"
## what i want is the surely correct length-values of each germinules and
## diaspores for each species. therefore i create new variables:

dat1 <-
reshape(direction = "wide", idvar = "Name",
timevar = "Länge_Bezug_Einschr",
v.names = "Länge_MW_mm",
drop = "Länge_Bezug",
data = dat)


## there's a warning that tells me i aggregated several variables that
## actually were varying - i solve this by pasting those variables
## into one character string per each species
## (beforehand i have to convert them to character vectors)
dat$Einführung <- as.character(dat$Einführung)
dat$Lebensform <- as.character(dat$Lebensform)
dat$Blattanatomie <- as.character(dat$Blattanatomie)
str(dat)

## make function for pasting multiple character strings of each species
## i want only unique values to be pasted, seperated by a comma:
my.paste <- function(x) {paste(unique(x), collapse = ", ")}
test.vector <- c("a", "a", "b", "c")
my.paste(test.vector)

print(dat2 <-
aggregate(cbind(dat$Einführung, dat$Lebensform, dat$Blattanatomie) ~
dat$Name, FUN = my.paste))

## reset column names:
colnames(dat2) <- c("Name", "Einführung", "Lebensform", "Blattanatomie")
head(dat2)
## re-bind dataframes, variables with varying entries per species
## falsely aggregated in dat1 will be discarded and instead the
## pasted character strings of dat2, containing all info will be used.
## also variables present in both sub-datasets will be discarded:
head(dat2)

dat3 <- data.frame(dat1[, - which(names(dat1)%in%
c("Einführung", "Lebensform", "Blattanatomie"))],
dat2)
head(dat3)
## check for consistency (should be 340 species):
## Name-variables of the two sub-datasets should match
head(dat3)
dim(dat3)
dat3$Name == dat3$Name.1


## now i will choose one of the length measurements in the
## dataset - actually i'm interested in propagation (shorter/lighter propagules
## will be distributed easierly than larger/heavy ones) - that said i want
## diaspore weights but i don't have the data for each species.
## as an escape i will use the largest value within the different variables.
## the logic behind this would be that if diaspore = germinule
## the lengths will be the same and a germinule as a part of a diaspore
## will never be longer.
## by choosing the largest value i thus only use the germinule's value
## if the diaspore's value does not exist. this is seldomly the case and the
## length class will very likely still be the same as it would be for the
## diaspore of that plant, if i had the data
## (subsequently i will classify lengths):

dat3$Length_Diasp_new <- rep(NA, nrow(dat3))
dat3[is.na(dat3)] <- 0
for(i in 1:nrow(dat3)) {dat3$Length_Diasp_new[i] <- max(dat3[i, c(6:9)])}

## classify by cut-function:
dat3$Length_Class <- as.ordered(cut(dat3$Length_Diasp_new, c(0, 2, 5, Inf), labels = FALSE))

## i have to convert some character strings to factors
str(dat3)
dat3$Einführung <- as.factor(dat3$Einführung)
dat3$Lebensform <- as.factor(dat3$Lebensform)
dat3$Blattanatomie <- as.factor(dat3$Blattanatomie)
str(dat3)

## number if factor levels is too high for some variables -
## good results will rather be yielded
## if there are not too many levels:

levels(dat3$Lebensform)
## (1) makrophanerophytes will be "Tree"
## (2) chamaephyts, hemiphanerophyts, nanophanerophyts &
## pseudophanerophyts will be "Shrub"
## (3) therophytes will be "Terr_Herbs_ann" (annual herbs)
## (4) geophytes and hemikryptophytes will be "Terr_Herb_per" (perennial herbs)
## (5) hydrophyts will be "Hydr_Herb"

## sometimes the classification in the original data is unclear -
## i will assign distinct lifeforms following the subsequent rules:
## "Chamaephyt, Therophyt" --> "Terr_Herb_ann"
## "Hemikryptophyt, Therophyt" --> "Terr_Herb_ann"
## "Geophyt, Hemikryptophyt" --> "Terr_Herb_per"
## "Hydrophyt, Geophyt" --> "Hydr_Herb"
## "Hydrophyt, Hemikryptophyt" --> "Hydr_Herb"
## "Makrophanerophyt, Nanophanerophyt" --> "Tree"

dat3$Lifeform <- rep(NA, nrow(dat3))

dat3$Lifeform[dat3$Lebensform == "Makrophanerophyt"] <- "Tree"
dat3$Lifeform[dat3$Lebensform == "Pseudophanerophyt"|
dat3$Lebensform == "Nanophanerophyt"|
dat3$Lebensform == "Hemiphanerophyt"|
dat3$Lebensform == "Chamaephyt"] <- "Scrub"
dat3$Lifeform[dat3$Lebensform == "Therophyt"] <- "Terr_Herb_ann"
dat3$Lifeform[dat3$Lebensform == "Geophyt"|
dat3$Lebensform == "Hemikryptophyt"] <- "Terr_Herb_per"
dat3$Lifeform[dat3$Lebensform == "Hydrophyt"] <- "Hydr_Herb"

dat3$Lifeform[dat3$Lebensform == "Chamaephyt, Therophyt"] <- "Terr_Herb_ann"
dat3$Lifeform[dat3$Lebensform == "Hemikryptophyt, Therophyt"] <- "Terr_Herb_ann"
dat3$Lifeform[dat3$Lebensform == "Geophyt, Hemikryptophyt"] <- "Terr_Herb_per"
dat3$Lifeform[dat3$Lebensform == "Hydrophyt, Geophyt"] <- "Hydr_Herb"
dat3$Lifeform[dat3$Lebensform == "Hydrophyt, Hemikryptophyt"] <- "Hydr_Herb"
dat3$Lifeform[dat3$Lebensform == "Makrophanerophyt, Nanophanerophyt"] <- "Tree"

## convert to factor:
dat3$Lifeform <- as.factor(dat3$Lifeform)
str(dat3)

levels(dat3$Einführung)
## "Begleiter"; "Begleiter, Verwilderte Nutzpflanze";
## "Begleiter, Verwilderte Zierpflanze"; "Transportbegleiter";
## "Saatgutbegleiter"; "Saatgutbegleiter, Verwilderte Nutzpflanze";
## "Saatgutbegleiter, Verwilderte Zierpflanze" <-- "accompanying"
## "spontan"; "unbekannt"; "spontan, Verwilderte Nutzpflanze" <-- "spontaneous"
## "Verwilderte Nutzpflanze"; "Verwilderte Nutzpflanze, Verwilderte Zierpflanze";
## "Verwilderte Zierpflanze" <-- "cultivated"

dat3$Introduction <- rep(NA, nrow(dat3))
dat3$Introduction[dat3$Einführung == "Begleiter"|
dat3$Einführung == "Begleiter, Verwilderte Nutzpflanze"|
dat3$Einführung == "Begleiter, Verwilderte Zierpflanze"|
dat3$Einführung == "Transportbegleiter"|
dat3$Einführung == "Saatgutbegleiter"|
dat3$Einführung == "Saatgutbegleiter, Verwilderte Nutzpflanze"|
dat3$Einführung == "Saatgutbegleiter, Verwilderte Zierpflanze"] <- "accompanying"
dat3$Introduction[dat3$Einführung == "spontan"|
dat3$Einführung == "unbekannt"|
dat3$Einführung == "spontan, Verwilderte Nutzpflanze"] <- "spontaneous"
dat3$Introduction[dat3$Einführung == "Verwilderte Nutzpflanze, Verwilderte Zierpflanze"|
dat3$Einführung == "Verwilderte Nutzpflanze"|
dat3$Einführung == "Verwilderte Zierpflanze"] <- "agric./ornam."

## convert to factor and check:
dat3$Introduction <- as.factor(dat3$Introduction)
str(dat3)

levels(dat3$Blattanatomie)
## leaf-anatomy: i will leaf this one out in the analysis
## as it is partly correlated with lifeforms

levels(dat3$Reproduktionstyp)
## "meist Samen, selten vegetativ" = mostly seed propagation
## "meist vegetativ, selten Samen" = mostly vegetative propagation
## "Samen und vegetativ" = both to same parts
## "Samen/Sporen" = only seeds
## "vegetativ" = only vegetative
## "meist Samen, selten vegetativ"; "Samen/Sporen" <-- "mostly seeds"
## "meist vegetativ, selten Samen"; "vegetativ" <-- "mostly vegetative"
## "Samen und vegetativ" <-- "both equally"
dat3$Reproduction <- rep(NA, nrow(dat3))
dat3$Reproduction[dat3$Reproduktionstyp == "meist Samen, selten vegetativ"|
dat3$Reproduktionstyp == "Samen/Sporen"] <- "mostly seeds"
dat3$Reproduction[dat3$Reproduktionstyp == "meist vegetativ, selten Samen"|
dat3$Reproduktionstyp == "vegetativ"] <- "mostly vegetativ"
dat3$Reproduction[dat3$Reproduktionstyp == "Samen und vegetativ"] <- "both equally"
dat3$Reproduction <- as.factor(dat3$Reproduction)

levels(dat3$Einbürg_grad_Max)
## this will be the dependent (outcome) variable
## i'll also try to use an outcome variable with
## only two levels:
## "Agriophyt" --> "invasive"
## "Epökophyt" & "Ephemerophyt" --> "not invasive"
dat3$Invasivness <- rep(NA, nrow(dat3))
dat3$Invasivness[dat3$Einbürg_grad_Max == "Epökophyt"|
dat3$Einbürg_grad_Max == "Ephemerophyt"] <- "not invasive"
dat3$Invasivness[dat3$Einbürg_grad_Max == "Agriophyt"] <- "invasive"
dat3$Invasivness <- as.factor(dat3$Invasivness)
dat3$Invasivness <- relevel(dat3$Invasivness , ref = "not invasive")

## another approach to this would be to define
## alien taxa that are established (in the sense
## that its population persist independently) and
## not established taxa.

dat3$Establ <- rep(NA, nrow(dat3))
dat3$Establ[dat3$Einbürg_grad_Max == "Epökophyt"|
dat3$Einbürg_grad_Max == "Agriophyt"] <- "establ."
dat3$Establ[dat3$Einbürg_grad_Max == "Ephemerophyt"] <- "not establ."
dat3$Establ <- as.factor(dat3$Establ)
dat3$Establ <- relevel(dat3$Establ, ref = "not establ.")
str(dat3)

print(ct <- ctree(Einbürg_grad_Max ~ Introduction + Lifeform + Length_Class +
Blühdauer_Mon + Reproduction, data = dat3,
control = ctree_control(mincriterion = 0.8, minsplit = 10L,
minbucket = 10L)))

plot(ct, gp = gpar(fontsize = 7))


## the first node divides annula herbs from the rest
## thus it might be more representive to recode lifeforms
## into a lifecycle variable:
dat3$Life_Cycle <- rep(NA, nrow(dat3))
dat3$Life_Cycle[grep("_ann", dat3$Lifeform)] <- "annual"
dat3$Life_Cycle[-grep("_ann", dat3$Lifeform)] <- "perennial"
dat3$Life_Cycle <- as.factor(dat3$Life_Cycle)

## re-run:
print(ct <- ctree(Einbürg_grad_Max ~ Introduction + Life_Cycle + Length_Class +
Blühdauer_Mon + Reproduction, data = dat3,
control = ctree_control(mincriterion = 0.8, minsplit = 10L,
minbucket = 10L)))

plot(ct, gp = gpar(fontsize = 7))

## check predictions:
data.frame(Species = dat3$Name,
Naturalization = dat3$Einbürg_grad_Max, predict(ct),
Correspondence = dat3$Einbürg_grad_Max == predict(ct))

table(dat3$Einbürg_grad_Max == predict(ct))

## interpreting the result:
##
## (1) there is a large probability for agriophytes (= naturalized/invasive)
## to occure in the group of perennial plants which
## exhibit both, vegetative and generative reproduction
##
## (2) ephemerophytes most likely are present if plants are annual
## and plant were introduced intentionally (ornamental, agricultural
## purpose. within annuals the group of plants that made there way
## without intentional help of men are likely to have more potential
## to get established (-> portion of agriophytes and epoekophytes (=
## naturalized/not-invasive) is larger here).
##
## also the other dependent nominal variables, invasivness & establishment
## could be investigated similarly..

To cite package ‘partykit’ in publications use:

  Torsten Hothorn and Achim Zeileis (2011). partykit: A Toolkit for Recursive
  Partytioning. R package version 0.1-0/r179.
  http://R-Forge.R-project.org/projects/partykit/
Read more »

Saturday, September 19, 2015

A Little Webscraping-Exercise...

In R it's quite easy to pull out anything from a webpage and I'll show a little exercise in doing so. Here I retrieve all blog addresses from R-bloggers by the function readLines() and some subsequent data processing.
















# get the page's html-code
web_page <- readLines("http://www.r-bloggers.com")

# extract relevant part of web page:
# missing line added on oct. 24th:
ul_tags <- grep("ul>", web_page)

pos_1 <- grep("Contributing Blogs", web_page) + 2
pos_2 <- ul_tags[which(ul_tags > pos_1)[1]] - 2

blog_list_1 <- web_page[pos_1:pos_2]

# extract 2nd element of sublists produced by stringsplit:
blog_list_2 <- unlist(lapply(strsplit(blog_list_1, "\""), "[[", 2))

# exclude elememts without propper address:
blog_list_3 <- blog_list_2[grep("http:", blog_list_2)]

# plot results:
len <- length(blog_list_3)
x <- rep(1:3, ceiling(len/3))[1:len]
y <- 1:len

par(mar = c(0, 5, 0, 5), xpd = T)
plot(x, y, ylab = "", xlab = "", type = "n",
bty = "n", axes = F)
text(x, y, blog_list_3, cex = 0.5)
Read more »