R/AtCtEt_boot.R

In ACEt: Estimating Dynamic Heritability and Twin Model Comparison

AtCtEt_boot <- function(res, model, data_m, data_d, knot_a, knot_c, knot_e, loc, B=100,alpha=0.05,m=500)
{
spline.main <- sp.spline.estimator(data_m, data_d, model, knot_a, knot_c, knot_e, loc, m)
spline.boots_a <- matrix(NA,m,B)
spline.boots_c <- matrix(NA,m,B)
spline.boots_e <- matrix(NA,m,B)
spline.boots_h <- matrix(NA,m,B)
b_c <- res$beta_c
b_a <- res$beta_a
b_e <- res$beta_e
order <- 3
if(res$n_beta_a==1)
{order <- 1}
B_des_a_m <- splineDesign(res$knots_a, x=data_m[,3], ord=order)	
B_des_a_d <- splineDesign(res$knots_a, x=data_d[,3], ord=order)
order <- 3
if(res$n_beta_c==1)
{order <- 1}
B_des_c_m <- splineDesign(res$knots_c, x=data_m[,3], ord=order)	
B_des_c_d <- splineDesign(res$knots_c, x=data_d[,3], ord=order)
order <- 3
if(res$n_beta_e==1)
{order <- 1}
B_des_e_m <- splineDesign(res$knots_e, x=data_m[,3], ord=order)	
B_des_e_d <- splineDesign(res$knots_e, x=data_d[,3], ord=order)
num_m <- nrow(data_m)
num_d <- nrow(data_d)

for(i in 1:B)
{
pheno_mr <- matrix(0, num_m,2)
for(j in 1:num_m)
{
	var_ct <- sum(B_des_c_m[j,]*b_c)
	var_at <- sum(B_des_a_m[j,]*b_a)
	var_et <- sum(B_des_e_m[j,]*b_e)
	#if(res$n_beta_a>1)
	#{var_at <- exp(var_at)}
	var_at <- exp(var_at)
	#if(res$n_beta_c>1)
	#{var_ct <- exp(var_ct)}
	var_ct <- exp(var_ct)
	#if(res$n_beta_e>1)
	#{var_et <- exp(var_et)}
	var_et <- exp(var_et)
	# var_ct <- var_c
	sigma <- matrix(c(var_at+var_ct+var_et,var_at+var_ct,var_at+var_ct,var_at+var_ct+var_et),2,2)

	pheno_mr[j,1:2] <- mvrnorm(1, rep(0,2), sigma)
}

pheno_dr <- matrix(0, num_d,2)	
		
for(j in 1:num_d)
{
	var_ct <- sum(B_des_c_d[j,]*b_c)
	var_at <- sum(B_des_a_d[j,]*b_a)
	var_et <- sum(B_des_e_d[j,]*b_e)
	#if(res$n_beta_a>1)
	#{var_at <- exp(var_at)}
	var_at <- exp(var_at)
	#if(res$n_beta_c>1)
	#{var_ct <- exp(var_ct)}
	var_ct <- exp(var_ct)
	#if(res$n_beta_e>1)
	#{var_et <- exp(var_et)}
	var_et <- exp(var_et)
	sigma <- matrix(c(var_at+var_ct+var_et,0.5*var_at+var_ct,0.5*var_at+var_ct,var_at+var_ct+var_et),2,2)

	pheno_dr[j,1:2] <- mvrnorm(1, rep(0,2), sigma)
}


spline.boots <- sp.spline.estimator(cbind(pheno_mr,data_m[,3]),cbind(pheno_dr, data_d[,3]),model=model, knot_a, knot_c, knot_e, loc, m=m)
spline.boots_a[,i] <- spline.boots$est_a
spline.boots_c[,i] <- spline.boots$est_c
spline.boots_e[,i] <- spline.boots$est_e
spline.boots_h[,i] <- spline.boots$est_a/(spline.boots$est_a+spline.boots$est_c+spline.boots$est_e)
}
# Result has m rows and B columns
#cis.lower_a <- 2*spline.main$est_a - apply(spline.boots_a,1,quantile,probs=1-alpha/2)
#cis.upper_a <- 2*spline.main$est_a - apply(spline.boots_a,1,quantile,probs=alpha/2)
#cis.lower_c <- 2*spline.main$est_c - apply(spline.boots_c,1,quantile,probs=1-alpha/2)
#cis.upper_c <- 2*spline.main$est_c - apply(spline.boots_c,1,quantile,probs=alpha/2)
#cis.lower_e <- 2*spline.main$est_e - apply(spline.boots_e,1,quantile,probs=1-alpha/2)
#cis.upper_e <- 2*spline.main$est_e - apply(spline.boots_e,1,quantile,probs=alpha/2)

# percentile method
cis.upper_a <- apply(spline.boots_a,1,quantile,probs=1-alpha/2)
cis.lower_a <- apply(spline.boots_a,1,quantile,probs=alpha/2)
cis.upper_c <- apply(spline.boots_c,1,quantile,probs=1-alpha/2)
cis.lower_c <- apply(spline.boots_c,1,quantile,probs=alpha/2)
cis.upper_e <- apply(spline.boots_e,1,quantile,probs=1-alpha/2)
cis.lower_e <- apply(spline.boots_e,1,quantile,probs=alpha/2)
cis.upper_h <- apply(spline.boots_h,1,quantile,probs=1-alpha/2)
cis.lower_h <- apply(spline.boots_h,1,quantile,probs=alpha/2)

# return(list(lower.ci_a=cis.lower_a,upper.ci_a=cis.upper_a,lower.ci_c=cis.lower_c,upper.ci_c=cis.upper_c, lower.ci_e=cis.lower_e,upper.ci_e=cis.upper_e, x=seq(from=res$min_t, to=res$max_t, length.out=m),boots_a=spline.boots_a, boots_c=spline.boots_c, boots_e=spline.boots_e))
return(list(lower.ci_a=cis.lower_a,upper.ci_a=cis.upper_a,lower.ci_c=cis.lower_c,upper.ci_c=cis.upper_c, lower.ci_e=cis.lower_e,upper.ci_e=cis.upper_e, lower.ci_h=cis.lower_h, upper.ci_h=cis.upper_h, x=seq(from=res$min_t, to=res$max_t, length.out=m)))

}


