Author Archives: tutor

Würzburg rf code

library(randomForest)
library(party)
library(nnet)
library(caret)
library(latticeExtra)
library(gridExtra)
library(relaimpo)
########################## regression example #############################
## global operations
data(readingSkills)

set.seed(123223)
sample.ind <- sample(1:nrow(readingSkills), 150, replace = FALSE)
train.set <- readingSkills[sample.ind, ]
pred.set <- readingSkills[-sample.ind, ]

## party example
## ————————————————————————
set.seed(42)
train.cf <- cforest(score ~ ., data = train.set,
control = cforest_unbiased(mtry = 2, ntree = 50))

fit.cf <- predict(train.cf, newdata = pred.set)

scat.cf <- xyplot(fit.cf ~ pred.set$score, pch = 16,
asp = 1, col = “black”,
xlim = c(23, 58), ylim = c(23, 58),
xlab = “observed score”, ylab = “predicted score”,
panel = function(x, y, …) {
panel.xyplot(x, y, …)
lm1 <- lm(y ~ x)
lm1sum <- summary(lm1)
r2 <- lm1sum$adj.r.squared
panel.text(labels =
bquote(italic(R)^2 ==
.(format(r2,
digits = 3))),
x = 27.5, y = 55)
panel.smoother(x, y, method = “lm”,
col = “black”,
col.se = “black”,
alpha.se = 0.3)
panel.abline(a = 0, b = 1, col = “red”, lty = 2)
})
scat.cf

set.seed(42)
varimp(train.cf)

set.seed(42)
varimp(train.cf, conditional = TRUE)

set.seed(42)
dotplot(sort(varimp(train.cf, conditional = TRUE)))
## ————————————————————————
## randomForest example
## ————————————————————————
set.seed(42)
train.rf <- randomForest(score ~ ., data = train.set,
ntree = 50,
mtry = 2,
importance = TRUE,
keep.forest = TRUE,
do.trace = 100)

fit.rf <- predict(train.rf, pred.set)

scat.rf <- xyplot(fit.rf ~ pred.set$score, pch = 16,
asp = 1, col = “black”,
xlim = c(23, 58), ylim = c(23, 58),
xlab = “observed score”, ylab = “predicted score”,
panel = function(x, y, …) {
panel.xyplot(x, y, …)
lm1 <- lm(y ~ x)
lm1sum <- summary(lm1)
r2 <- lm1sum$adj.r.squared
panel.text(labels =
bquote(italic(R)^2 ==
.(format(r2,
digits = 3))),
x = 27.5, y = 55)
panel.smoother(x, y, method = “lm”,
col = “black”,
col.se = “black”,
alpha.se = 0.3)
panel.abline(a = 0, b = 1, col = “red”, lty = 2)
})
scat.rf

importance(train.rf, scale = FALSE)

varImpPlot(train.rf, scale = FALSE)
## ————————————————————————
## neural network (caret – nnet) example
## ————————————————————————
set.seed(42)
train.nn <- train(score ~ ., data = train.set, method = “nnet”,
preProcess = “range”, linout = TRUE, maxit = 100)

fit.nn <- predict(train.nn, pred.set, type = “raw”)

scat.nn <- xyplot(fit.nn ~ pred.set$score, pch = 16,
asp = 1, col = “black”,
xlim = c(23, 58), ylim = c(23, 58),
xlab = “observed score”, ylab = “predicted score”,
panel = function(x, y, …) {
panel.xyplot(x, y, …)
lm1 <- lm(y ~ x)
lm1sum <- summary(lm1)
r2 <- lm1sum$adj.r.squared
panel.text(labels =
bquote(italic(R)^2 ==
.(format(r2,
digits = 3))),
x = 27.5, y = 55)
panel.smoother(x, y, method = “lm”,
col = “black”,
col.se = “black”,
alpha.se = 0.3)
panel.abline(a = 0, b = 1, col = “red”, lty = 2)
})
scat.nn

varImp(train.nn)

dotPlot(varImp((train.nn)))
## ————————————————————————
## linear model example
## ————————————————————————
train.lm <- lm(score ~ ., data = train.set)

