using_clipp.R
In clipp: Calculating Likelihoods by Pedigree Paring

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)

## -----------------------------------------------------------------------------
library(clipp)
MAF <- 0.1
geno_freq <- geno_freq_monogenic(p_alleles = c(1 - MAF, MAF))
trans <- trans_monogenic(n_alleles = 2)

## -----------------------------------------------------------------------------
geno_freq_monogenic(p_alleles = c(1 - MAF, MAF), annotate = TRUE)

## -----------------------------------------------------------------------------
trans_monogenic(n_alleles = 2, annotate = TRUE)

## -----------------------------------------------------------------------------
data("dat_small", "penet_small", "dat_large", "penet_large")
head(dat_small)
head(penet_small)

## -----------------------------------------------------------------------------
monozyg_small <- list(c("ora024", "ora027"), c("aey063", "aey064"))

## -----------------------------------------------------------------------------
pedigree_loglikelihood(dat_small, geno_freq, trans, penet_small, 
                       monozyg = monozyg_small, sum_loglik = FALSE, ncores = 2)

## ---- eval = FALSE------------------------------------------------------------
#  system.time(ll <- pedigree_loglikelihood(dat_large, geno_freq, trans, penet_large))
#  #>    user  system elapsed
#  #>   10.64    0.15   10.83
#  ll
#  #> [1] -18020.99

## -----------------------------------------------------------------------------
head(dat_small)

## -----------------------------------------------------------------------------
penet.fn <- function(i, logodds1, logodds2) {
  prob1 <- 1/(1 + exp(-logodds1))
  prob2 <- 1/(1 + exp(-logodds2))
  penet.i <- c(prob1, prob1, prob2)
  if (dat_small$aff[i] == 0)  penet.i <- 1 - penet.i
  return(penet.i)
}

## ---- eval = FALSE------------------------------------------------------------
#  library(stats4)
#  minusll <- function(logodds1 = 0, logodds2 = 0) {
#    penet <- t(sapply(1:nrow(dat_small), penet.fn, logodds1, logodds2))
#    loglik <- pedigree_loglikelihood(dat_small, geno_freq, trans, penet,
#                                     monozyg = monozyg_small, ncores = 2)
#    return(-loglik)
#  }
#  minusll()
#  #> [1] 705.6238
#  fit <- mle(minusll)
#  summary(fit)
#  #> Maximum likelihood estimation
#  #>
#  #> Call:
#  #> mle(minuslogl = minusll)
#  #>
#  #> Coefficients:
#  #>           Estimate Std. Error
#  #> logodds1 -1.359962  0.1054336
#  #> logodds2 -1.313467  6.8597364
#  #>
#  #> -2 log L: 1030.901

## -----------------------------------------------------------------------------
head(dat_small)

## -----------------------------------------------------------------------------
penet.fn <- function(i, logodds1, logodds2) {
  prob1 <- 1/(1 + exp(-logodds1))
  prob2 <- 1/(1 + exp(-logodds2))
  penet.i <- c(prob1, prob1, prob2)
  if (dat_small$aff[i] == 0)  penet.i <- 1 - penet.i
  if (dat_small$geno[i] == "1/1")  penet.i[-1] <- 0     ###
  if (dat_small$geno[i] == "1/2")  penet.i[-2] <- 0     ###
  if (dat_small$geno[i] == "2/2")  penet.i[-3] <- 0     ###
  return(penet.i)
}

## ---- eval = FALSE------------------------------------------------------------
#  minusll()
#  #> [1] 788.5003
#  fit <- mle(minusll)
#  summary(fit)
#  #> Maximum likelihood estimation
#  #>
#  #> Call:
#  #> mle(minuslogl = minusll)
#  #>
#  #> Coefficients:
#  #>           Estimate Std. Error
#  #> logodds1 -1.357596 0.08405912
#  #> logodds2 -1.395321 0.61632248
#  #>
#  #> -2 log L: 1196.65
#  

## -----------------------------------------------------------------------------
head(dat_small)
genotype_probabilities(target = "ora008", dat_small, geno_freq, trans, 
                       penet_small, monozyg_small)

## -----------------------------------------------------------------------------
penet.fn <- function(i) {
  penet.i <- rep(1, 3)
  if (dat_small$geno[i] == "1/1")  penet.i[-1] <- 0
  if (dat_small$geno[i] == "1/2")  penet.i[-2] <- 0
  if (dat_small$geno[i] == "2/2")  penet.i[-3] <- 0
  return(penet.i)
}
penet <- t(sapply(1:nrow(dat_small), penet.fn))
genotype_probabilities(target = "ora008", dat_small, geno_freq, trans, 
                       penet, monozyg_small)