.Rhistory

knitr::opts_chunk$set(echo = FALSE)
knitr::include_graphics("progit2.png")
knitr::include_graphics("progit2.png")
knitr::include_graphics("progit2.png")
knitr::include_graphics("progit2.png")
ls()
ls
ls()
x <- 3
ls()
knitr::include_graphics("signup.png")
knitr::include_graphics("yourprofile.png")
knitr::include_graphics("your-profile.png")
knitr::include_graphics("my-profile.png")
knitr::include_graphics("new-repo.png")
knitr::include_graphics("new-repo.png")
knitr::include_graphics("new-repo.png")
knitr::include_graphics("fever.png")
knitr::include_graphics("fever.png")
knitr::include_graphics("fever.PNG")
knitr::include_graphics("fever.jpg")
knitr::include_graphics("fever.jpg")
knitr::include_graphics("fever.jpg")
knitr::include_graphics("Octocat.png")
?.subset
install.packages(c("investr", "scam"))
library(investr)
library(mgcv)
library(scam)
### GAM/SCAM/GP/Kriging
#' Inverse Prediction Function
#'
#' Reverse predicts x for a given y.
#'
#' @param object: A fitted model object.
#' objects).
#'
#' @param y: The response value of interest.
#'
#' @param x.name: Character string giving the name of the (single) independent
#' variable.
#'
#' @param ... Additional arguments to be passed onto \code{uniroot}.
#'
#' @export
inv.predict <- function(object, y, x.name, ...) {
UseMethod("inv.predict")
}
#' @rdname inv.predict
#' @export
inv.predict.default <- function(object, y, x.name, interval = FALSE,
lower, upper, level = 0.95,...) {
x0.est <- uniroot(lower = lower, upper = upper, extendInt = "yes", ..., f = function(x) {
predFit(object, newdata = setNames(data.frame(x), x.name)) - y
})$root
res <- if (interval) {
lwr <- uniroot(lower = lower, upper = x0.est, extendInt = "yes", ..., f = function(x) {
predFit(object, newdata = setNames(data.frame(x), x.name),
interval = "confidence")[, "upr"] - y
})$root
upr <- uniroot(lower = x0.est, upper = upper, extendInt = "yes", ..., f = function(x) {
predFit(object, newdata = setNames(data.frame(x), x.name),
interval = "confidence")[, "lwr"] - y
})$root
lwr <- min(c(lwr, upr))
upr <- max(c(lwr, upr))
c("estimate" = x0.est, "lower" = lwr, "upper" = upr)
} else {
x0.est
}
res
}
#' GAM prediction
#' @rdname predFit
#' @export
predFit.gam <- function(object, newdata, type = c("link", "response"),
interval = c("none", "confidence", "prediction"),
level = 0.95, ...) {
type <- match.arg(type)
interval <- match.arg(interval)
res <- if (interval == "none") {
predict.gam(object, newdata = newdata, type = type, ...)
} else if (interval == "confidence") {
pred <- predict.gam(object, newdata = newdata, se.fit = TRUE, type = "link",
...)
out <- cbind("fit" = pred$fit,
"lwr" = pred$fit - pred$se.fit * stats::qnorm((level + 1) / 2),
"upr" = pred$fit + pred$se.fit * stats::qnorm((level + 1) / 2))
if (type == "response") {
out <- apply(out, MARGIN = 2, FUN = function(x) {
stats::family(object)$linkinv(x)
})
}
out
} else {
stop("Prediction intervals are currently not supported for GAMs.")
}
res
}
#' SCAM prediction
#' @rdname predFit
#' @export
predFit.scam <- function(object, newdata, type = c("link", "response"),
interval = c("none", "confidence", "prediction"),
level = 0.95, ...) {
type <- match.arg(type)
interval <- match.arg(interval)
res <- if (interval == "none") {
predict.scam(object, newdata = newdata, type = type, ...)
} else if (interval == "confidence") {
pred <- predict.scam(object, newdata = newdata, se.fit = TRUE, type = "link",
...)
out <- cbind("fit" = pred$fit,
"lwr" = pred$fit - pred$se.fit * stats::qnorm((level + 1) / 2),
"upr" = pred$fit + pred$se.fit * stats::qnorm((level + 1) / 2))
if (type == "response") {
out <- apply(out, MARGIN = 2, FUN = function(x) {
stats::family(object)$linkinv(x)
})
}
out
} else {
stop("Prediction intervals are currently not supported for SCAMs.")
}
res
}
# ptm <- proc.time()
N <- 50 ### How many individuals are in each permutation
Hstar <- 2 ### Assume we know how many total haplotypes are possible
perm <- 10000 ### How many permutations are we considering
specs <- 1:N
### Set up a container to hold the identity of each individual from each permutation
pop <- array(dim = c(perm, N))
### Create an ID for each haplotype
haps <- as.character(1:Hstar)
### Assign probabilities of occurrence to each haplotype, ensure they sum to 1
### This is where we assume we "know" the distribution of haplotypes
### Here, I have assumed they all occur with equal frequency, but you can change this to assume some dominant ones and some
###  rare ones, whatever you want
# probs <- rep(1/Hstar, Hstar)
probs <- c(0.5, 0.5)
# probs <- runif(H) # random values from U(0, 1) distribution
# probs <- sort(probs/sum(probs), decreasing = TRUE)
### Generate permutations, we assume each permutation has N individuals, and we sample those individuals' haplotypes
### from our probabilities
for(i in 1:perm){
pop[i,] <- sample(haps, size = N, replace = TRUE, prob = probs)
}
### Make a matrix to hold the 1:N individuals from each permutation
HAC.mat <- array(dim = c(perm, N))
for(j in specs){
for(i in 1:perm){
ind.index <- sample(specs, size = j, replace = FALSE) ## which individuals will we sample
hap.plot <- pop[i, ind.index]	## pull those individuals from a permutation
HAC.mat[i,j] <- length(unique(hap.plot))  ## how many haplotypes did we get for a given sampling intensity (j) from each         ### permutation (i)
}
}
### Calculate the mean and CI for number of haplotypes at each sampling intensity (j)
means <- apply(HAC.mat, MARGIN = 2, mean)
lower <- apply(HAC.mat, MARGIN = 2, function(x) quantile(x, 0.025))
upper <- apply(HAC.mat, MARGIN = 2, function(x) quantile(x, 0.975))
### Plot the curve and frequency barplot
par(mfrow = c(1, 2))
plot(specs, means, type = "n", xlab= "Specimens sampled", ylab = "Unique haplotypes", ylim = c(1, H))
polygon(x = c(specs, rev(specs)), y = c(lower, rev(upper)), col = "gray")
lines(specs, means, lwd = 2)
# HAC.bar <- barplot(N*probs, xlab = "Unique haplotypes", ylab = "Specimens sampled", names.arg = 1:H)
d <- data.frame(specs, means)
plot(specs, means, type = "n", xlab= "Specimens sampled", ylab = "Unique haplotypes")
polygon(x = c(specs, rev(specs)), y = c(lower, rev(upper)), col = "gray")
lines(specs, means, lwd = 2)
d <- data.frame(specs, means)
HAC.mod <- lm(log(means) ~ log(specs), data = d)
HAC.nls <- nls(means ~ a*specs^b, start=list(a = exp(coef(HAC.mod)[1]), b = coef(HAC.mod)[2]), data = d)
AIC(HAC.nls)
### Plot curve
x <- specs
curve(coef(HAC.nls)[1]*x^coef(HAC.nls)[2], add = TRUE, from = 1, to = N, col = "red", lwd = 2)
### Calculate N*
HAC.invest <- invest(HAC.nls, y0 = Hstar, interval = "inversion", upper = 1000000)
HAC.invest
# GAMs
HAC.cr <- gam(means ~ s(specs, bs = "cr", k = 10), optimizer = c("outer", "bfgs"), data = d)
HAC.cr$aic
HAC.cr <- inv.predict(HAC.cr, y = Hstar, x.name = "specs", interval = TRUE, lower = 1, upper = 1000000)
HAC.cr
plot(HAC.cr)
Hstar
library(pdp)
partial(HAC.cr, pred.var = "specs", plot = TRUE)
plot(d)
lines(xx, yy, type = "l")
xx <- seq(from = min(d$specs), to = max(d$specs))
yy <- predict(HAC.dr, newdata = data.frame(specs = xx))
yy <- predict(HAC.cr, newdata = data.frame(specs = xx))
plot(d)
lines(xx, yy, type = "l")
plot(d, ylim = c(1, 2))
lines(xx, yy, type = "l")
plot(d, ylim = c(1, 2.5))
lines(xx, yy, type = "l")
plot(d, ylim = c(0.8, 2.2))
lines(xx, yy, type = "l")
inv.predict(HAC.cr, y = 1.6, x.name = "specs", interval = TRUE, lower = 1, upper = 1000000)
inv.predict(HAC.cr, y = 1.6, x.name = "specs", interval = TRUE)
inv.predict(HAC.cr, y = 1.6, x.name = "specs", interval = TRUE, lower = 0, upper = 10)
xx <- seq(from = min(d$specs), to = max(d$specs) + 100000, length = 500)
yy <- predict(HAC.cr, newdata = data.frame(specs = xx))
plot(d, ylim = c(0.8, 2.2))
lines(xx, yy, type = "l")
partial(HAC.nls, pred.var = "specs", plot = T)
xx <- seq(from = min(d$specs), to = max(d$specs), length = 500)
yy <- predict(HAC.cr, newdata = data.frame(specs = xx))
plot(d, ylim = c(0.8, 2.2))
lines(xx, yy, type = "l")
Hstar
abline(h = Hstar, lty = 2, col = "red")
HAC.gps <- gam(means ~ s(specs, bs = "gp", k = 10, m = 1), optimizer = c("outer", "bfgs"), data = d) # spherical
plot(HAC.gps)
HAC.mpi <- scam(means ~ s(specs, bs = "mpi", k = 10), data = d) # mpi
plot(HAC.mpi)
HAC.mpi <- inv.predict(HAC.mpi, y = Hstar, x.name = "specs", interval = TRUE, lower = 1, upper = 1000000)
xx <- seq(from = min(d$specs), to = max(d$specs), length = 500)
yy <- predict(HAC.mpi, newdata = data.frame(specs = xx))
plot(d, ylim = c(0.8, 2.2))
lines(xx, yy, type = "l")
abline(h = Hstar, lty = 2, col = "red")  # looks like multiple solutions
HAC.mpi <- inv.predict(HAC.mpi, y = Hstar, x.name = "specs", interval = TRUE, lower = 1, upper = 10)
devtools::install_github("hadley/emo")
add
source('C:/Users/greenweb/Desktop/testrepo/add.R', echo=TRUE)
add(1, 3.5)
library(ggplot2)
ggplot(cars, aes_string(x = "dist", y = "speed") + geom_point()
)
ggplot(cars, aes_string(x = "dist", y = "speed")) + geom_point()
file.exists("~/.ssh/id_rsa.pub")