R/loglik_frailty.R

In FamEvent: Family Age-at-Onset Data Simulation and Penetrance Estimation

loglik_frailty<- function(X, Y, theta, cuts=NULL, nbase, data, design, base.dist, frailty.dist, agemin, vec=FALSE)
{
  if(base.dist=="lognormal") bparms <- c(theta[1], exp(theta[2]))
  else bparms <- exp(theta[1:nbase])
  
  nX <- dim(X)[2]
  xbeta <- c(X%*%theta[(nbase+1):(nbase+nX)])
  kappa <- exp(theta[length(theta)])
  
  time0 <- Y[,1] - agemin
  cuts0 <- cuts - agemin
  status <- Y[,2]
  
  
  wt <- data$weight
  wt_fam <- wt[!duplicated(data$famID.byuser)]
  
  bhaz <- hazards(base.dist, time0, bparms, cuts=cuts0)
  bcumhaz <- cumhaz(base.dist, time0, bparms, cuts=cuts0)
  H <- bcumhaz*exp(xbeta)
  logh <- log(bhaz) + xbeta
  loglik <-  wt * (status*logh)
  df <- data$df[!duplicated(data$famID.byuser)]
  s <- aggregate(H, list(data$famID.byuser), sum)[,2]
  logdL <- wt_fam*log( dlaplace(frailty.dist, g=s, d=df, k=kappa) )
  
  
  # Ascertainment correction (design = pop, pop+)
  # families with no proband, logasc=0
  # families with affected proband, logasc = log(1-P(T>a))
  # families with unaffected proband, logasc = log(P(T>a))
  
  if(design=="noasc"){
    sloglik <- sum(loglik[loglik!=-Inf], na.rm=T) + sum(logdL[logdL!=-Inf], na.rm=T) 
    loglik[!duplicated(data$famID.byuser)] <- loglik[!duplicated(data$famID.byuser)] + logdL
  } else{
    
    ip <- data$proband == 1
    ip_fam <- aggregate(data$proband, list(data$famID.byuser), sum)[,2] # indicates if family has proband or not
    i_ap <- ifelse(data$proband == 1 & data$time < data$currentage, 1, 0)[ip] ### indicates if proband is affected ###
    
    cagep <- data$currentage[ip]-agemin
    xbeta.p <- xbeta[ip]
    bcumhaz.p <- cumhaz(base.dist, cagep, bparms, cuts=cuts0)
    laplace.p <- laplace(frailty.dist, bcumhaz.p*exp(xbeta.p), kappa)
    logasc.p <- ifelse(i_ap==1, log(1-laplace.p), log(laplace.p))
    logasc <- wt_fam*ifelse(ip_fam==0, 0, logasc.p)
    #ip_fam[ip_fam!=0] <- i_ap
    
    #print(c(sum(ip), sum(ip_fam==1)))
    logasc[logasc == -Inf] <- 0
    sloglik <- sum(loglik[loglik!=-Inf], na.rm=T) + sum(logdL[logdL!=-Inf], na.rm=T) - sum(logasc[logasc!=-Inf], na.rm=T)
    #loglik[ip] <- loglik[ip] + logdL[ip_fam==1] - logasc[ip_fam==1]
    loglik[!duplicated(data$famID.byuser)] <- loglik[!duplicated(data$famID.byuser)] + logdL
    loglik[ip] <- loglik[ip] - logasc[ip_fam==1]
  }
  
  #print(c(theta, -sloglik))
  #print(c(theta, -sloglik, sum(logdL), sum(logasc)))
  if(vec) return(-loglik)
  else  return(-sloglik)
  
}