Online Trick

# purpose: crossfade 2 jpeg images
# packages: jpeg
# arguments: img1 (path.to.img1), img2 (path.to.img2),
#            outpath(defaults to current directory), outname, 
#            frames
# output: png


require(jpeg)

jpgfader <- function(img1 = NA, img2 = NA, outpath = NA, frames = NA, outname = NA){

    if(is.na(outpath)) {outpath <- path.expand("~")}

# stop if images are missing
    if(is.na(img1)|is.na(img2)) stop(cat("\nAt least one image is missing!\n"))
    if(is.na(outname)) {outname = "img.1.2"}

    if(is.na(frames)) {frames = 10}
# read 2 jpegs, assuming same size!
    pic.1 <- readJPEG(img1)
    pic.2 <- readJPEG(img2)

# warn if images dont have same size:
    if(sum(dim(pic.1) != dim(pic.2))>1) warning(cat("\nImages do not have same dimensions!"))

# create new array with 4 dimensions, the new dimension
# representing alpha:

    by = 1/(frames-1)
    alpha = seq(0, 1, by)
    n = length(alpha)

    for(j in n:1){

        pic.2.a <- array(data = c(as.vector(pic.2),
                                  rep(alpha[j], dim(pic.1)[1]*dim(pic.1)[2])),
                         dim = c(dim(pic.1)[1], dim(pic.1)[2], 4))

# assign output file name:
        pic.out <- paste(outpath, "\\", outname, ".", j ,".png", sep = "")

# and open device:
        png(pic.out, width = dim(pic.1)[2], height = dim(pic.1)[1])

# plot-parameters:
        par(mar = rep(0, 4), oma = rep(0, 4), new = F)

# print img.a to plot region:
        plot(1:2,
             xlim = c(0, dim(pic.1)[2]), ylim = c(0, dim(pic.1)[1]),
             xlab="", ylab="", type = "n",
             yaxs ="i", xaxs = "i")
        rasterImage(pic.1, 0, 0, dim(pic.1)[2], dim(pic.1)[1])

# overplotting with new alpha-pic,
# starting with full transparency, decreasing in steps, showing pic.2
# finally:
        rasterImage(pic.2.a, 0, 0, dim(pic.1)[2], dim(pic.1)[1])
        dev.off()
    }
}

# Example, with 2 images, one system.file and one altered 
# version of it:

# make black jpg and save to default system folder
Rlogo <- readJPEG(system.file("img", "Rlogo.jpg", package="jpeg"))
Rlogo[] <- 0
jpeg(path.expand("~\\Rlogo_Black.jpg"), dim(Rlogo)[2], dim(Rlogo)[1])
par(mar = rep(0, 4), oma = rep(0, 4))

# save black image:
plot(1:2,
     xlim = c(0, 1), ylim = c(0, 1),
     xlab="", ylab="", type = "n",
     yaxs ="i", xaxs = "i")
rasterImage(Rlogo, 0, 0, 1, 1)
dev.off()

# function call:
jpgfader(img1 = system.file("img", "Rlogo.jpg", package="jpeg"),
         img2 = path.expand("~/Rlogo_black.jpg"),
         outname = "img12",
         outpath = path.expand("~"),
         frames = 10)

# see the images:
browseURL(path.expand("~"))

# remove files:
# files <- dir(path.expand("~"), full.names = T)
# file.remove(c(files[grep("img12.", files)],
                path.expand("~/Rlogo_black.jpg"))) 

library(lattice)
df <- read.csv("PATH/TO/Downloads/lattice_data.csv",
               header = T, sep = ";")
mypch <- rep(21, 3)
mycol <- c(rgb(0.2, 0.8, 0.9), rgb(0.8, 0.2, 0.9), rgb(0.9, 0.2, 0.2))

stripplot(mean ~ group1 | item_parc, groups = group2, data = df,
          type = c("a", "p"), strip = strip.custom(strip.names = c(F, T)),
          ylab = "Score-Mittelwert", layout = c(2, 4),
          scales = list(y = list(cex = 0.9),
                        x = list(cex = 1.1, labels = c("männl.", "weibl."),
                        tck = c(1, 0))),
          par.strip.text = list(cex = 0.9),
          panel = function(...) {
             panel.stripplot(..., col = mycol, pch = mypch, fill = mycol)
          },
          index.cond = list(c(7,8,5,6,3,4,1,2)),
          key = list(space = "bottom", text = list(c("Alter < 35 J.",
                                                     "Alter 35-50 J.",
                                                     "Alter > 50 J.")),
                     points = list(col = mycol, cex = 0.8, pch = mypch, fill = mycol),
                     rep = F))