fit.lm <- predict(train.lm, newdata = pred.set)

scat.lm <- xyplot(fit.lm ~ pred.set$score, pch = 16,
asp = 1, col = “black”,
xlim = c(23, 58), ylim = c(23, 58),
xlab = “observed score”, ylab = “predicted score”,
panel = function(x, y, …) {
panel.xyplot(x, y, …)
lm1 <- lm(y ~ x)
lm1sum <- summary(lm1)
r2 <- lm1sum$adj.r.squared
panel.text(labels =
bquote(italic(R)^2 ==
.(format(r2,
digits = 3))),
x = 27.5, y = 55)
panel.smoother(x, y, method = “lm”,
col = “black”,
col.se = “black”,
alpha.se = 0.3)
panel.abline(a = 0, b = 1, col = “red”, lty = 2)
})
scat.lm

calc.relimp(train.lm,type = c(“lmg”, “last”, “first”, “pratt”),
rela = TRUE)

# Bootstrap Measures of Relative Importance (1000 samples)
boot <- boot.relimp(train.lm, b = 1000,
type = c(“lmg”, “last”, “first”, “pratt”),
rank = TRUE, diff = TRUE, rela = TRUE)

booteval.relimp(boot)

plot(booteval.relimp(boot, sort = TRUE))
## ————————————————————————
###########################################################################
########################## classification example #########################

## global

data(iris)

set.seed(42)
sample.ind <- sample(1:nrow(iris), 100, replace = FALSE)
train.set <- iris[sample.ind, ]
pred.set <- iris[-sample.ind, ]
## party example rf
## ————————————————————————
set.seed(42)
train.cf <- cforest(Species ~ ., data = train.set,
control = cforest_unbiased(mtry = 2, ntree = 50))

fit.cf <- predict(train.cf, newdata = pred.set)

ftable(fit.cf ~ pred.set$Species)

set.seed(42)
varimp(train.cf)

set.seed(42)
varimp(train.cf, conditional = TRUE)

set.seed(42)
dotplot(sort(varimp(train.cf, conditional = TRUE)))
## ————————————————————————

## party example ctree
# ————————————————————————-
train.ct <- ctree(Species ~ ., data = train.set)

plot(train.ct)

fit.ct <- predict(train.ct, newdata = pred.set)

result <- ftable(fit.ct ~ pred.set$Species)

grid.newpage()
grid.table(as.data.frame(result))

out <- data.frame(expand.grid(attr(result, “row.vars”)),
unclass(result))
rownames(out) <- out[, 1]
out <- out[, -1]
colnames(out) <- rownames(out)

grid.newpage()
gt <- tableGrob(out,
gpar.rowtext = gpar(col = “black”, cex = 1,
fontface = “bold”))

h <- grobHeight(gt)
w <- grobWidth(gt)

col.title <- textGrob(“Predicted”, y = unit(0.5, “npc”) + 0.55 * h,
vjust = 0, gp = gpar(fontface = “italic”))
row.title <- textGrob(“Observed”, x = unit(0.5, “npc”) – 0.52 * w,
y = unit(0.5, “npc”) – 0.55 * h,
vjust = 0, hjust = 0, rot = 90,
gp = gpar(fontface = “italic”))

gtr <- gTree(children = gList(gt, col.title, row.title))

grid.draw(gtr)

Würzburg map code

library(sp)
library(RColorBrewer)
library(raster)
library(rgdal)
library(OpenStreetMap)
library(latticeExtra)
library(plyr)
library(colorspace)
library(gridExtra)
library(scales)

# library(munsell)
# library(dichromat)
# library(labeling)

library(Rtography)
# library(devtools)
# install_github(‘Rtography’, ‘tim-salabim’)
# library(Rtography)

# rasterOptions(tmpdir = “/media/HDEXTENSION/”)
### set global options
map.zoom <- 12
map.type = ‘bing’ #’opencyclemap’ #
#launchMapHelper()
### read plot data
#plots <- read.csv(‘/media/windows/tappelhans/uni/marburg/kili/plots/kili_plots_64.csv’)
plots <- read.csv(‘C:/tappelhans/uni/marburg/kili/plots/kili_plots_64.csv’)

