### R code from vignette source 'genphenManual.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: genphenManual.Rnw:363-392
###################################################
require(genphen)
require(ggplot2)
require(knitr)
require(ggrepel)
require(reshape2)
require(ape)
require(xtable)
options(xtable.floating = FALSE)

# genotype inputs:
data(genotype.saap)
# phenotype inputs:
data(phenotype.saap)

# run genphen
c.out <- runGenphen(genotype = genotype.saap[, 80:85],
                    phenotype = phenotype.saap,
                    phenotype.type = "continuous",
                    mcmc.chains = 2,
                    mcmc.iterations = 2000,
                    mcmc.warmup = 500,
                    mcmc.cores = 4,
                    hdi.level = 0.95,
                    stat.learn.method = "rf",
                    cv.iterations = 300)

# scores
c.convergence <- c.out$convergence
c.score <- c.out$scores


###################################################
### code chunk number 2: genphenManual.Rnw:404-420
###################################################
df <- data.frame(genotype.saap[, 80:85], stringsAsFactors = FALSE)
df$phenotype <- phenotype.saap
df <- melt(data = df, id.vars = "phenotype")
colnames(df) <- c("phenotype", "site", "genotype")
df$site <- gsub(pattern = "X", replacement = '', x = df$site)

ggplot(data = df)+
  facet_wrap(facets = ~site, ncol = 6, scales = "free_x")+
  geom_violin(aes(x = genotype, y = phenotype))+
  geom_point(aes(x = genotype, y = phenotype, col = genotype), size = 1,
             shape = 21, position = position_jitterdodge(jitter.width = 0.2,
                                                         jitter.height = 0.1,
                                                         dodge.width = 0.5))+
  scale_color_discrete(name = "genotype")+
  theme_bw(base_size = 16)+
  theme(legend.position="none")


###################################################
### code chunk number 3: genphenManual.Rnw:435-447
###################################################
c.score$label <- paste(c.score$site, c.score$mutation,
                       round(x = c.score$bc, digits=2), sep=':')
ggplot(data = c.score)+
  geom_errorbar(aes(x = ca, ymin = cohens.d.L, ymax = cohens.d.H))+
  geom_point(aes(x = ca, y = cohens.d, fill = kappa), shape = 21, size = 3)+
  geom_text_repel(aes(x = ca, y = cohens.d, label = label), size = 4)+
  theme_bw(base_size = 16)+
  xlab(label = "CA")+
  ylab(label = "Cohen's d (with 95% HDI)")+
  scale_x_continuous(limits = c(0, 1))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  scale_fill_gradientn(colours = terrain.colors(n = 10), limits = c(0, 1))


###################################################
### code chunk number 4: genphenManual.Rnw:453-458
###################################################
print(xtable(c.out$debug.scores[, c("site", "mutation", "cohens.d",
                                    "cohens.d.hdi", "bc", "ca",
                                    "ca.hdi", "kappa", "kappa.hdi"), ],
             align = rep(x = "c", times = 10, digits = 2)),
      include.rownames = FALSE, size = "small")


###################################################
### code chunk number 5: genphenManual.Rnw:471-473
###################################################
print(xtable(c.convergence, align = rep(x = "c", times = 10, digits = 2)),
      include.rownames = FALSE, size = "small")


###################################################
### code chunk number 6: genphenManual.Rnw:484-499
###################################################
# compute the phylogenetic bias
bias <- runPhyloBiasCheck(genotype = genotype.saap, input.kinship.matrix = NULL)

# extract kinship matrix
kinship.matrix <- bias$kinship.matrix

# extract the bias associated with mutations of the sites which were included
# in the association analysis
mutation.bias <- bias$bias.mutations
mutation.bias <- mutation.bias[mutation.bias$site %in% 80:85, ]
mutation.bias <- mutation.bias[mutation.bias$mutation %in% c.score$mutation, ]

# show the bias table
print(xtable(mutation.bias, align = rep(x = "c", times = 3+1, digits = 2)),
      include.rownames = FALSE, size = "small")


###################################################
### code chunk number 7: genphenManual.Rnw:511-528
###################################################
color.a <- character(length = nrow(genotype.saap))
color.a[1:length(color.a)] <- "gray"
color.a[which(genotype.saap[, 82] == "h")] <- "orange"
color.a[which(genotype.saap[, 82] == "q")] <- "blue"

color.b <- character(length = nrow(genotype.saap))
color.b[1:length(color.b)] <- "gray"
color.b[which(genotype.saap[, 84] == "a")] <- "orange"
color.b[which(genotype.saap[, 84] == "d")] <- "blue"

c.hclust <- hclust(as.dist(kinship.matrix), method = "average")

par(mfrow = c(1, 2), mar = c(0,0,1,0) + 0.1)
plot(as.phylo(c.hclust), tip.color = color.a, cex = 0.6, type = "fan",
     main = "B = 0.15")
plot(as.phylo(c.hclust), tip.color = color.b, cex = 0.6, type = "fan",
     main = "B = 0.43")


###################################################
### code chunk number 8: genphenManual.Rnw:541-561
###################################################
# genotype inputs:
data(genotype.saap)
# phenotype inputs:
data(dichotomous.phenotype.saap)