## comparing flowering time in indigenous and alien species by
## usage of a quantile-quantile plot (qqplot), the ks-test, by
## testing shifts in medians (wilcox test) and by testing
## difference of means (t-test, as we deal with integer and
## not a continous variable it is not the most appropiate choice)
## as well as by a chi-square test:

dat <- read.csv("E:\\R\\Data\\flowering_alien_vs_indigen.csv",
                sep = ";")

library(lattice)

histogram(~ Flowering|Status, data = dat, col = "gray60", layout = c(1, 2),
          xlab = list("Months of flowering"),
          ylab = list("Percentage of total"),
          scales = list(y = list(alternating = F)),
          strip = strip.custom(factor.levels = c("alien", "indigenous")))

qqplot(dat$Flowering[dat$Status == "indigen"],
       dat$Flowering[dat$Status == "Neophyt"])
abline(a = 0, b = 1, lty = 3)

ks.test(dat$Flowering[dat$Status == "indigen"],
        dat$Flowering[dat$Status == "Neophyt"])

wilcox.test(Flowering ~ Status, data = dat)

t.test(Flowering ~ Status, data = dat)

## Note that in the two-sample case the estimator for the
## difference in location parameters does not estimate
## the difference in medians (a common misconception) but
## rather the median of the difference between a sample
## from x and a sample from y.

## as we deal with a limited number of classes (1-12 months)
## and sample size is big enough my favourite would be a
## chi-square test:

m <- table(dat$Status, dat$Flowering)

(Xsq <- chisq.test(m))  # Prints test summary
Xsq$observed   # observed counts (same as M)
Xsq$expected   # expected counts under the null
Xsq$residuals  # Pearson residuals
Xsq$stdres     # standardized residuals

I present code for a stacked graph with common axes only for panels of the same row and with axis titles at different sides. This admittedly took me days (because i had not much of a clue how to use lattice), but eventually I did it and maybe someone can use this for his/her own purpose:

library(lattice)

y1 <- rnorm(120,100,10)
y2 <- rnorm(120,10,1)

facs <- expand.grid(Sites=rep(c("Site I","Site II",
                       "Site III"),20),Treatment = c("A","B"))

trellis.par.set(clip = list(panel = "off",strip = "off"),
   layout.widths = list(right.padding = 10,
   left.padding = 8))

PanFun<-function(...){
   panel.dotplot(...)
     if(is.element(panel.number(),1:2))
       {at<-pretty(c(min(y1),max(y1)))
        panel.axis("right",at = at,outside = T,
        labels = F,half = F)}
     if(panel.number() == 3)
       {at<-pretty(c(min(y1),max(y1)))
        panel.axis("right",at = at,outside = T,
        labels = T,half = F)}
     if(panel.number() == 4)
       {at<-pretty(c(min(y2),max(y2)))
        panel.axis("left",at = at,outside = T,
        labels = T,half = F)}
     if(is.element(panel.number(),5:6))
       {at<-pretty(c(min(y2),max(y2)))
        panel.axis("left",at = at,outside = T,
        labels = F,half = F)}}

dotplot(y1+y2 ~ Treatment|Sites,
        outer = TRUE,data = facs,
        scales = list(
             x = list(rot = 0, tck=c(1,0)),
             y = list(draw=F,relation="free",
        limits = list(range(y1),range(y1),range(y1),
                    range(y2),range(y2),range(y2)))),
        ylab = list("response one (y1)", y = 0.75, x = -4),
        ylab.right = list("response two (y2)", y = 0.25, 
                          rot = 270, x = 4),
        xlab = c("Site 1", "Site 2","Site 3"),
        strip = FALSE, panel=PanFun)

library(lme4)
library(mvtnorm)
library(multcomp)
library(lattice)

########################################################################
# setting up fake data, for details on sampling design see:            #
# cichini et al. (2011)                                                #
# http://www.springerlink.com/content/t4u305045t214882                 #
########################################################################

# 3 successional stages with 44 units and 1 succ. stage wit 10 units:
N = 44*3+10

