Error in integrate : evaluation of function gave a result of wrong length after Vectorize

Question

I am trying to simulate survival data from a Weibull PH model with a time varying covariate.

• I wrote a function invS that integrates the hazard function containing the time-varying component.
• Then I used the uniroot function with invS as a function to integrate. For me to do that, I wrote a function invS that returns the integral of the hazard function + log survival probability. In invS, I used integrate function.

I pass this invS function to uniroot with lower and upper bounds defined.

### VALUES TO TEST ####

survival_probability= runif(N,0,1) #Uniform (0,1)

raneff <- mvrnorm(N, mu=c(0,0), Sigma = vcov)



# ORIGINAL SETUP # The covariate in the longitudinal model is not related to the covariates in the survival model

# ORIGINAL SETUP # the baseline covariates in the survival model are binomial whereas the covariate in the longitudinal model is normally distributed



#baseline_covariates_w1 = rbinom(N, 1, 0.2) #can change to other distribution

#baseline_covariates_w2 = rbinom(N, 1, 0.8) #can change to other distribution too

#X2 choose one or the other.

#X2 <- rnorm(N, 0, 1) #original univariate problem



#X2 <- rbinom(N, 1, 0.5) #to reflect randomly assigned treatment

vcov_covariates <-matrix(c(1,0,0.5, 0,1,0.5, 0.5,0.5,1),3,3)

covariate_vector <- mvrnorm(N, mu = c(0,0,0), vcov_covariates )



baseline_covariates_w1 <- covariate_vector[,1]

baseline_covariates_w2 <- covariate_vector[,2]

X2 <- covariate_vector[,3]

alpha = 0.2#correlation between the processes



lambda = 0.2

#beta = c(0.5,3,1)

beta = c(0.5, 3, 1,1)

#beta = c(0.5, 3, 1,1,1,1,1) #this is for the multivariate option

gamma = c(1,1)


getSurvTimeWeibullModel2 <-function(sigma.t, alpha, surv_prob, beta, random_effects, gamma, X2, baseline_cov, N, max.FUtime, method="Weibull-PH"){

#' u/description Get's survival time from Weibull distribution

#' u/param sigma.t scale parameter for the Weibull baseline risk function;

#' u/param alpha association parameter

#' u/param surv_prob double between 0 and 1 N by 1 vector

#' u/param beta vector of fixed effects Should contain beta_0, beta_1, beta_2 and beta_3

#' u/param random_effects # a matrix of N by 2

#' u/param gamma #fixed effects for baseline cov in survival model pby 1

#' u/param X2 # vector of 1 x N or N by 1

#' u/param baseline_cov Nxp

#' u/param N number of individuals

#' u/return double, survival time a vector of N by 1

#'

#compute m_i(t) = ((\beta_0 + b_{i0}) + (\beta_1 + b_{i1})t_{ij} + (\beta_2 X_2i) + \beta_3X_2it_{ij} )

invS <- function (t, u, i) {

timedependentMeasurement <- function(v){

#code up the error free measurement of the current value association of marker (alpha*mi(t))

return(alpha*((beta[1] + random_effects[i, 1]) + (beta[2] + random_effects[i,2])*v + beta[3]*X2[i] + beta[4]*X2[i]*v))

}

h <- function (s) {

TD.i <- timedependentMeasurement(s)

if(method=="weibull-PH"){

exp(log(sigma.t) + (sigma.t - 1) * log(s) +

baseline_cov[,i]%*%(gamma) + TD.i)

}

}

integrate(Vectorize(h), lower=0, upper = t)$value + log(surv_prob)



}

trueTimes <- numeric(N)

for (i in 1:N) {

Root <- uniroot(invS, interval = c(1e-05, max.FUtime), u = surv_prob[i],

i = i)$root

trueTimes[i] <- Root

}

return(trueTimes)

}

getSurvTimeWeibullModel2(sigma.t = 1, alpha, surv_prob=survival_probability, beta, random_effects=raneff, gamma, X2, baseline_cov=cbind(baseline_covariates_w1,baseline_covariates_w2), N=N, max.FUtime=100)

I see in other thread to use Vectorize (Error in integrate : evaluation of function gave a result of wrong length) but I still have the same error. Appreciate if anyone could help!