mikelove

npg <- 20
mu <- c(200,400,400,800,800,1600)
cond <- rep(rep(c("A","B"),each=npg),3)
geno <- rep(c("X","Y","Z"),each=2*npg)
table(cond, geno)
counts <- rnbinom(6*npg, mu=rep(mu,each=npg), size=1/.01)
d <- data.frame(counts, cond, geno)
library(ggplot2)
plotit <- function(d) {
  ggplot(d, aes(x=cond, y=counts, group=geno)) + 

mikelove

my_test <- function() {
  tests <- list.files("tests/testthat","test_.*")
  for (i in seq_along(tests)) {
    message(paste0(i," / ",length(tests),": ",tests[i]))
    test(filter=sub("test_(.*).R","\\1",tests[i]))
  }
}

mikelove

# run linear model for each unique level of 'cyl' and return R^2 
library(purrr)
mtcars %>% 
  split(.$cyl) %>%  
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%  
  map_dbl("r.squared")

# in base R this might look like
mtcars$cyl <- factor(mtcars$cyl)

mikelove

library(purrr)
library(dplyr)

# some functions

# just a convenience function, gives back random assignments
# conceptually like: sample(labels, size=n, replace=TRUE, prob=prob)
random_group <- function(n, probs) {
  probs <- probs / sum(probs)
  g <- findInterval(seq(0, 1, length = n), c(0, cumsum(probs)),

mikelove

# the task here is just to compare the mean and median as estimators
# for the location of the distribution. I want to try 100 reps for different
# sample size, and also for different types of distribution (normal vs t)

nrep <- 100
d <- expand.grid(n=c(3,5,10,20), type=c("normal","t"))
d <- d[rep(seq_len(nrow(d)),each=nrep),]

res <- lapply(seq_len(nrow(d)), function(i) {
         if (d$type[i] == "normal") {

mikelove

nrep <- 100
d <- expand.grid(n=c(3,5,10,20), type=c("normal","t"))
d <- d[rep(seq_len(nrow(d)),each=nrep),]

simulate <- function(n, type) {
  if (type == "normal") {
    dat <- rnorm(n)
  } else {
      dat <- rt(n, df=3)
    }

mikelove

d <- expand.grid(n=c(3,5,10,20), type=c("normal","t"))
d <- d[rep(seq_len(nrow(d)),each=nrep),]

simulate <- function(n, type) {
  if (type == "normal") {
    dat <- rnorm(n)
  } else {
      dat <- rt(n, df=3)
    }
  dat

mikelove

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
g <- keys(txdb, "GENEID")
df <- select(txdb, keys=g, keytype="GENEID", columns="TXID")
ebt <- exonsBy(txdb, by="tx")
set.seed(1)
random.genes <- sample(g, 500, replace=FALSE)
res <- sapply(random.genes, function(gene) {
         txs <- df$TXID[df$GENEID == gene]
         if (length(txs) == 1) return(NA)

mikelove

echo 22 > /sys/class/backlight/acpi_video0/brightness

mikelove

n <- 50
m <- matrix(1:(50000*n),ncol=n)
f <- factor(rep(1:25000,each=2))
system.time({ z <- do.call(rbind, by(m, f, colSums)) })
# 16.3 seconds

library(dplyr)
d <- as.data.frame(cbind(f,m))
system.time({ d %>% group_by(f) %>% summarize_each(funs(sum)) })
# 0.137 seconds