## author: Hae Kyung Im
## profile likelihood for estimating effect of gene expression on drug effect
## IC50 ~ gene.exp with trio covariance structure
## 
## profile restricted likelihood

## input Y (e.g. IC50), M (3m by q matrix)
## m = number of trios
## q number of explanatory variables including constant 1
## rho is the correlation between parent and child

pRL = function(rho, Y, M, n, m, q)
{
##   rho2 = rho^2
##   iRho3 = matrix( c(1-rho2, -rho, rho2, -rho, 1, -rho, rho2, -rho, 1-rho2),
##     3,3, byrow=T)
##   iRho3 = iRho3 / (1-2*rho2) 
##   iRho = kronecker(diag(m), iRho3)

  iRho = iRhofun(rho, m)
  logdetRho = m * log(1 - 2*rho^2)

  Wrho = t(M) %*% iRho %*% M
  iWrho = solve(Wrho)

  Qrho = diag(n) - M %*% iWrho %*% t(M) %*% iRho

  res = - (n-q)/2 * log( t(Y) %*% iRho %*% Qrho %*% Y ) - 1/2 * logdetRho +
    - 1/2 * log( det(Wrho) )

  return(res)
}

## restricted likelihood. Use this for calculating standard errors of estimates rhohat and sig2hat

RL = function(rho, sig2, Y, M, n, m, q)
{
  iRho = iRhofun(rho, m)
  logdetRho = m * log(1 - 2*rho^2)

  Wrho = t(M) %*% iRho %*% M
  iWrho = solve(Wrho)

  Qrho = diag(n) - M %*% iWrho %*% t(M) %*% iRho

  res = - (n-q)/2 * log(sig2) - 1/2 * logdetRho - 1/2 * log( det(Wrho) ) +
    t(Y) %*% iRho %*% Qrho %*% Y/ (2*sig2)
  return(res)
}

## ener function to minimize ener = - loglik

ener = function(rho,Y,M,n,m,q)
  {
    - pRL(rho,Y,M,n,m,q)
  }

## sig2hat maximizes RL for given rho

sig2fun = function(rho, Y, M, n, m, q)
  {

##   rho2 = rho^2
##   iRho3 = matrix( c(1-rho2, -rho, rho2, -rho, 1, -rho, rho2, -rho, 1-rho2),
##     3,3, byrow=T)
##   iRho3 = iRho3 / (1-2*rho2) 
##   iRho = kronecker(diag(m), iRho3)

    iRho = iRhofun(rho,m)
    Wrho = t(M) %*% iRho %*% M
    iWrho = solve(Wrho)
    Qrho = diag(n) - M %*% iWrho %*% t(M) %*% iRho
    res = t(Y) %*% iRho %*% Qrho %*% Y / (n-q)

  return(res)
  }

## betahat maximizes likelihood given sig2 and rho

betafun = function(rho, Y, M, n, m, q)
  {
    iRho = iRhofun(rho, m)  
    Wrho = t(M) %*% iRho %*% M
    iWrho = solve(Wrho)
    res = iWrho %*% t(M) %*% iRho %*% Y 
  }

iRhofun = function(rho,m)
{
  rho2 = rho^2
  iRho3 = matrix( c(1-rho2, -rho, rho2, -rho, 1, -rho, rho2, -rho, 1-rho2),
    3,3, byrow=T)
  iRho3 = iRho3 / (1-2*rho2) 
  res = kronecker(diag(m), iRho3)
  return(res)
}