sp.spline.estimator <- function(data_m, data_d, model, knot_a, knot_c, knot_e, loc, m) {
# Fit spline to data, with cross-validation to pick lambda
fit <- AtCtEt(data_m, data_d, model, knot_a, knot_c, knot_e, loc)
T_m <- rep(data_m[,3], each=2)
T_d <- rep(data_d[,3], each=2)
eval.grid <- seq(from=min(T_m, T_d), to=max(T_m, T_d), length.out=m)
order <- 3
if(fit$n_beta_a==1)
{order <- 1}
bb_a <- splineDesign(fit$knots_a, x = eval.grid, ord=order, outer.ok = TRUE)
order <- 3
if(fit$n_beta_c==1)
{order <- 1}
bb_c <- splineDesign(fit$knots_c, x = eval.grid, ord=order, outer.ok = TRUE)
order <- 3
if(fit$n_beta_e==1)
{order <- 1}
bb_e <- splineDesign(fit$knots_e, x = eval.grid, ord=order, outer.ok = TRUE)

est_a <- bb_a%*%fit$beta_a
#if(fit$n_beta_a>1)
#{est_a <- exp(est_a)}
est_a <- exp(est_a)

est_c <- bb_c%*%fit$beta_c
#if(fit$n_beta_c>1)
#{est_c <- exp(est_c)}
est_c <- exp(est_c)

est_e <- bb_e%*%fit$beta_e
#if(fit$n_beta_e>1)
#{est_e <- exp(est_e)}
est_e <- exp(est_e)

return(list(e = fit$beta_e, c = fit$beta_c, a = fit$beta_a, est_c = est_c,est_a = est_a, est_e = est_e))

}