#plots <- subset(plots, Categories == “Helichrysum”)
#plots <- droplevels(plots)

coordinates(plots) <- ~ Easting + Northing
proj4string(plots) <- CRS(“+proj=utm +ellps=clrk80 +zone=37 +units=m +south”)

plots.ll <- spTransform(plots, CRS(“+proj=longlat”))
### define extent of map based on points and download
#bbPoints <- bbox(plots.ll)
xmin <- bbox(plots.ll)[1, 1] – 0.2
ymin <- bbox(plots.ll)[2, 1] – 0.2
xmax <- bbox(plots.ll)[1, 2] + 0.2
ymax <- bbox(plots.ll)[2, 2] + 0.2

ul <- c(ymax, xmin)
lr <- c(ymin, xmax)

omap.merc <- openmap(ul, lr, type = map.type, zoom = map.zoom)

# plot(omap.merc)

omap.utm <- openproj(omap.merc,
“+proj=utm +ellps=clrk80 +zone=37 +units=m +south”)
### use downloaded map for sp raster layout definition
omap_cols <- matrix(omap.utm$tiles[[1]]$colorData,
nrow = omap.utm$tiles[[1]]$xres,
ncol = omap.utm$tiles[[1]]$yres)
attr(omap_cols, “class”) <- c(“ggmap”, “raster”)
attr(omap_cols, “bb”) <- data.frame(ll.lat = omap.utm$tiles[[1]]$bbox$p2[2],
ll.lon = omap.utm$tiles[[1]]$bbox$p1[1],
ur.lat = omap.utm$tiles[[1]]$bbox$p1[2],
ur.lon = omap.utm$tiles[[1]]$bbox$p2[1])

bbMap <- attr(omap_cols, ‘bb’)
latCenter <- with(bbMap, ll.lat + ur.lat)/2
lonCenter <- with(bbMap, ll.lon + ur.lon)/2
height <- with(bbMap, ur.lat – ll.lat)
width <- with(bbMap, ur.lon – ll.lon)

## Use sp.layout of spplot: a list with the function and its
## arguments
sp.raster <- list(‘grid.raster’, omap_cols,
x=lonCenter, y=latCenter,
width=width, height=height,
default.units=’native’)
### define colours for study plot locations
n <- length(unique(plots$Categories))

clrs.hcl <- function(n) {
hcl(h = seq(0, 270, length.out = n),
c = 60, l = 50, fixup = TRUE)
}

 

### create spplot object and Rtographize
pts <- spplot(plots["Categories"], col.regions = clrs.hcl(n),
sp.layout = list(sp.raster))

pts <- Rtographize(pts)
### add map objects using viewports
#png(“KiLi_plots_hel.png”, width = 1024 * 3, height = 768 * 3, res = 300)
grid.newpage()
print(pts)

downViewport(trellis.vpname(name = “figure”))
vp1 <- viewport(x = unit(0.99, “npc”), y = unit(0.01, “npc”),
width = 0.25, height = 0.1, just = c(“right”, “bottom”),
name = “inset_scalebar”, gp = gpar(cex = 0.5))

pushViewport(vp1)
#grid.rect(gp = gpar(fill = “transparent”, col = “transparent”))
gridScaleBar(pts, addParams = list(noBins = 4, scale.fact = 1.5,
vpwidth = as.numeric(vp1$width)))
upViewport(1)

vp2 <- viewport(x = unit(0.01, “npc”), y = unit(0.01, “npc”),
width = 0.05, height = 0.15, just = c(“left”, “bottom”),
name = “inset_northarrow”, gp = gpar(cex = 0.8))

pushViewport(vp2)
gridNorthArrow()
upViewport(1)

vp3 <- viewport(x = unit(0.98, “npc”), y = unit(0.98, “npc”),
width = 0.2, height = 0.025 * n,
just = c(“right”, “top”),
name = “inset_legend”)

pushViewport(vp3)
gridMapLegend(labs = plots$Categories, clrs = clrs.hcl(n),
type = “circles”)
upViewport(0)

#dev.off()