# Copyright (C) 2017 Francesca Gasperoni <francesca.gasperoni@polimi.it>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Computing standard errors
#
# @param ncovs number of covariates
# @param w: vector of frailties, Kx1
# @param p: vector of proportions, Kx1
# @param D: vector of death per group, Hx1
# @param sumhaz: cumulative baseline
# @param dsumhaz: first derivative of cumulative master
# @param d2sumhaz: second derivative of cumulative master
# @param lamsum:
exact_sec_der <- function( ncovs, w, p, D, sumhaz, dsumhaz, d2sumhaz, lamsum)
{
K <- length( w )
H <- length( D )
#to build
d2p <- array(dim=c(H, K-1, K-1))
d2w <- array(dim=c(H, K, K))
d2cov <- array(dim=c(H, ncovs, ncovs))
d2pw <- array(dim=c(H, K-1, K))
d2pcov <- array(dim=c(H, K-1, ncovs))
d2wcov <- array(dim=c(H, K, ncovs))
exact_der <- matrix(0, nrow=(K-1+K+ncovs), ncol=(K-1+K+ncovs))
if( K > 1 )
{
w_sum_haz = matrix( unlist( lapply( 1:K, function(i) w[i] * sumhaz )), nrow = H, ncol= K)
denominator = rowSums(exp( matrix( unlist( lapply( 1:K, function(i) ( log(p[i]) + D*log(w[i])) )), nrow = H, ncol= K) -w_sum_haz))
numerator1 = rowSums(exp( matrix( unlist( lapply( 1:K, function(i) ( log(p[i]) + (D+1)*log(w[i])) )), nrow = H, ncol= K) -w_sum_haz))
numerator2 = rowSums(exp( matrix( unlist( lapply( 1:K, function(i) ( log(p[i]) + (D+2)*log(w[i])) )), nrow = H, ncol= K) -w_sum_haz))
}else{
denominator = numerator1 = numerator2 = exp(-sumhaz)
}
for( j in 1:H)
{
if(K>=2){
for( g in 1:(K-1) )
{
#Pure term 1 #dim(dpdp) = (K-1)x(K-1)
for( l in 1:(K-1) )
{
d2p[ j, g, l] = -(exp(D[j]*log(w[g])-w[g]*sumhaz[j])-exp(D[j]*log(w[K])-w[K]*sumhaz[j]))*(exp(D[j]*log(w[l])-w[l]*sumhaz[j])-exp(D[j]*log(w[K])-w[K]*sumhaz[j]))/(denominator[j])^2
}
#Mixed term 1-2 #dim(d2pw) = (K-1)xK
for(l in 1:K)
{
d2pw[ j, g, l ] = -(exp(D[j]*log(w[g])-w[g]*sumhaz[j])-exp(D[j]*log(w[K])-w[K]*sumhaz[j]))*(p[l]*exp((D[j]-1)*log(w[l])-w[l]*sumhaz[j])*(D[j]-w[l]*sumhaz[j]))/(denominator[j])^2+
+ ifelse( l == g, exp(D[j]*log(w[l])-w[l]*sumhaz[j])*(D[j]-w[l]*sumhaz[j])/denominator[j], 0) +
+ ifelse( l == K, exp(D[j]*log(w[l])-w[l]*sumhaz[j])*(D[j]-w[l]*sumhaz[j])/denominator[j], 0)
}
#Mixed term 1-3 #dim(d2pcov) = (K-1)xncovs
for(r in 1:ncovs)
{
d2pcov[ j, g, r ] = dsumhaz[j,r] * ( (exp((D[j]+1)*log(w[K])-w[K]*sumhaz[j]) - exp((D[j]+1)*log(w[g])-w[g]*sumhaz[j]))/denominator[j] +
-(exp(D[j]*log(w[K])-w[K]*sumhaz[j]) - exp(D[j]*log(w[g])-w[g]*sumhaz[j]))*numerator1[j]/(denominator[j])^2 )
}
}
}
#Pure term 2 #dim(d2w) = KxK
for( g in 1:K )
{
for( l in 1:K )
{
d2w[ j, g, l] = -p[g]*p[l]*exp((D[j]-1)*log(w[g]*w[l])-(w[g]+w[l])*sumhaz[j])*(D[j]-w[g]*sumhaz[j])*(D[j]-w[l]*sumhaz[j])/(denominator[j])^2+
+ ifelse(g!=l, 0, p[g]*exp((D[j]-2)*log(w[g])-w[g]*sumhaz[j])*((D[j]-1-w[g]*sumhaz[j])*(D[j]-w[g]*sumhaz[j])-w[g]*sumhaz[j])/denominator[j])
}
}
#Pure term 3 #dim(d2cov) = ncovs x ncovs
for(r in 1:ncovs)
{
for(s in 1:ncovs)
{
d2cov[ j, r, s ] = lamsum[ j, r, s ]+ numerator2[j]/denominator[j]*(dsumhaz[j,r]*dsumhaz[j,s])+
- numerator1[j]/denominator[j]*d2sumhaz[j,r,s] +
- (numerator1[j]/denominator[j])^2*(dsumhaz[j,r]*dsumhaz[j,s])
}
}
#Mixed term 2-3 #dim(d2wcov) = Kxncovs
for(g in 1:K)
{
for(r in 1:ncovs)
{
d2wcov[ j, g, r ] = p[g] * exp( (D[j]-1)*log(w[g]) - w[g]*sumhaz[j]) * dsumhaz[j,r] * ( (-w[g]*(D[j]+1-w[g]*sumhaz[j]))/denominator[j]+
+ ((D[j]-w[g]*sumhaz[j])*numerator1[j])/(denominator[j])^2)
}
}
if( K > 2 ){
exact_dertmp = rbind( cbind(d2p[j, , ], d2pw[j, , ], d2pcov[j, , ]),
cbind(t(d2pw[j, , ]), d2w[j, , ], d2wcov[j, , ]),
cbind(t(d2pcov[j, , ]), t(d2wcov[j, , ]), d2cov[j, , ]) )
}else if( K == 2 ){
exact_dertmp = rbind( c(d2p[j, , ], d2pw[j, , ], d2pcov[j, , ]),
cbind(d2pw[j, , ], d2w[j, , ], d2wcov[j, , ]),
cbind(d2pcov[j, , ], t(d2wcov[j, , ]), d2cov[j, , ]) )
}else if( K == 1 ){
exact_dertmp = rbind( c(d2w[j, , ], d2wcov[j, , ]), cbind( d2wcov[j, , ], d2cov[j, , ] ) )
}
exact_der = exact_der + exact_dertmp
}
return(exact_der)
}
# @param ncovs number of covariates
# @param w: vector of frailties, Kx1
# @param p: vector of proportions, Kx1
# @param D: vector of death per group, Hx1
# @param sumhaz: Lambda_0*exp(X_ij* beta), H x 1
# @param dsumhaz: first derivative of cumulative baseline haz, H x ncovs
# @param d2sumhaz: second derivative of cumulative baseline haz, H x ncovs
# @param lamroutlam0: lambda_0r/lambda_0 + X_ijr, H x ncovs
# @param deltaijcovr:
# @param dsumhazred:
exact_first_der <- function( ncovs, w, p, D, sumhaz, dsumhaz, lamroutlam0, deltaijcovr, dsumhazred )
{
K <- length( w ) #number of latent populations
H <- length( D ) #number of groups
#to build
dp <- rep( 0, K-1 )
dw <- rep( 0, K )
dcov <- rep( 0, ncovs )
dcov2 <- rep( 0, ncovs )
exact_derivative1 <- exact_derivative2 <- rep( 0, K-1+K+ncovs )
if( K > 1 )
{
w_sum_haz <- matrix( unlist( lapply( 1:K, function(i) w[i] * sumhaz )), nrow = H, ncol= K)
denominator <- rowSums(matrix( unlist( lapply( 1:K, function(i) p[i]*w[i]^D )), nrow = H, ncol= K)*exp(-w_sum_haz))
numerator1 <- rowSums(matrix( unlist( lapply( 1:K, function(i) p[i]*w[i]^(D+1) )), nrow = H, ncol= K)*exp(-w_sum_haz))
numerator2 <- rowSums(matrix( unlist( lapply( 1:K, function(i) p[i]*w[i]^(D+2) )), nrow = H, ncol= K)*exp(-w_sum_haz))
}else{
denominator <- numerator1 <- numerator2 <- exp(-sumhaz)
}
for( j in 1:H )
{
for(g in 1:(K-1))
{
dp[g] <- (w[g]^(D[j])*exp(-w[g]*sumhaz[j]) - w[K]^(D[j])*exp(-w[K]*sumhaz[j]))/denominator[j]
}
for(g in 1:K)
{
dw[g] <- (p[g] * w[g]^(D[j]-1) * exp(-w[g]*sumhaz[j]) * (D[j]-w[g]*sumhaz[j]))/denominator[j]
}
for(r in 1:ncovs)
{
dcov[r] <- -numerator1[j]*dsumhaz[j,r]/denominator[j] + lamroutlam0[j,r]
dcov2[r] <- -numerator1[j]*dsumhazred[j,r]/denominator[j] + deltaijcovr[j,r]
}
exact_derivative1 <- exact_derivative1 + c( dp, dw, dcov)
exact_derivative2 <- exact_derivative2 + c( dp, dw, dcov2)
}
return( list( exact_derivative1, exact_derivative2 ) )
}
# @param ncovs: number of covariates
# @param alpha: matrix of posterior density of Z, HxK
# @param w: vector of frailties, Kx1
# @param p: vector of p, Kx1
# @param D: number of events in each group, Hx1
# @param sumhaz: Lambda_0*exp(X_ij* beta), H x 1
# @param lamroutlam0: lambda_0r/lambda_0 + X_ijr, H x ncovs
# @param dsumhaz: first derivative of cumulative baseline haz, H x ncovs
Iy <- function( ncovs, alpha, w, p, D, sumhaz, lamroutlam0, dsumhaz)
{
K <- dim(alpha)[2] #number of latent populations
H <- dim(alpha)[1] #number of groups
#to build
dpdp <- array(dim=c(H, K-1))
dwdw <- array(dim=c(H, K))
dcovdcov <- array(dim=c(H, ncovs))
expectedStexpectedSj <- matrix(0, nrow=(K-1+K+ncovs), ncol=(K-1+K+ncovs))
#each array is group specific
for( j in 1:H)
{
#Pure term 1 #dim(dpdp) = (K-1)xq
dpdp[ j, ] <- alpha[ j, 1:(K-1) ]/p[1:(K-1)] - alpha[ j, K ]/p[K]
#Pure term 2 #dim(dwdw) = K x K
dwdw[ j, ] <- alpha[ j, 1:K]*(D[j]/w[1:K]-sumhaz[j])
#Pure term 3 #dim(dcovdcov) = ncovs x ncovs
for(r in 1:ncovs) {
dcovdcov[ j, r ] <- alpha[ j, ] %*% t(lamroutlam0[j,r]-w*dsumhaz[j,r])
}
expectedStexpectedSj <- expectedStexpectedSj + c( dpdp[j,], dwdw[j, ], dcovdcov[j, ] ) %*% t(c( dpdp[j,], dwdw[j, ], dcovdcov[j, ] ))
}
return(expectedStexpectedSj)
}
# @param ncovs: number of covariates
# @param alpha: matrix of posterior density of Z, HxK
# @param w: vector of frailties, Kx1
# @param p: vector of p, Kx1
# @param D: number of events in each group, Hx1
# @param sumhaz: Lambda_0*exp(X_ij* beta), H x 1
# @param lamroutlam0: lambda_0r/lambda_0 + X_ijr, H x ncovs
# @param dsumhaz: first derivative of cumulative baseline haz, H x ncovs
Ixy = function( ncovs, alpha, w, p, D, sumhaz, lamroutlam0, dsumhaz )
{
K <- dim(alpha)[2] #number of latent populations
H <- dim(alpha)[1] #number of groups
#to build
dpdp <- array(dim=c(H, K-1, K-1))
dwdw <- array(dim=c(H, K, K))
dcovdcov <- array(dim=c(H, ncovs, ncovs))
dpdw <- array(dim=c(H, K-1, K))
dpdcov <- array(dim=c(H, K-1, ncovs))
dwdcov <- array(dim=c(H, K, ncovs))
expectedStSgroupj <- matrix(0, nrow=(K-1+K+ncovs), ncol=(K-1+K+ncovs))
for( j in 1:H)
{
#Pure term 1 #dim(dpdp) = (K-1)x(K-1)
dpdp[ j, , ] <- alpha[ j, K ]/p[K]^2
if( K > 2){
diag( dpdp[ j, , ] ) <- diag( dpdp[ j, , ] ) + alpha[ j, 1:(K-1) ]/p[1:(K-1)]^2
}else if( K <= 2){
dpdp[ j, , ] <- alpha[ j, K ]/p[K]^2 + alpha[ j, 1:(K-1) ]/p[1:(K-1)]^2
}
#Pure term 2 #dim(dwdw) = K x K
if(K>1){
dwdw[ j, , ] <- diag(alpha[ j, 1:K]*(D[j]/w[1:K]-sumhaz[j])^2)#( lamwoutlam0[ j, 1:K ] - dsumhazw[j,1:K] )^2 )
}else if(K==1){
dwdw[ j, , ] <- alpha[ j, 1:K]*(D[j]/w[1:K]-sumhaz[j])^2#( lamwoutlam0[ j, 1:K ] - dsumhazw[j,1:K] )^2 )
}
#Pure term 3 #dim(dcovdcov) = ncovs x ncovs
for(r in 1:ncovs) {
for(s in 1:ncovs) {
dcovdcov[ j, r, s ] <- alpha[j,] %*% t( (lamroutlam0[j,r]-w*dsumhaz[j,r])*(lamroutlam0[j,s]-w*dsumhaz[j,s]) )
}
}
#Mixed term 1-2 #dim(dpdw) = (K-1)xK
dpdwtmp <- rep(0,(K-1))
dpdwtmp[K-1] <- -alpha[ j,K]/p[K] * (D[j]/w[K]-sumhaz[j])
if(K>2){
dpdw[j, , ] <- cbind( diag(alpha[ j,1:(K-1)]/p[1:(K-1)] * (D[j]/w[1:(K-1)]-sumhaz[j])), dpdwtmp )
}else{
dpdw[j, , ] <- cbind( alpha[ j,1:(K-1)]/p[1:(K-1)] * (D[j]/w[1:(K-1)]-sumhaz[j]), dpdwtmp )
}
#Mixed term 1-3 #dim(dpdcov) = K-1 x ncovs
for (r in 1:ncovs)
{
dpdcov[j, , r] <- alpha[j,1:(K-1)]/p[1:(K-1)] * (lamroutlam0[j,r]-w[1:(K-1)]*dsumhaz[j,r]) - alpha[j,K]/p[K]*(lamroutlam0[j,r]-w[K]*dsumhaz[j,r])
}
#Mixed term 2-3 #dim(dpdcov) = K x ncovs
for (r in 1:ncovs)
{
dwdcov[j, , r] <- alpha[j,1:K] * (lamroutlam0[j,r]-w[1:K]*dsumhaz[j,r]) * (D[j]/w[1:K]-sumhaz[j])
}
if(K>2){
expectedStSgroupjtmp <- rbind( cbind(dpdp[j, , ], dpdw[j, , ], dpdcov[j, , ]),
cbind(t(dpdw[j, , ]), dwdw[j, , ], dwdcov[j, , ]),
cbind(t(dpdcov[j, , ]), t(dwdcov[j, , ]), dcovdcov[j, , ]) )
}else if( K==2 ){
expectedStSgroupjtmp <- rbind( c(dpdp[j, , ], dpdw[j, , ], dpdcov[j, , ]),
cbind(dpdw[j, , ], dwdw[j, , ], dwdcov[j, , ]),
cbind(dpdcov[j, , ], t(dwdcov[j, , ]), dcovdcov[j, , ]) )
}else if(K==1){
expectedStSgroupjtmp <- rbind( c(dwdw[j, , ], dwdcov[j, , ]), cbind( dwdcov[j, , ], dcovdcov[j, , ] ) )
}
expectedStSgroupj <- expectedStSgroupj + expectedStSgroupjtmp
}
return(expectedStSgroupj)
}
# @param ncovs: number of covariates
# @param alpha: matrix of posterior density of Z, HxK
# @param w: vector of frailties, Kx1
# @param p: vector of proportion, Kx1
# @param D: number of events in each group, Hx1
# @param d2sumhaz: second derivative of cumulative baseline haz, H x ncovs x ncovs
# @param lamsum: second derivative of baseline haz, H x ncovs x ncovs
# @param dsumhaz: first derivative of cumulative baseline haz, H x ncovs
Ix <- function( ncovs, alpha, w, p, D, d2sumhaz, lamsum, dsumhaz)
{
K <- dim(alpha)[2] #number of latent populations
H <- dim(alpha)[1] #number of groups
#pure terms
d2pp <- array( dim = c( H, K-1, K-1 ) )
d2ww <- array( dim = c( H, K, K ) )
d2covcov <- array( dim = c( H, ncovs, ncovs ) )
#mixed terms
d2wcov <- array( dim = c( H, K , ncovs ) )
#result
expectedBgroupj <- matrix(0, nrow = (K-1+K+ncovs), ncol = (K-1+K+ncovs) )
for(j in 1:H)
{
#Pure term 1 #dim(d2pp) = (K-1)x(K-1)
d2pp[ j, , ] <- - alpha[ j, K ]/p[K]^2
if(K>2){
diag( d2pp[ j, , ] ) <- diag( d2pp[ j, , ] ) - alpha[ j, 1:(K-1) ]/p[1:(K-1)]^2
}else{
d2pp[ j, , ] <- - alpha[ j, K ]/p[K]^2 - alpha[ j, 1:(K-1) ]/p[1:(K-1)]^2
}
#Pure term 2 #dim(d2pp) = K x K
if(K>1){
d2ww[ j, , ] <- diag( -alpha[ j, 1:K ]*D[j]/w[1:K]^2 )
}else if(K==1){
d2ww[ j, , ] <- -alpha[ j, 1:K ]*D[j]/w[1:K]^2
}
#Pure term 3 #dim(d2covcov) = ncovs x ncovs
for(r in 1:ncovs)
{
for(s in 1:ncovs)
{
d2covcov[ j, r, s ] <- t(alpha[ j, 1:K]) %*% (lamsum[ j, r, s]- w[1:K] * d2sumhaz[ j , r, s])
}
}
#Mixed term 1 #dim(d2wcov) = K x ncovs
for(k in 1:K)
{
for(r in 1:ncovs)
{
d2wcov[ j, k, r ] <- - alpha[j, k] * dsumhaz[ j, r ] #(lamwbetaoutlam02[j, k, r]-d2sumhazwbeta[j, k, r])
}
}
if(K>=2)
{
expectedBgroupjtmp <- rbind( cbind( d2pp[ j, , ], matrix( 0, nrow = K-1, ncol = K ), matrix( 0, nrow = K-1, ncol = ncovs ) ),
cbind( matrix( 0, nrow = K, ncol = K-1 ), d2ww[j, , ], d2wcov[j, , ] ),
cbind( matrix( 0, nrow = ncovs, ncol = K-1 ), t(d2wcov[j, , ]), d2covcov[j, , ]) )
}else{
expectedBgroupjtmp <- rbind( c( d2ww[j, , ], d2wcov[j, , ] ),
cbind( d2wcov[j, , ], d2covcov[j, , ]) )
}
expectedBgroupj <- expectedBgroupj + expectedBgroupjtmp
}
return( expectedBgroupj )
}
#computation of the information matrix
information_matrix <- function( expectedS, expectedB, expectedSexpectedS )
{
info <- -expectedB-expectedS+expectedSexpectedS
cov <- solve(info)
se <- sqrt(diag(cov))
return( list( cov, se ) )
}
hazard_derivatives <- function( X, risk_index, w_off, beta_hat, N, YY)
{
ncovs <- ncol(X)
nt <- dim(risk_index)[2]
dhaz <- dcumhaz <- matrix( nrow = nt, ncol = ncovs )
d2haz <- d2cumhaz <- array( dim = c( nt, ncovs, ncovs ) )
#nt x ncovs
hazbase <- as.numeric(exp( beta_hat %*% t(X) ))
dYY <- t(risk_index) %*% (w_off * hazbase * X)
dhaz <- - dYY * N/YY^2
if( ncovs > 1 )
{
nobs <- nrow(X)
d2YYtemp <- array( t( apply( X, 1, function(x) x%*%t(x)) ), dim=c( nobs, ncovs, ncovs ) )
#nt x ncovs x ncovs
dYYdYY <- array( t( apply( dYY[,1:ncovs], 1, function(x) x%*%t(x)) ), dim=c( nt, ncovs, ncovs ) )
L <- lapply( seq_len( dim(risk_index)[2] ), function(i) apply(as.vector(w_off*exp(beta_hat %*% t(X)))*d2YYtemp*risk_index[,i],MARGIN = 2:3,sum))
d2YY <- array( unlist(L), dim = c(nrow(L[[1]]), ncol(L[[1]]), length(L)))
d2YY <- aperm( d2YY, c( 3, 1, 2 ))
}else{
d2YYtemp <- X^2
dYYdYY <- dYY^2
d2YY <- t(risk_index) %*% ( w_off*exp( beta_hat * X )*d2YYtemp ) # nt x ncovs
}
# constructing array for point-wise operation
d2haz <- array( N/YY^3*(2*dYYdYY-d2YY*YY), dim = c( nt, ncovs, ncovs ))
if(ncovs>1)
{
dcumhaz <- matrix( apply( dhaz, MARGIN = 2, cumsum ), ncol = ncovs) #sum by column
d2cumhaz <- array( apply( d2haz, MARGIN = 2:3, cumsum ), dim = c( nt, ncovs, ncovs ) )
}else{
dcumhaz <- matrix( cumsum( dhaz ), ncol = ncovs)
d2cumhaz <- array( cumsum( d2haz ), dim = c( nt, ncovs, ncovs ))
}
lis <- list()
lis[['dhaz']] <- dhaz
lis[['d2haz']] <- d2haz
lis[['dcumhaz']] <- dcumhaz
lis[['d2cumhaz']] <- d2cumhaz
return(lis)
}
input_se_functions <- function( groups, time, status, X, beta_hat, haz, cumhaz, dhaz, dcumhaz, d2haz, d2cumhaz )
{
H <- length( unique( groups ) )
ncovs <- length( beta_hat )
lamroutlam0 <- matrix( nrow = H, ncol = ncovs )
sumhaz <- rep( 0, H )
dsumhaz <- array( dim = c( H, ncovs ) )
d2sumhaz <- lamsum <- array( dim = c( H, ncovs, ncovs ) )
deltaijcovr <- matrix( 0, nrow = H, ncol = ncovs )
dsumhazred <- matrix( 0, nrow = H, ncol = ncovs )
for( j in 1:H )
{
hospj <- groups==j
ebz <- exp( as.matrix( X[hospj,,drop=FALSE] ) %*% beta_hat )
ti <- match(time[hospj], cumhaz$time)
lam0t <- cumhaz$hazard[ti]
sumhaz[j] <- sum(lam0t * ebz)
for (r in 1:ncovs){
covr <- X[hospj,r]
lamrt <- dcumhaz[ti, r]
#dsumhaz is the second part of III term
dsumhaz[j,r] <- sum(ebz*(covr*lam0t + lamrt))
#new lines for III term
lamroutlam0tmp <- status[hospj]*( dhaz[ti,r]/haz[ti] + covr )
lamroutlam0[j,r] <- sum(lamroutlam0tmp[haz[ti]>0])
deltaijcovrtmp <- status[hospj]*covr
deltaijcovr[j,r] <- sum(deltaijcovrtmp[haz[ti]>0])
dsumhazred[j,r] <- sum(ebz*(covr*lam0t))
for (s in 1:ncovs){
covs <- X[hospj,s]
lamst <- dcumhaz[ti, s]
lamrst <- d2cumhaz[ti, r, s]
d2sumhaz[j,r,s] <- sum(ebz*(covr*covs*lam0t + covr*lamst + covs*lamrt + lamrst))
lamtmp <- status[hospj] * (d2haz[ti,r,s]*haz[ti] - dhaz[ti,r]*dhaz[ti,s])/haz[ti]^2
lamsum[j,r,s] <- sum(lamtmp[haz[ti]>0]) # exclude zero hazard
}
}
}
lis <- list()
lis[['lamroutlam0']] <- lamroutlam0
lis[['lamsum']] <- lamsum
lis[['deltaijcovr']] <- deltaijcovr
lis[['sumhaz']] <- sumhaz
lis[['dsumhazred']] <- dsumhazred
lis[['dsumhaz']] <- dsumhaz
lis[['d2sumhaz']] <- d2sumhaz
return(lis)
}
log_likelihood <- function( groups, time, status, X, p, w, alpha, beta_hat, haz, cumhaz)
{
H <- dim(alpha)[1]
K <- dim(alpha)[2]
ncovs <- ncol(X)
llik = 0
for( j in 1:H )
{
hospj <- groups==j
hospj_lam <- groups==j & status==1
ti <- match(time[hospj], cumhaz$time)
ti_lam <- match(time[hospj_lam], cumhaz$time)
for( k in 1:K )
{
llik_LAM <- ifelse( ncovs > 0,
sum( cumhaz$hazard[ti]*w[k]*exp( as.matrix( X[hospj,,drop=FALSE] ) %*% beta_hat ) ),
sum( cumhaz$hazard[ti]*w[k] ) )
#w_off[hospj][1] because w_off[hospj] are all equal
llik_lam <- ifelse( ncovs > 0,
sum( log( haz[ti_lam]*w[k]*exp( as.matrix( X[hospj_lam,,drop=FALSE] ) %*% beta_hat ) ) ),
sum( log( haz[ti_lam]*w[k] ) ) )
llik_p <- ifelse( K==1, 1, p[ k ] )
llik <- llik + alpha[j,k]*(llik_lam-llik_LAM) + alpha[j,k]*log( llik_p )
}
}
return(llik)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.