AtDtEt_boot <- function(res, model, data_m, data_d, knot_a, knot_d, knot_e, loc, B=100,alpha=0.05,m=500)
{
  spline.main <- sp.spline.estimator_ade(data_m, data_d, model, knot_a, knot_d, knot_e, loc, m)
  spline.boots_a <- matrix(NA,m,B)
  spline.boots_d <- matrix(NA,m,B)
  spline.boots_e <- matrix(NA,m,B)
  spline.boots_h <- matrix(NA,m,B)
  b_d <- res$beta_d
  b_a <- res$beta_a
  b_e <- res$beta_e
  order <- 3
  if(res$n_beta_a==1)
  {order <- 1}
  B_des_a_m <- splineDesign(res$knots_a, x=data_m[,3], ord=order)	
  B_des_a_d <- splineDesign(res$knots_a, x=data_d[,3], ord=order)
  order <- 3
  if(res$n_beta_d==1)
  {order <- 1}
  B_des_d_m <- splineDesign(res$knots_d, x=data_m[,3], ord=order)	
  B_des_d_d <- splineDesign(res$knots_d, x=data_d[,3], ord=order)
  order <- 3
  if(res$n_beta_e==1)
  {order <- 1}
  B_des_e_m <- splineDesign(res$knots_e, x=data_m[,3], ord=order)	
  B_des_e_d <- splineDesign(res$knots_e, x=data_d[,3], ord=order)
  num_m <- nrow(data_m)
  num_d <- nrow(data_d)
  
  for(i in 1:B)
  {
    pheno_mr <- matrix(0, num_m,2)
    for(j in 1:num_m)
    {
      var_dt <- sum(B_des_d_m[j,]*b_d)
      var_at <- sum(B_des_a_m[j,]*b_a)
      var_et <- sum(B_des_e_m[j,]*b_e)
      #if(res$n_beta_a>1)
      #{var_at <- exp(var_at)}
      var_at <- exp(var_at)
      #if(res$n_beta_c>1)
      #{var_ct <- exp(var_ct)}
      var_dt <- exp(var_dt)
      #if(res$n_beta_e>1)
      #{var_et <- exp(var_et)}
      var_et <- exp(var_et)
      # var_ct <- var_c
      sigma <- matrix(c(var_at+var_dt+var_et,var_at+var_dt,var_at+var_dt,var_at+var_dt+var_et),2,2)
      
      pheno_mr[j,1:2] <- mvrnorm(1, rep(0,2), sigma)
    }
    
    pheno_dr <- matrix(0, num_d,2)	
    
    for(j in 1:num_d)
    {
      var_dt <- sum(B_des_d_d[j,]*b_d)
      var_at <- sum(B_des_a_d[j,]*b_a)
      var_et <- sum(B_des_e_d[j,]*b_e)
      #if(res$n_beta_a>1)
      #{var_at <- exp(var_at)}
      var_at <- exp(var_at)
      #if(res$n_beta_c>1)
      #{var_ct <- exp(var_ct)}
      var_dt <- exp(var_dt)
      #if(res$n_beta_e>1)
      #{var_et <- exp(var_et)}
      var_et <- exp(var_et)
      sigma <- matrix(c(var_at+var_dt+var_et,0.5*var_at+0.25*var_dt,0.5*var_at+0.25*var_dt,var_at+var_dt+var_et),2,2)
      
      pheno_dr[j,1:2] <- mvrnorm(1, rep(0,2), sigma)
    }
    
    
    spline.boots <- sp.spline.estimator_ade(cbind(pheno_mr,data_m[,3]),cbind(pheno_dr, data_d[,3]),model=model, knot_a, knot_d, knot_e, loc, m=m)
    spline.boots_a[,i] <- spline.boots$est_a
    spline.boots_d[,i] <- spline.boots$est_d
    spline.boots_e[,i] <- spline.boots$est_e
    spline.boots_h[,i] <- (spline.boots$est_a+spline.boots$est_d)/(spline.boots$est_a+spline.boots$est_d+spline.boots$est_e)
  }
  
  # percentile method
  cis.upper_a <- apply(spline.boots_a,1,quantile,probs=1-alpha/2)
  cis.lower_a <- apply(spline.boots_a,1,quantile,probs=alpha/2)
  cis.upper_d <- apply(spline.boots_d,1,quantile,probs=1-alpha/2)
  cis.lower_d <- apply(spline.boots_d,1,quantile,probs=alpha/2)
  cis.upper_e <- apply(spline.boots_e,1,quantile,probs=1-alpha/2)
  cis.lower_e <- apply(spline.boots_e,1,quantile,probs=alpha/2)
  cis.upper_h <- apply(spline.boots_h,1,quantile,probs=1-alpha/2)
  cis.lower_h <- apply(spline.boots_h,1,quantile,probs=alpha/2)
  
  # return(list(lower.ci_a=cis.lower_a,upper.ci_a=cis.upper_a,lower.ci_c=cis.lower_c,upper.ci_c=cis.upper_c, lower.ci_e=cis.lower_e,upper.ci_e=cis.upper_e, x=seq(from=res$min_t, to=res$max_t, length.out=m),boots_a=spline.boots_a, boots_c=spline.boots_c, boots_e=spline.boots_e))
  return(list(lower.ci_a=cis.lower_a,upper.ci_a=cis.upper_a,lower.ci_d=cis.lower_d,upper.ci_d=cis.upper_d, lower.ci_e=cis.lower_e,upper.ci_e=cis.upper_e, lower.ci_h=cis.lower_h, upper.ci_h=cis.upper_h, x=seq(from=res$min_t, to=res$max_t, length.out=m)))
  
}


