Nothing
###############################################################################
# Functions for 1st step
###############################################################################
fitClayton <- function(data,
varcor = FALSE,
optimx.method,
MAXFUN = 1e8,
ini_rho = 0.5,
# ini_kTau = 0.5,
twoStage = TRUE) {
OPTIMX_CONTROL <- list(
maxit = MAXFUN,
dowarn = FALSE,
starttests = FALSE,
kkt = FALSE
)
# Initial values of parameters: marginal estiamtes
Mpars <- margFits(data)
INIpars <- c(
log(copula::iRho(claytonCopula(), ini_rho)),
#log(copula::iTau(claytonCopula(), ini_kTau)),
log(Mpars$lambdaS),
log(Mpars$rhoS),
Mpars$alpha,
log(Mpars$lambdaT),
log(Mpars$rhoT),
Mpars$beta
)
nTrials <- (length(INIpars) - 1) / 6
# Possibly fix baseline hazard parameters to their marginal values
if (twoStage) {
fixed <- c(
NA,
# copula dependence parameter
log(Mpars$lambdaS),
# baseline scales for S
log(Mpars$rhoS),
# baseline shapes for S
rep(NA, nTrials),
# alphas
log(Mpars$lambdaT),
# baseline scales for T
log(Mpars$rhoT),
# baseline shapes for T
rep(NA, nTrials)
) # betas
} else {
fixed <- rep(NA, length(INIpars))
}
# library('optimx')
RES <- optimx(
INIpars,
mloglik,
FIXpars = fixed,
data = data,
family = 'clayton',
method = optimx.method,
control = OPTIMX_CONTROL
)
if (twoStage) {
RES[, which(!is.na(fixed))] <- fixed[!is.na(fixed)]
}
if (varcor) {
# library('optextras')
attr(RES, 'grad') <-
grnd(as.numeric(RES[, 1:attr(RES, 'npar')]),
mloglik,
data = data,
family = 'clayton')
attr(RES, 'hessian') <- optimHess(
RES[, 1:length(INIpars)],
mloglik,
data = data,
family = 'clayton',
control = OPTIMX_CONTROL
)
# library('msm')
tranform <- eval(parse(
text = paste(
# theta
'list(~ exp(x1),\n',
# lambdaS and rhoS
paste0('~ exp(x', 1 + 1:(2 * nTrials), '),', collapse = '\n'),
# alpha
paste0('~ x', 1 + 2 * nTrials + 1:nTrials, ',', collapse = '\n'),
# lambdaT and rhoT
paste0('~ exp(x', 1 + 3 * nTrials + 1:(2 * nTrials), '),', collapse = '\n'),
# beta
paste0('~ x', 1 + 5 * nTrials + 1:nTrials, collapse = ',\n'),
')'
)
))
VarCor <- try(deltamethod(
tranform,
mean = as.numeric(RES[1:length(INIpars)]),
cov = solve(attr(RES, 'hessian')),
ses = FALSE
))
}
RES <- list(
theta = as.numeric(exp(RES[, 1])),
lambdaS = as.numeric(exp(RES[, 1 + 1:nTrials])),
rhoS = as.numeric(exp(RES[, 1 + nTrials + 1:nTrials])),
alpha = as.numeric(RES[, 1 + 2 * nTrials + 1:nTrials]),
lambdaT = as.numeric(exp(RES[, 1 + 3 * nTrials + 1:nTrials])),
rhoT = as.numeric(exp(RES[, 1 + 4 * nTrials + 1:nTrials])),
beta = as.numeric(RES[, 1 + 5 * nTrials + 1:nTrials]),
optimxRES = RES
)
if (varcor)
RES <- c(RES, list(VarCor = VarCor))
return(RES)
}
kTau.clay <- function(x) {
# library('copula')
copula::tau(claytonCopula(x$theta))
}
###############################################################################
fitPlackett <- function(data,
varcor = FALSE,
optimx.method,
MAXFUN = 1e8,
ini_rho = 0.5,
# ini_kTau = 0.5,
twoStage = TRUE) {
OPTIMX_CONTROL <- list(
maxit = MAXFUN,
dowarn = FALSE,
starttests = FALSE,
kkt = FALSE
)
# Initial values of parameters: marginal estimates
Mpars <- margFits(data)
INIpars <- c(
log(copula::iRho(plackettCopula(), ini_rho)),
#log(copula::iTau(plackettCopula(), ini_kTau)),
log(Mpars$lambdaS),
log(Mpars$rhoS),
Mpars$alpha,
log(Mpars$lambdaT),
log(Mpars$rhoT),
Mpars$beta
)
nTrials <- (length(INIpars) - 1) / 6
# Possibly fix baseline hazard parameters to their marginal values
if (twoStage) {
fixed <- c(
NA,
# copula dependence parameter
log(Mpars$lambdaS),
# baseline scales for S
log(Mpars$rhoS),
# baseline shapes for S
rep(NA, nTrials),
# alphas
log(Mpars$lambdaT),
# baseline scales for T
log(Mpars$rhoT),
# baseline shapes for T
rep(NA, nTrials)
) # betas
} else {
fixed <- rep(NA, length(INIpars))
}
# library('optimx')
RES <- optimx(
INIpars,
mloglik,
FIXpars = fixed,
data = data,
family = 'plackett',
method = optimx.method,
control = OPTIMX_CONTROL
)
if (twoStage) {
RES[, which(!is.na(fixed))] <- fixed[!is.na(fixed)]
}
if (varcor) {
# library('optextras')
attr(RES, 'grad') <-
grnd(as.numeric(RES[, 1:attr(RES, 'npar')]),
mloglik,
data = data,
family = 'plackett')
attr(RES, 'hessian') <- optimHess(
RES[1:length(INIpars)],
mloglik,
data = data,
family = 'plackett',
control = OPTIMX_CONTROL
)
# library('msm')
tranform <- eval(parse(
text = paste(
# theta
'list(~ exp(x1),\n',
# lambdaS and rhoS
paste0('~ exp(x', 1 + 1:(2 * nTrials), '),', collapse = '\n'),
# alpha
paste0('~ x', 1 + 2 * nTrials + 1:nTrials, ',', collapse = '\n'),
# lambdaT and rhoT
paste0('~ exp(x', 1 + 3 * nTrials + 1:(2 * nTrials), '),', collapse = '\n'),
# beta
paste0('~ x', 1 + 5 * nTrials + 1:nTrials, collapse = ',\n'),
')'
)
))
VarCor <- try(deltamethod(
tranform,
mean = as.numeric(RES[1:length(INIpars)]),
cov = solve(attr(RES, 'hessian')),
ses = FALSE
))
}
RES <- list(
theta = as.numeric(exp(RES[, 1])),
lambdaS = as.numeric(exp(RES[, 1 + 1:nTrials])),
rhoS = as.numeric(exp(RES[, 1 + nTrials + 1:nTrials])),
alpha = as.numeric(RES[, 1 + 2 * nTrials + 1:nTrials]),
lambdaT = as.numeric(exp(RES[, 1 + 3 * nTrials + 1:nTrials])),
rhoT = as.numeric(exp(RES[, 1 + 4 * nTrials + 1:nTrials])),
beta = as.numeric(RES[, 1 + 5 * nTrials + 1:nTrials]),
optimxRES = RES
)
if (varcor)
RES <- c(RES, list(VarCor = VarCor))
return(RES)
}
kTau.plack <- function(x) {
# library('copula')
copula::tau(plackettCopula(x$theta))
}
###############################################################################
fitHougaard <- function(data,
varcor = FALSE,
optimx.method,
MAXFUN = 1e8,
ini_rho = 0.5,
# ini_kTau = 0.5,
twoStage = TRUE) {
OPTIMX_CONTROL <- list(
maxit = MAXFUN,
dowarn = FALSE,
starttests = FALSE,
kkt = FALSE
)
# Initial values of parameters: marginal estiamtes
Mpars <- margFits(data)
INIpars <- c(
log(1 / copula::iRho(gumbelCopula(), ini_rho)),
#log(1 - ini_kTau),
log(Mpars$lambdaS),
log(Mpars$rhoS),
Mpars$alpha,
log(Mpars$lambdaT),
log(Mpars$rhoT),
Mpars$beta
)
nTrials <- (length(INIpars) - 1) / 6
# Possibly fix baseline hazard parameters to their marginal values
if (twoStage) {
fixed <- c(
NA,
# copula dependence parameter
log(Mpars$lambdaS),
# baseline scales for S
log(Mpars$rhoS),
# baseline shapes for S
rep(NA, nTrials),
# alphas
log(Mpars$lambdaT),
# baseline scales for T
log(Mpars$rhoT),
# baseline shapes for T
rep(NA, nTrials)
) # betas
} else {
fixed <- rep(NA, length(INIpars))
}
# library('optimx')
RES <- optimx(
INIpars,
mloglik,
FIXpars = fixed,
data = data,
family = 'hougaard',
method = optimx.method,
control = OPTIMX_CONTROL
)
if (twoStage) {
RES[, which(!is.na(fixed))] <- fixed[!is.na(fixed)]
}
if (varcor) {
# library('optimx')
attr(RES, 'grad') <-
grnd(as.numeric(RES[, 1:attr(RES, 'npar')]),
mloglik,
data = data,
family = 'hougaard')
attr(RES, 'hessian') <- optimHess(
RES[1:length(INIpars)],
mloglik,
data = data,
family = 'hougaard',
control = OPTIMX_CONTROL
)
# library('msm')
tranform <- eval(parse(
text = paste(
# theta
'list(~ exp(-exp(x1)),\n',
# lambdaS and rhoS
paste0('~ exp(x', 1 + 1:(2 * nTrials), '),', collapse = '\n'),
# alpha
paste0('~ x', 1 + 2 * nTrials + 1:nTrials, ',', collapse = '\n'),
# lambdaT and rhoT
paste0('~ exp(x', 1 + 3 * nTrials + 1:(2 * nTrials), '),', collapse = '\n'),
# beta
paste0('~ x', 1 + 5 * nTrials + 1:nTrials, collapse = ',\n'),
')'
)
))
VarCor <- try(deltamethod(
tranform,
mean = as.numeric(RES[1:length(INIpars)]),
cov = solve(attr(RES, 'hessian')),
ses = FALSE
))
}
RES <- list(
theta = as.numeric(exp(-exp(RES[, 1]))),
lambdaS = as.numeric(exp(RES[, 1 + 1:nTrials])),
rhoS = as.numeric(exp(RES[, 1 + nTrials + 1:nTrials])),
alpha = as.numeric(RES[, 1 + 2 * nTrials + 1:nTrials]),
lambdaT = as.numeric(exp(RES[, 1 + 3 * nTrials + 1:nTrials])),
rhoT = as.numeric(exp(RES[, 1 + 4 * nTrials + 1:nTrials])),
beta = as.numeric(RES[, 1 + 5 * nTrials + 1:nTrials]),
optimxRES = RES
)
if (varcor)
RES <- c(RES, list(VarCor = VarCor))
return(RES)
}
kTau.houg <- function(x) {
1 - x$theta
}
###############################################################################
margFits <- function(data) {
tS <- data$timeS
tT <- data$timeT
cS <- data$statusS
cT <- data$statusT
# library('eha')
# Surrogate S
sreg <- survreg(Surv(timeS, statusS) ~
trt:trialref + strata(trialref) + trialref,
data = data)
lS <- exp(-(sreg$coeff['(Intercept)'] +
c(0, sreg$coeff[(names(sreg$coeff) != '(Intercept)') &
!grepl('trt', names(sreg$coeff))])) / sreg$scale)
rS <- 1 / sreg$scale
bS <- -sreg$coeff[grep('trt', names(sreg$coef))] / sreg$scale
# sS <- exp(-lS[data$trialref] * tS^(rS[data$trialref]) *
# exp(lS[data$trialref] * data$trt))
# l_fS <- log(lS[data$trialref] * exp(lS[data$trialref] * data$trt) *
# rS[data$trialref] *
# tS^(rS[data$trialref]-1) *
# sS)
# True endpoint T
treg <- survreg(Surv(timeT, statusT) ~
trt:trialref + strata(trialref) + trialref,
data = data)
lT <- exp(-(treg$coeff['(Intercept)'] +
c(0, treg$coeff[(names(treg$coeff) != '(Intercept)') &
!grepl('trt', names(treg$coeff))])) / treg$scale)
rT <- 1 / treg$scale
bT <- -treg$coeff[grep('trt', names(treg$coef))] / treg$scale
# sT <- exp(-lT[data$trialref] * tT^(rT[data$trialref]) *
# exp(lT[data$trialref] * data$trt))
# l_fT <- log(lT[data$trialref] * exp(lT[data$trialref] * data$trt) *
# rT[data$trialref] *
# tT^(rT[data$trialref]-1) *
# sT)
RES <- list(
lambdaS = as.numeric(lS),
rhoS = as.numeric(rS),
alpha = as.numeric(bS),
lambdaT = as.numeric(lT),
rhoT = as.numeric(rT),
beta = as.numeric(bT)
)
return(RES)
}
mloglik <- function(pars,
FIXpars = rep(NA, length(pars)),
data,
family = c('clayton', 'plackett', 'hougaard')) {
family <- match.arg(family)
# If any parameter have a fixed value, use this latter
pars[!is.na(FIXpars)] <- FIXpars[!is.na(FIXpars)]
# pars contains
# - one value for log(theta)
# - nTrials values for log(lambdaS)
# - nTrials values for log(rhoS)
# - nTrials values for alphaS
# - nTrials values for log(lambdaT)
# - nTrials values for log(rhoT)
# - nTrials values for betaT
theta <- exp(pars[1])
nTrials <- (length(pars) - 1) / 6
lambdaS <- exp(pars[1 + 1:nTrials])
rhoS <- exp(pars[1 + nTrials + 1:nTrials])
alpha <- pars[1 + 2 * nTrials + 1:nTrials]
lambdaT <- exp(pars[1 + 3 * nTrials + 1:nTrials])
rhoT <- exp(pars[1 + 4 * nTrials + 1:nTrials])
beta <- pars[1 + 5 * nTrials + 1:nTrials]
ll <- 0
tS <- data$timeS
tT <- data$timeT
cS <- data$statusS
cT <- data$statusT
lS <-
lambdaS[data$trialref] * exp(alpha[data$trialref] * data$trt)
rS <- rhoS[data$trialref]
sS <- exp(-lS * tS ^ rS)
l_fS <- log(lS * rS * tS ^ (rS - 1) * sS)
lT <- lambdaT[data$trialref] * exp(beta[data$trialref] * data$trt)
rT <- rhoT[data$trialref]
sT <- exp(-lT * tT ^ rT)
l_fT <- log(lT * rT * tT ^ (rT - 1) * sT)
if (family == 'clayton') {
### S_ST(s, t) = (S_S(s)^(-theta) + S_T(t)^(-theta) - 1)^(-1/theta)
l_Sst <- -log(sS ^ (-theta) + sT ^ (-theta) - 1) / theta
# Both censored
sel <- !cS & !cT
if (any(sel)) {
ll <- ll + sum(l_Sst[sel])
}
# S censored only
sel <- !cS & cT
if (any(sel)) {
ll <- ll + sum(l_fT[sel]) +
(1 + theta) * sum(l_Sst[sel] - log(sT[sel]))
}
# T censored only
sel <- cS & !cT
if (any(sel)) {
ll <- ll + sum(l_fS[sel]) +
(1 + theta) * sum(l_Sst[sel] - log(sS[sel]))
}
# both measured
sel <- cS & cT
if (any(sel)) {
ll <- ll + sum(sel) * log(1 + theta) +
sum(l_fS[sel] + l_fT[sel]) +
(1 + 2 * theta) * sum(l_Sst[sel]) +-(1 + theta) * sum(log(sS[sel] * sT[sel]))
}
} else if (family == 'plackett') {
if (theta == 1) {
l_Sst <- log(sS) + log(sT)
if (any(!cS & !cT)) {
ll <- ll + sum(l_Sst[!cS & !cT])
}
if (any(!cS & cT)) {
ll <- ll + sum(log(sT[!cS & cT]) + l_fS[!cS & cT])
}
if (any(cS & !cT)) {
ll <- ll + sum(log(sS[cS & !cT]) + l_fT[cS & !cT])
}
if (any(cS & cT)) {
ll <- ll + sum(l_fS[cS & cT] + l_fT[cS & cT])
}
} else {
### S_ST(s, t) = (Q - sqrt(R))/(2 *(theta - 1))
### Q = 1 + (theta - 1) * (S_S(s)+S_T(t))
### R = Q^2 - 4 * S_S(s) * S_T(t) * theta * (theta - 1)
Q <- 1 + (theta - 1) * (sS + sT)
R <- Q ^ 2 - 4 * theta * (theta - 1) * sS * sT
l_Sst <- log((Q - R ^ (1 / 2)) / (2 * (theta - 1)))
# Both censored
sel <- !cS & !cT
if (any(sel)) {
ll <- ll + sum(l_Sst[sel])
}
# S censored only
sel <- !cS & cT
if (any(sel)) {
ll <- ll - sum(sel) * log(2) +
sum(log(1 - (Q[sel] - 2 * theta * sS[sel]) * R[sel] ^ (-1 / 2)) + l_fT[sel])
}
# T censored only
sel <- cS & !cT
if (any(sel)) {
ll <- ll - sum(sel) * log(2) +
sum(log(1 - (Q[sel] - 2 * theta * sT[sel]) * R[sel] ^ (-1 / 2)) + l_fS[sel])
}
# both measured
sel <- cS & cT
if (any(sel)) {
ll <- ll + sum(sel) * log(theta) + sum(-3 / 2 * log(R[sel]) +
log(Q[sel] - 2 * (theta - 1) * sS[sel] * sT[sel]) + l_fS[sel] + l_fT[sel])
}
}
} else if (family == 'hougaard') {
### S_ST(s, t) = exp(-Q^theta)
### Q = (-log u)^1/theta + (-log v)^1/theta
Q = (-log(sS)) ^ (1 / theta) + (-log(sT)) ^ (1 / theta)
l_Sst <- -Q ^ theta
# Both censored
sel <- !cS & !cT
if (any(sel)) {
ll <- ll + sum(l_Sst[sel])
}
# S censored only
sel <- !cS & cT
if (any(sel)) {
ll <- ll + sum(l_fT[sel]) +
(1 / theta - 1) * sum(log(-log(sT[sel]))) +-sum(log(sT[sel])) + sum(l_Sst[sel]) +
(theta - 1) * sum(log(Q[sel]))
}
# T censored only
sel <- cS & !cT
if (any(sel)) {
ll <- ll + sum(l_fS[sel]) +
(1 / theta - 1) * sum(log(-log(sS[sel]))) +-sum(log(sS[sel])) + sum(l_Sst[sel]) +
(theta - 1) * sum(log(Q[sel]))
}
# both measured
sel <- cS & cT
if (any(sel)) {
ll <- ll + sum(l_fS[sel] + l_fT[sel]) +
(1 / theta - 1) * sum(log(-log(sS[sel])) + log(-log(sT[sel]))) +-sum(log(sS[sel]) + log(sT[sel])) + sum(l_Sst[sel]) +
(theta - 2) * sum(log(Q[sel])) +
sum(log(1 / theta - 1 + Q[sel] ^ theta))
}
}
return(-ll)
}
###############################################################################
# Functions for 2nd step
###############################################################################
copuSurr <-
function(data,
family = c('clayton', 'plackett', 'hougaard'),
varcor1 = FALSE,
optimx.method,
INIrho = 0.5,
# INIkTau = 0.5,
twoStage = TRUE) {
startTime <- Sys.time()
family <- tolower(family)
family <- match.arg(family)
margPars <- margFits(data)
if (family == 'clayton') {
step1 <- fitClayton(
data,
varcor = varcor1,
optimx.method = optimx.method,
ini_rho = INIrho,
#ini_kTau = INIkTau
twoStage = twoStage
)
kTau <- kTau.clay(step1)
} else if (family == 'plackett') {
step1 <- fitPlackett(
data,
varcor = varcor1,
optimx.method = optimx.method,
ini_rho = INIrho,
#ini_kTau = INIkTau
twoStage = twoStage
)
kTau <- kTau.plack(step1)
} else if (family == 'hougaard') {
step1 <- fitHougaard(
data,
varcor = varcor1,
optimx.method = optimx.method,
ini_rho = INIrho,
#ini_kTau = INIkTau
twoStage = twoStage
)
kTau <- kTau.houg(step1)
}
step1est <- data.frame(
trteff = c(step1$alpha, step1$beta),
ep = rep(c('S', 'T'), each = length(step1$alpha)),
trial = rep(1:length(step1$alpha), 2),
weights = as.vector(table(data$trialref))
)
step2 <- lm(beta ~ alpha,
data = as.data.frame(step1[c('alpha', 'beta')]))
alpha <- mean(step1$alpha)
beta <- mean(step1$beta)
R2 <- summary(step2)$r.squared
RES <- list(
unadj = list(
kTau = kTau,
alpha = alpha,
beta = beta,
step1 = step1,
step2 = step2,
R2 = R2,
ranef = NULL,
optimxRES = step1$optimxRES,
VarCor2 = NULL
)
)
if (varcor1)
RES$unadj <- c(RES$unadj, list(VarCor1 = step1$VarCor))
RES$unadj$runTime <- Sys.time() - startTime
############################################################################
# 2nd step using the mvmeta package #
# #
# Gasparrini A, Armstrong B, Kenward MG. Multivariate meta-analysis for #
# non-linear abd other multi-parameter associations. Stat Med 2012; 31. #
# #
# the same as performed in SAS by #
# Burzykowski T,Molenberghs G, Buyse M, Geys H, Renard D. Validation of #
# surrogate and points in multiple randomized clinical trials with fail- #
# ure time end points. Appl Staist 2001; 50. #
# and #
# van Houwelingen mstrong B, Kenward MG. Multivariate meta-analysis for #
# non-linear abd other multi-parameter associations. Stat Med 2012; 31. #
############################################################################
poss <- cumsum(sapply(step1, length))
ind_AlphaBeta <- which(names(poss) %in% c('alpha', 'beta'))
ind_AlphaBeta <- mapply(':',
poss[ind_AlphaBeta - 1] + 1,
poss[ind_AlphaBeta])
myS <- with(step1est,
cbind(
alpha_Var = pmax(diag(step1$VarCor[ind_AlphaBeta[, 1],
ind_AlphaBeta[, 1]]), 1e-2),
albe_Covar = diag(step1$VarCor[ind_AlphaBeta[, 1]
, ind_AlphaBeta[, 2]]),
beta_Var = pmax(diag(step1$VarCor[ind_AlphaBeta[, 2],
ind_AlphaBeta[, 2]]), 1e-2)
))
# library('mvmeta')
bvmodel <- mvmeta(
cbind(alpha, beta),
S = myS,
data = as.data.frame(step1[c('alpha', 'beta')]),
control = list(hessian = TRUE)
)
alpha <- bvmodel$coefficients[1]
beta <- bvmodel$coefficients[2]
R2 <- bvmodel$Psi[1, 2] ^ 2 / prod(diag(bvmodel$Psi))
# -------------------------------------------------------------------------- #
# Regression parameters ####
PARS <- c(
alpha = bvmodel$coefficients[, 'alpha'],
beta = bvmodel$coefficients[, 'beta'],
Psi.aa = bvmodel$Psi['alpha', 'alpha'],
Psi.ab = bvmodel$Psi['alpha', 'beta'],
Psi.bb = bvmodel$Psi['beta', 'beta']
)
# library(msm)
bvmodel$Psi.vcov <- deltamethod(
g = list( ~ x1 ^ 2, ~ x1 * x2, ~ x2 ^ 2 + x3 ^ 2),
# for this tranformation, see the Remark on page 612 of
# van Houwelingen, Arends, Stijnen, Stat Med 2002; 21:589-624
mean = bvmodel$par,
cov = solve(-bvmodel$hessian),
ses = FALSE
)
VCOV <- diag(5)
VCOV[1:2, 1:2] <- bvmodel$vcov
VCOV[3:5, 3:5] <- bvmodel$Psi.vcov
rownames(VCOV) <- colnames(VCOV) <- names(PARS)
# The following parameter estimates are given by equations (18.23), page 331
# Burzykowski and Buyse, An alternative measure for surrogate
# endpoint validation, in Burzykowski, Molenberghs, Buyse, 2005
# The evaluation of surrogate endpoints, New York Springer
bvmodel$gamma <-
c(gamma0 = as.numeric(PARS['beta'] - PARS['alpha'] *
PARS['Psi.ab'] / PARS['Psi.aa']),
gamma1 = as.numeric(PARS['Psi.ab'] / PARS['Psi.aa']))
bvmodel$gamma.vcov <- deltamethod(
g = list( ~ x2 - x1 * x4 / x3, ~ x4 / x3),
mean = PARS,
cov = VCOV,
ses = FALSE
)
bvmodel$sigma <- PARS['Psi.bb'] * (1 - R2)
# -------------------------------------------------------------------------- #
runTime.Adj <- Sys.time() - startTime
RES$adj <- list(
kTau = kTau,
alpha = alpha,
beta = beta,
step1 = step1,
step2 = bvmodel,
R2 = R2,
ranef = blup(bvmodel),
optimxRES = step1$optimxRES,
VarCor2 = bvmodel$Psi
)
if (varcor1)
RES$adj <- c(RES$adj, list(VarCor1 = step1$VarCor))
RES$adj$runTime <- Sys.time() - startTime
return(RES)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.