cen <- per <- cbind(n = rep(NA, N), X = rep(NA, N))

# i simulate an effect between factor-levels cen and per
# all other fixed and random factors and interactions will have
# no effect.
# X = incidents, n = no. of observations:
for(i in 1:N){cen[i,"X"]<-sum(rbinom(5, 1, 0.35)); cen[i, "n"] <- 5}
for(i in 1:N){per[i,"X"]<-sum(rbinom(8, 1, 0.65)); per[i, "n"] <- 8}

dat <- data.frame(rbind(cen, per)) 

# probabilities:
dat$p <- dat$X/dat$n
stage <- as.factor(rep(c(rep("A",44), rep("B",44), rep("C",44), rep("D",10)),2))
plotA <- paste("A",rep(1:11,c(4,4,4,4,4,4,4,4,4,4,4)),sep="")
plotB <- paste("B",rep(1:13,c(3,4,2,4,4,4,4,4,3,4,4,2,2)),sep="")
plotC <- paste("C",rep(1:12,c(3,4,4,4,4,4,4,4,3,4,4,2)),sep="")
plotD <- paste("D",rep(1:3,c(2,4,4)),sep="")
plot <- as.factor(rep(c(plotA, plotB, plotC, plotD),2))
subplot <-  as.factor(rep(1:142,2))
pos <- as.factor(c(rep("cen", N),rep("per", N)))

df <- data.frame(cbind(dat, pos = pos, stage = stage, plot = plot, subplot = subplot))
str(df)

### visualize data:
stripplot(jitter(p,10) ~ pos | plot, data = df, groups = subplot,
strip = strip.custom(strip.names = c(F,T)),
ylab="p",
layout=c(12,4),
type = c("p", "a"),col=1,cex=0.5,
scales = list(y = list(cex=0.5),
              x = list(rot=90,cex=0.6,labels=c("Centers","Periphery"))),
              par.strip.text=list(cex=0.7),
     page = function(...) {
          ltext(x = 0.39,
                y = 0.835,
                lab = c("CENTERS vs. PERIPHERY\npanels = plots;\nrows = successional stages:\nA: pioneer, B: early, C: late, D: alpine gr.land; \nlines = paired samples (= subplots) "),
                cex = 1,adj=0)
      })

gmod <- glmer(cbind(dat$X, dat$n - dat$X) ~ stage * pos + (1|plot/subplot), family = binomial)

# comparisons -
# i want to know whether there are differences between center vs. periphery
# within each stage, and whether this effect differs across
# stages (interactions).
# as i simulated an effect only for centre vs. periphery and this
# effect was simulated
# to be the same across all stages, there should not be any sign.
# interactions:

c1 <- rbind("A: Center vs. Peri. Effect" = c(0,0,0,0,1,0,0,0),
            "B: Center vs. Peri. Effect" = c(0,0,0,0,1,1,0,0),
            "C: Center vs. Peri. Effect" = c(0,0,0,0,1,0,1,0),
            "D: Center vs. Peri. Effect" = c(0,0,0,0,1,0,0,1),
            "Center vs. Peri. Effect, A vs. B" = c(0,0,0,0,0,1,0,0),
            "Center vs. Peri. Effect, A vs. C" = c(0,0,0,0,0,0,1,0),
            "Center vs. Peri. Effect, A vs. D" = c(0,0,0,0,0,0,0,1),
            "Center vs. Peri. Effect, B vs. C" = c(0,0,0,0,0,-1,1,0),
            "Center vs. Peri. Effect, B vs. D" = c(0,0,0,0,0,-1,0,1),
            "Center vs. Peri. Effect, C vs. D" = c(0,0,0,0,0,0,-1,1))

summary(glht(gmod, c1))
###########################################################################

# EDIT:
# for investigation of contrasts you can use this functions, i.e.:
stage <- relevel(stage, "B")
model.matrix(~ stage * pos, contrasts = list(stage = "contr.treatment", pos = "contr.treatment"))

## data:
dat <- data.frame(list(structure(list(flstat = structure(c(2L, 1L, 2L, 1L),
.Label = c("allo", "auto"), class = "factor"),
reprod = structure(c(1L, 1L, 2L, 2L),
.Label = c("non-veg", "veg"), class = "factor"),
X = c(872L, 423L, 1872L, 232L)),
.Names = c("flstat", "reprod", "X"),
class = "data.frame", row.names = c(NA, -4L))))