# run genphen
d.out <- runGenphen(genotype = genotype.saap[, 68:79],
                    phenotype = dichotomous.phenotype.saap,
                    phenotype.type = "dichotomous",
                    mcmc.chains = 2,
                    mcmc.iterations = 2000,
                    mcmc.warmup = 500,
                    mcmc.cores = 4,
                    hdi.level = 0.95,
                    stat.learn.method = "rf",
                    cv.iterations = 300)

# scores
d.convergence <- d.out$convergence
d.score <- d.out$scores


###################################################
### code chunk number 9: genphenManual.Rnw:573-588
###################################################
df <- data.frame(genotype.saap[, 68:79], stringsAsFactors = FALSE)
df$phenotype <- dichotomous.phenotype.saap
df <- melt(data = df, id.vars = "phenotype")
colnames(df) <- c("phenotype", "site", "genotype")
df$site <- gsub(pattern = "X", replacement = '', x = df$site)

ggplot(data = df)+
  facet_wrap(facets = ~site, ncol = 6, scales = "free_x")+
  geom_point(aes(x = genotype, y = phenotype, col = genotype), size = 1,
             shape = 21, position = position_jitterdodge(jitter.width = 0.2,
                                                         jitter.height = 0.1,
                                                         dodge.width = 0.5))+
  scale_color_discrete(name = "genotype")+
  theme_bw(base_size = 16)+
  theme(legend.position="none")


###################################################
### code chunk number 10: genphenManual.Rnw:603-616
###################################################
d.score$label <- paste(d.score$site, d.score$mutation,
                       round(x = d.score$bc, digits = 2),
                       sep = ':')
ggplot(data = d.score)+
  geom_errorbar(aes(x = ca, ymin = absolute.d.L, ymax = absolute.d.H))+
  geom_point(aes(x = ca, y = absolute.d, fill = kappa), shape = 21, size = 3)+
  geom_text_repel(aes(x = ca, y = absolute.d, label = label), size = 4)+
  theme_bw(base_size = 16)+
  xlab(label = "CA")+
  ylab(label = "Absolute d (with 95% HDI)")+
  scale_x_continuous(limits = c(0, 1))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  scale_fill_gradientn(colours = terrain.colors(n = 10), limits = c(0, 1))


###################################################
### code chunk number 11: genphenManual.Rnw:625-627
###################################################
print(xtable(d.convergence, align = rep(x = "c", times = 8, digits = 2)),
      include.rownames = FALSE, size = "small")


###################################################
### code chunk number 12: genphenManual.Rnw:657-697
###################################################
set.seed(seed = 551155)
g1 <- replicate(n=10^4, expr = as.character(rbinom(n=30, size = 1, prob = 0.5)))
g2 <- replicate(n=10^4, expr = as.character(rbinom(n=30, size = 1, prob = 0.5)))
gen <- rbind(g1, g2)
phen <- c(rnorm(n = 30, mean = 3, sd = 3), rnorm(n = 30, mean = 5, sd = 3))

# run diagnostics
diag.pre <- runDiagnostics(genotype = gen,
                           phenotype = phen,
                           phenotype.type = "continuous",
                           rf.importance.trees = 50000,
                           with.anchor.points = FALSE,
                           mcmc.chains = 2,
                           mcmc.iterations = 1500,
                           mcmc.warmup = 500,
                           mcmc.cores = 2,
                           hdi.level = 0.95,
                           anchor.points = c(1:10))

# run diagnostics
diag.post <- runDiagnostics(genotype = gen,
                           phenotype = phen,
                           phenotype.type = "continuous",
                           rf.importance.trees = 50000,
                           with.anchor.points = TRUE,
                           mcmc.chains = 2,
                           mcmc.iterations = 1500,
                           mcmc.warmup = 500,
                           mcmc.cores = 2,
                           hdi.level = 0.99,
                           anchor.points = c(1:5, 101:105, 301:305,
                                             501:505, 1001:1005,
                                             2001:2005, 5001:5005,
                                             9001:9005))

pre.scores <- diag.pre$importance.scores
post.scores <- diag.post$scores
post.scores$significant <- ifelse(test = post.scores$cohens.d.L <= 0
                                  & post.scores$cohens.d.H >= 0,
                                  yes = FALSE, no = TRUE)


###################################################
### code chunk number 13: genphenManual.Rnw:707-714
###################################################
ggplot(data = pre.scores)+
  geom_step(aes(x = importance.rank, y = importance))+
  xlab("Rank")+
  ylab("Importance")+
  ggtitle("Importance distribution")+
  theme_bw(base_size = 16)+
  scale_fill_gradientn(colours = terrain.colors(n = 10), limits = c(0, 1))


###################################################
### code chunk number 14: genphenManual.Rnw:732-744
###################################################
ggplot(data = post.scores)+
  geom_errorbar(aes(x = anchor.point, ymin = cohens.d.L, ymax = cohens.d.H),
                col = "gray", width = 0.1)+
  geom_point(aes(x = anchor.point, y = cohens.d, col = significant), size = 2)+
  geom_hline(yintercept = 0, linetype = "dashed", col = "gray")+
  theme_bw(base_size = 16)+
  xlab(label = "Rank")+
  ylab(label = "Cohen's d (with 99% HDI)")+
  scale_x_continuous(trans = "log10", breaks = c(1, 10, 100, 1000, 10000),
                     labels = c(1, 10, 100, 1000, 10000))+
  annotation_logticks(sides = "b")+
  theme(legend.position = "top")