./blupsurv/R/bivrec.r

1. fortranwrappers/fmkstderr

fortranwrappers/fmkstderr [ Functions ]

NAME

    fmkstderr --- compute standard error

FUNCTION

Compute the standard error entirely within FORTRAN. This is a wrapper for the FORTRAN routine fmksterr.

SYNOPSIS

2322 fmkstderr <- function(m, Ji, b, datamat, datamatd, alphars, alpharsd, 
2323                     betahat, betadhat, as, asd, frailtyoutput, dispersionoutput, 
2324                     K, Kd, univariate = FALSE)

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
    betahat            regression coefficient estimates for event 
    betadhat           regression coefficient estimates for event 2
    as                 matrix of discretization breakpoints for event 1
    asd                matrix of discretization breakpoints for event 2
    frailtyoutput  output from fupdatefrailties4
    dispersionoutput  output from one of the dispersion estimation routines
    K                  number of discretization intervals for event 1
    Kd                 number of discretization intervals for event 2
    univariate     boolean indicating whether the process is univariate
  OUTPUT
    stderr         a vector of standard errors for regression and hazard params

SOURCE

2327 {
2328     pi <- frailtyoutput$pqrs$pi; qi <- frailtyoutput$pqrs$qi;
2329     ri <- frailtyoutput$pqrs$ri; si <- frailtyoutput$pqrs$si;
2330     wij <- frailtyoutput$pqrs$wij;zij <- frailtyoutput$pqrs$zij;
2331     wi <- frailtyoutput$pqrs$wi
2332     sigma2hat <- dispersionoutput$sigma2hat;
2333     sigma2dhat <- dispersionoutput$sigma2dhat
2334     nu2hat <- dispersionoutput$nu2hat; nu2dhat <- dispersionoutput$nu2dhat
2335     thetahat <- dispersionoutput$thetahat
2336    
2337     ncovs1 <- length(betahat)
2338     ncovs2 <- length(betadhat)
2339     
2340     # Prepare data for submission to Fortran function fmkstderr
2341     Z <- matrix(datamat[, -(1:8)], dim(datamat)[1], dim(datamat)[2] - 8)
2342     d1 <- dim(Z)[1]
2343     index1 <- datamat[, c("i", "j", "k", "r", "smin", "smax")]
2344     times1 <- datamat[, "time"]
2345     Zd <- matrix(datamatd[, -(1:8)], dim(datamatd)[1], dim(datamatd)[2] - 8)
2346     d2 <- dim(Zd)[1]
2347     index2 <- datamatd[, c("i", "j", "k", "r", "smin", "smax")]
2348     times2 <- datamatd[, "time"]
2349     nr = dim(alphars)[1]
2350     ns = dim(alphars)[2]
2351     nsd = dim(alpharsd)[2]
2352     Kcum <- cumsum(c(1, K))
2353     Kdcum <- cumsum(c(1, Kd))
2354     np <- sum(K) + ncovs1
2355     npd <- sum(Kd) + ncovs2
2356     S <- matrix(0, np + npd, np + npd)
2357     # FORTRAN call
2358     out <- try(.Fortran("fmkstderr", Smat = as.double(S), B = as.double(b),
2359         index1 = as.integer(index1), index2 = as.integer(index2),
2360         Z = as.double(Z), Zd = as.double(Zd),
2361         alphars = as.double(alphars), alpharsd = as.double(alpharsd),
2362         as = as.double(as), asd = as.double(asd),
2363         betahat = as.double(betahat), betadhat = as.double(betadhat),
2364         times1 = as.double(times1), times2 = as.double(times2),
2365         pi = as.double(pi), qi = as.double(qi),
2366         ri = as.double(ri), si = as.double(si),
2367         wi = as.double(wi), wij = as.double(wij), zij = as.double(zij),
2368         sig2 = as.double(sigma2hat), sig2d = as.double(sigma2dhat),
2369         nu2 = as.double(nu2hat), nu2d = as.double(nu2dhat), theta = as.double(thetahat),
2370         ncovs1 = as.integer(ncovs1), ncovs2 = as.integer(ncovs2),
2371         nr = as.integer(nr), ns = as.integer(ns), nsd = as.integer(nsd),
2372         np = as.integer(np), npd = as.integer(npd),
2373         d1 = as.integer(d1), d2 = as.integer(d2),
2374         m = as.integer(m), Ji = as.integer(Ji), maxj = as.integer(max(Ji)),
2375         Kcum = as.integer(Kcum), Kdcum = as.integer(Kdcum), DUP = FALSE))
2376     if(!univariate){
2377         # extract standard error from the output
2378         x <- matrix(out$B, np + npd, ncovs1 + ncovs2)
2379         stderr <- sqrt(x[b == 1])
2380     }else{
2381         # For univariates, the matrix has to be solved. This could be sped
2382         # up in the same way, but I didn't have time.
2383         S <- matrix(out$Smat, np + npd, np + npd)[1:np, 1:np]
2384         stderr <- sqrt(diag(solve(S))[np - (ncovs1 - 1):0])
2385     }
2386     return(stderr)
2387 }

TABLE OF CONTENTS

fortranwrappers/fmkstderr [ Functions ]