## proportion of species with vegetative reproduction
p_allo <- dat$X[4] / (dat$X[2] + dat$X[4])
p_auto <- dat$X[3] / (dat$X[1] + dat$X[3])
p_allo
p_auto

## restructure data for glm:
dat1 <- dat[rep(1:4, dat$X), 1:2]
head(dat1)
dat1$inc <- ifelse(dat1$reprod == "non-veg", 0, 1)

## glm:
summary(gmod <- glm(inc ~ flstat, data = dat1, family = binomial))

## intercept = logit(p_allo):
print(est_p_allo <- plogis(gmod$coef[1]))

## intercept + b = logit(p_allo+p_auto):
print(est_p_auto <- plogis(gmod$coef[1] + gmod$coef[2]))

## alternatively test difference in two proportions with prop.test():
ptest_diff <- prop.test(x = c(dat$X[1], dat$X[2]),
                        n = c(dat$X[1] + dat$X[3], dat$X[2] + dat$X[4]))

## only for one proportion prop.test gives you the confidence
## intervals of p.
## (you could also extract the glm-standard errors and calculate
## the conf.int. for this purpose..):
ptest_auto <- prop.test(x = dat$X[3], n = dat$X[1] + dat$X[3])
ptest_allo <- prop.test(x = dat$X[4], n = dat$X[2] + dat$X[4])

## plot with confidence intervals from prop.test
## (see methods in ?prop.test):

## coordinates for plotting confidence interval bars:
y0_al <- ptest_allo$conf[1]
y1_al <- ptest_allo$conf[2]
y0_au <- ptest_auto$conf[1]
y1_au <- ptest_auto$conf[2]

library(grid)
library(lattice)

## panel function for suppressing tck at top and right side,
## drawing bar with confidence interval,
## plotting glm-estimates (the crosses)
mpanel = function(...) {grid.segments(x0 = c(0.2725, 1 - 0.2725),
                                      x1 = c(0.2725, 1 - 0.2725),
                                      y0 = c(y0_al, y0_au),
                                      y1 = c(y1_al, y1_au))
                        panel.points(x = c(0.8, 2.2),
                                    y = c(est_p_allo, est_p_auto), pch = 4)
                        panel.abline(h = c(p_allo, p_auto), lty = 15,
                                     col = "grey70")
                        panel.text(x = 1.5, y = 0.9, cex = 1.2,
                                   "Species With\nVegetative Reproduction");
                        panel.xyplot(...)}

xyplot(c(p_allo, p_auto) ~ as.factor(c("Alien", "Indigenous")), type = "b",
       ylab = "Prop. +/- CIs\nX = GLM-Estimates",
       xlab = "", ylim = c(0, 1),
       panel = mpanel, pch = 16,
       scales = list(alternating = 1, tck = c(1, 0)))

library(vegan)

### data on 35 species of Oribatid mites
### data description at: http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/mite.html
data(mite)

### ...and environmental variables
data(mite.env)

### the factor which you wish to use for labeling the displayed sites
### in the ordination diagramm:
fac<-mite.env$Topo
sol<-metaMDS(mite)

windows(9,5)
par(mfrow=c(1,2))

### graph with "spider":
plot(sol,type="n")
points(sol, display = "sites", cex=0.65, select=which(fac=="Hummock"),
       pch = 21,col="black", bg="black")
points(sol, display = "sites", cex=0.65, select=which(fac=="Blanket"),
       pch = 21,col="black", bg="white")
ordihull(sol,group=fac,show.groups="Hummock")
ordihull(sol,group=fac,show.groups="Blanket",lty=3)
orditorp(sol, dis = "sp", pcex=0,air=0.85,col="grey35",cex=0.8 ,font=3)

### graph with "hull":
plot(sol,type="n")
points(sol, display = "sites", cex=0.65, select=which(fac=="Hummock"),
       pch = 21,col="black", bg="black")
points(sol, display = "sites", cex=0.65, select=which(fac=="Blanket"),
       pch = 21,col="black", bg="white")
ordispider(sol,group=fac,show.groups="Hummock")
ordispider(sol,group=fac,show.groups="Blanket",lty=3)
orditorp(sol, dis = "sp", pcex=0,air=0.85,col="grey35",cex=0.8 ,font=3)

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..