TABLE OF CONTENTS

estimation/updateregprof [ Functions ]

NAME

    updateregprof --- update regression coefficients (profile likelihood)

FUNCTION

Update regression coefficients beta1 and beta1 by maximizing the condition profile loglikelihood given all the other parameters. Give the regression coefficients, the baseline hazard coefficients are then updated.

SYNOPSIS

1876 updateregprof <- function(m, Ji, datamat, datamatd, alphars, alpharsd, as, asd, 
1877     betahat, betadhat, K, Kd, frailtyoutput, dispersionoutput, smooth)

INPUTS

    m              number of clusters
    Ji             cluster sizes
    datamat        data matrix generated by makedatamat for event 1
    datamatd       data matrix generated by makedatamat for event 2
    alphars            matrix of baseline hazard parameters for event 1
    alpharsd           matrix of baseline hazard parameters for event 2
    as                 matrix of discretization breakpoints for event 1
    asd                matrix of discretization breakpoints for event 2
    betahat            regression coefficient estimates for event 
    betadhat           regression coefficient estimates for event 2
    K                  number of discretization intervals for event 1
    Kd                 number of discretization intervals for event 2
    frailtyoutput  output from fupdatefrailties4
    dispersionoutput  output from one of the dispersion estimation routines
    smooth         boolean, indicating whether hazard estimates should be smoothed

OUTPUTS

    betahat            regression coefficient estimates for event 
    betadhat           regression coefficient estimates for event 2
    alphars            matrix of baseline hazard parameters for event 1
    alpharsd           matrix of baseline hazard parameters for event 2
    as                 matrix of discretization breakpoints for event 1
    asd                matrix of discretization breakpoints for event 2
    loglik         loglikelihood ( = loglik1 + loglik2 )
    loglik1        profile loglikelihood of betahat
    loglik2        profile loglikelihood of betadhat

SOURCE

1880 {
1881     
1882     # Prepare vectors and matrices for direct submission into Fortran, 
1883     # without requiring type conversions
1884     Uijhatframe <- frailtyoutput$Uijframe
1885     fUijmat <- as.double(Uijhatframe$Uijmat)  # Frailty matrices
1886     fVijmat <- as.double(Uijhatframe$Vijmat)    
1887     fas = as.double(as)    # Baseline hazard parameters
1888     fasd = as.double(asd)
1889     # Covariates for event 1
1890     Z <- as.double(matrix(datamat[, -(1:8)], dim(datamat)[1], dim(datamat)[2] - 8)) 
1891     ncovs <- as.integer(length(betahat))   # number of covariates
1892     d1 <- as.integer(dim(datamat)[1])     # data matrix dimension (event 1)
1893     index <- as.integer(datamat[, c("i", "j", "k", "r", "smin", "smax")])   
1894     delta <- as.double(datamat[, "delta"])     # recurrent event indicators
1895     times <- as.double(datamat[, "time"])     # recurrent event times
1896     # Covariates for event 2
1897     Zd <- as.double(matrix(datamatd[, -(1:8)], dim(datamatd)[1], dim(datamatd)[2] - 8)) 
1898     d2 <- as.integer(dim(datamatd)[1])     # data matrix dimension (event 2)
1899     indexd <- as.integer(datamatd[, c("i", "j", "k", "r", "smin", "smax")])  
1900     deltad <- as.double(datamatd[, "delta"])  # terminal event indicators
1901     timesd <- as.double(datamatd[, "time"])   # terminal event times
1902     nr = as.integer(dim(as)[1])    # Maximum value of r (number of strata)
1903     ns = as.integer(dim(as)[2])    # Maximum value of s (recurrent) + 1
1904     nsd = as.integer(dim(asd)[2])  # Maximum value of s (terminal)  + 1 
1905 
1906     # Update regression parameters betahat, betadhat using the profile likelihood. 
1907     # Use the built-in  optim method with the BFGS algorithm.
1908     # The loglikelihood is computed via fmkproflik2and the gradient 
1909     # via fmkprofgr2. Conditionally on the frailties, the processes are 
1910     # independent and can be fit separately 
1911     opt <- try(optim(betahat, fn = fmkproflik2, gr = fmkprofgr2,
1912         control = list(fnscale=-1), method = "BFGS",
1913         m = m, Ji = Ji, index = index, Z = Z, delta = delta, times = times, as = fas, 
1914         Uijmat = fUijmat, ncovs = ncovs, d = d1, nr = nr, ns = ns))
1915     optd <- try(optim(betadhat, fn = fmkproflik2, gr = fmkprofgr2, 
1916     control = list(fnscale=-1), method = "BFGS", 
1917         m = m, Ji = Ji, index = indexd, Z = Zd, delta = deltad, times = timesd, 
1918         as = fasd, Uijmat = fVijmat, ncovs = ncovs, d = d2, nr = nr, ns = nsd))
1919     # Check for errors and quit if it didn't converge
1920     if(inherits(opt, "try - error") | inherits(optd, "try - error")){
1921         warning("Inner loop did not converge")
1922           browser()
1923         return(0)
1924     }    
1925     if(opt$convergence != 0 | optd$convergence != 0){
1926         warning("Inner loop did not converge")
1927           browser()
1928         return(0)
1929     }       
1930     betahat <- opt$par
1931     betadhat <- optd$par
1932     loglik1 <- opt$value
1933     loglik2 <- optd$value
1934     loglik <- loglik1 + loglik2
1935     # Given the updated regression parameters, compute the 
1936     # corresponding baseline hazard parameters
1937     alphars <- fmakealphars2(m, Ji, datamat, betahat, as, fUijmat, smooth)
1938     alpharsd <- fmakealphars2(m, Ji, datamatd, betadhat, asd, fVijmat, smooth)
1939     
1940     # Format for output
1941     return(list(betahat = betahat, betadhat = betadhat,
1942                 alphars = alphars, alpharsd = alpharsd,
1943                 as = as, asd = asd,
1944                 loglik = loglik, loglik1 = loglik1, loglik2 = loglik2))  
1945 }