cheuerde

# Claas Heuer, October 2015
#
# Using Sagemath (http://www.sagemath.org/) for calculus
# and analytical maximum likelihood solutions
#
# Examples taken from: https://en.wikipedia.org/wiki/Maximum_likelihood

#########
### a binomial example
#########

cheuerde

# Claas Heuer, October 2015
#
# This is a multivariate extension of the approach described in 
# Kang et.al (2008) and implemented by Zhou and Stephens (2014)
# in 'GEMMA'. The diagonalization of V in a mixed model with 
# one (!) random effect can easily be extended to the multivariate 
# case using Asreml or MCMCglmm.
# This paper made me do this again: http://biorxiv.org/content/biorxiv/early/2015/09/18/027201.full.pdf 
# In that paper they made a kind of unfair comparison to Asreml by not using the
# transformed model, which has de facto nothing to do with Asreml. Here I show how

cheuerde

# this is the conditional model, but does not seem to work well
library(BGLR)

data(wheat)
L <- t(chol(wheat.A))
y <- wheat.Y[,1]

# regression

cheuerde

# taken from here:
# http://stats.stackexchange.com/a/135725/90936

set.seed(7)

#create samples
sample.1 <- rnorm(8, 100, 3)
sample.2 <- rnorm(10, 103, 7)

#we need a pooled data set for estimating parameters in the prior.

cheuerde

# Claas Heuer, September 2015

glht.table <- function(x) {

# I took this from somewehre, but cant remember the source (probably SO))
  pq <- summary(x)$test
  mtests <- cbind(pq$coefficients, pq$sigma, pq$tstat, pq$pvalues)
  error <- attr(pq$pvalues, "error")
  pname <- switch(x$alternativ, less = paste("Pr(<", ifelse(x$df ==0, "z", "t"), ")", sep = ""), 
  greater = paste("Pr(>", ifelse(x$df == 0, "z", "t"), ")", sep = ""), two.sided = paste("Pr(>|",ifelse(x$df == 0, "z", "t"), "|)", sep = ""))

cheuerde

# Claas Heuer, September 2015
#
# Multinomial class probabilities given pairwise probabilities
# from bianry calssifiers.
# Reference: Probability Estimates for Multi-class Classification by Pairwise Coupling, Wu et al. (2003)

# number of classes
K <- 4

# some random probabilities from binary contrasts

cheuerde

# Claas Heuer, September 2015
#
# Simple and stupid pedigree simulation. just for software tests


ped_sim <- function(n = 1000, n_founders = 100, generations  = 5) {

  ped <- data.frame(id = 1:n, sire = as.integer(NA), dam = as.integer(NA), sex = as.integer(NA), gen = as.integer(NA))

# the founders

cheuerde

  # Credit: Taken from:  http://stackoverflow.com/questions/1358003/tricks-to-manage-the-available-memory-in-an-r-session
# improved list of objects
.ls.objects <- function (pos = 1, pattern, order.by,
                        decreasing=FALSE, head=FALSE, n=5) {
    napply <- function(names, fn) sapply(names, function(x)
                                         fn(get(x, pos = pos)))
    names <- ls(pos = pos, pattern = pattern)
    obj.class <- napply(names, function(x) as.character(class(x))[1])
    obj.mode <- napply(names, mode)
    obj.type <- ifelse(is.na(obj.class), obj.mode, obj.class)

cheuerde

# Claas Heuer, September 2015
#
# References: 
# 1) Rohan Fernando, personal communication (2013)
# 2) Wang, Rutledge and Gianola, 1993 GSE
# 3) Linear Models for the Prediction of Animal Breeding Values, Mrode


'gibbs_ginv' <-function(y, Z = NULL, Ainv, scale_a = 0, df_a = -2, scale_e = 0, df_e = -2, niter = 1000, burnin = 500, h2 = 0.3) {

cheuerde

# Claas Heuer, September 2015

# JDBC
sudo apt-get install libpostgresql-jdbc-java

# get MSSQL driver from Microsoft

wget http://download.microsoft.com/download/D/6/A/D6A241AC-433E-4CD2-A1CE-50177E8428F0/1033/sqljdbc_3.0.1301.101_enu.tar.gz

# unpack