sp.spline.estimator_ade <- function(data_m, data_d, model, knot_a, knot_d, knot_e, loc, m) {
  # Fit spline to data, with cross-validation to pick lambda
  fit <- AtDtEt(data_m, data_d, model, knot_a, knot_d, knot_e, loc)
  T_m <- rep(data_m[,3], each=2)
  T_d <- rep(data_d[,3], each=2)
  eval.grid <- seq(from=min(T_m, T_d), to=max(T_m, T_d), length.out=m)
  order <- 3
  if(fit$n_beta_a==1)
  {order <- 1}
  bb_a <- splineDesign(fit$knots_a, x = eval.grid, ord=order, outer.ok = TRUE)
  order <- 3
  if(fit$n_beta_d==1)
  {order <- 1}
  bb_d <- splineDesign(fit$knots_d, x = eval.grid, ord=order, outer.ok = TRUE)
  order <- 3
  if(fit$n_beta_e==1)
  {order <- 1}
  bb_e <- splineDesign(fit$knots_e, x = eval.grid, ord=order, outer.ok = TRUE)
  
  est_a <- bb_a%*%fit$beta_a
  #if(fit$n_beta_a>1)
  #{est_a <- exp(est_a)}
  est_a <- exp(est_a)
  
  est_d <- bb_d%*%fit$beta_d
  #if(fit$n_beta_c>1)
  #{est_c <- exp(est_c)}
  est_d <- exp(est_d)
  
  est_e <- bb_e%*%fit$beta_e
  #if(fit$n_beta_e>1)
  #{est_e <- exp(est_e)}
  est_e <- exp(est_e)
  
  return(list(e = fit$beta_e, d = fit$beta_d, a = fit$beta_a, est_d = est_d,est_a = est_a, est_e = est_e))
  
}