R/AtCtEtp.R

Defines functions AtCtEtp

Documented in AtCtEtp

AtCtEtp <-
function(data_m, data_d, knot_a=8, knot_c=8, knot_e=8, eps = 0.1, mod=c('d','d','d'), robust=0)
{

if((is.vector(mod)==FALSE) | (length(mod)!=3) )
{stop('The \'mod\' argument must be a vector of length 3.')}

if(!(mod[1] %in% c('d','c','l')))
{stop('The \'mod\' argument for the A component must be \'d\'(dynamic), \'c\'(constant) or \'l\'(linear).')}

if(!(mod[2] %in% c('d','c','l')))
{stop('The \'mod\' argument for the C component must be \'d\'(dynamic), \'c\'(constant) or \'l\'(linear).')}

if(!(mod[3] %in% c('d','c','l')))
{stop('The \'mod\' argument for the E component must be \'d\'(dynamic), \'c\'(constant) or \'l\'(linear).')}
  
if((knot_a<3)|(knot_c<3)|(knot_e<3)|(knot_a>40)|(knot_c>40)|(knot_e>40))
{stop('The number of knots must be an integer within [3, 40].')}


num_m <- nrow(data_m)*2
num_d <- nrow(data_d)*2
pheno_m <- matrix(NA, num_m, 1)
pheno_d <- matrix(NA, num_d, 1)
pheno_m[seq(from=1, to=num_m, by=2),1] <- data_m[,1]
pheno_m[seq(from=2, to=num_m, by=2),1] <- data_m[,2]
pheno_d[seq(from=1, to=num_d, by=2),1] <- data_d[,1]
pheno_d[seq(from=2, to=num_d, by=2),1] <- data_d[,2]
T_m <- rep(data_m[,3],each=2)
T_d <- rep(data_d[,3],each=2)

mag <- var(pheno_m)
init_max <- log(mag)
init_min <- log(mag) - abs(log(mag))*1.2
limit <- 12
limit_e <- 10
# low_var <- mag/100000
low_var <- 1e-06
upp_var <- 100
var_ran_up <- 3
var_ran_lo <- 1

eps <- eps*2

order <- 3

penal_a <- 2
penal_c <- 2
penal_e <- 2
if(mod[1]=='c')
{penal_a <- 1}
if(mod[2]=='c')
{penal_c <- 1}
if(mod[3]=='c')
{penal_e <- 1}

t_int <- max(c(T_m,T_d))-min(c(T_m,T_d))

l_m_1 <- (max(c(T_m,T_d))-T_m)/t_int
l_m_2 <- (T_m-min(c(T_m,T_d)))/t_int
l_d_1 <- (max(c(T_m,T_d))-T_d)/t_int
l_d_2 <- (T_d-min(c(T_m,T_d)))/t_int

if(mod[1]=='d')
{
delta_a <- matrix(0, knot_a+order-2-penal_a, knot_a+order-2)
for(i in 1:nrow(delta_a))
{
	if(penal_a==2)
	{delta_a[i, i:(i+2)] <- c(1,-2,1)}else{
		delta_a[i, i:(i+1)] <- c(1,-1)
	}
}
D_a <- t(delta_a)%*%delta_a
knots_a <- seq(from=min(T_m, T_d), to=max(T_m, T_d), length.out=knot_a)
interval_a <- knots_a[2] - knots_a[1]
knots_a <- c(c(min(T_m, T_d)-interval_a*2,min(T_m, T_d)-interval_a), knots_a)
knots_a <- c(knots_a, c(max(T_m, T_d)+interval_a,max(T_m, T_d)+interval_a*2))
B_des_a_m <- splineDesign(knots_a, x=T_m, ord=order)
B_des_a_d <- splineDesign(knots_a, x=T_d, ord=order)
ei_a <- eigen(D_a)
B_des_a_m <- B_des_a_m%*%ei_a$vectors
B_des_a_d <- B_des_a_d%*%ei_a$vectors
D_a <- diag(c(ei_a$values[1:(length(ei_a$values)-2)],0,0))
}else{
if(mod[1]=='l')
{
D_a <- matrix(0,2,2)
B_des_a_m <- matrix(NA, num_m, 2)
B_des_a_m[,1] <- l_m_1
B_des_a_m[,2] <- l_m_2
B_des_a_d <- matrix(NA, num_d, 2)
B_des_a_d[,1] <- l_d_1
B_des_a_d[,2] <- l_d_2
knots_a <- c(min(T_m, T_d),max(T_m, T_d))
}else{
D_a <- matrix(0,1,1)
B_des_a_m <- matrix(1, num_m, 1)
B_des_a_d <- matrix(1, num_d, 1)
knots_a <- c(min(T_m, T_d))
}
}

if(mod[2]=='d')
{
delta_c <- matrix(0, knot_c+order-2-penal_c, knot_c+order-2)
for(i in 1:nrow(delta_c))
{
	if(penal_c==2)
	{delta_c[i, i:(i+2)] <- c(1,-2,1)}else{
		delta_c[i, i:(i+1)] <- c(1,-1)
	}
}
D_c <- t(delta_c)%*%delta_c

knots_c <- seq(from=min(T_m, T_d), to=max(T_m, T_d), length.out=knot_c)
interval_c <- knots_c[2] - knots_c[1]
knots_c <- c(c(min(T_m, T_d)-interval_c*2,min(T_m, T_d)-interval_c), knots_c)
knots_c <- c(knots_c, c(max(T_m, T_d)+interval_c,max(T_m, T_d)+interval_c*2))
B_des_c_m <- splineDesign(knots_c, x=T_m, ord=order)
B_des_c_d <- splineDesign(knots_c, x=T_d, ord=order)
ei_c <- eigen(D_c)
B_des_c_m <- B_des_c_m%*%ei_c$vectors
B_des_c_d <- B_des_c_d%*%ei_c$vectors
D_c <- diag(c(ei_c$values[1:(length(ei_c$values)-2)],0,0))
}else{
if(mod[2]=='l')
{
	B_des_c_m <- matrix(NA, num_m, 2)
	B_des_c_m[,1] <- l_m_1
	B_des_c_m[,2] <- l_m_2
	B_des_c_d <- matrix(NA, num_d, 2)
	B_des_c_d[,1] <- l_d_1
	B_des_c_d[,2] <- l_d_2
	D_c <- matrix(0,2,2)
	knots_c <- c(min(T_m, T_d),max(T_m, T_d))
}else{
	B_des_c_m <- matrix(1, num_m, 1)
	B_des_c_d <- matrix(1, num_d, 1)
	D_c <- matrix(0,1,1)
	knots_c <- c(min(T_m, T_d))
}
}

if(mod[3]=='d')
{
delta_e <- matrix(0, knot_e+order-2-penal_e, knot_e+order-2)
for(i in 1:nrow(delta_e))
{
	if(penal_e==2)
	{delta_e[i, i:(i+2)] <- c(1,-2,1)}else{
		delta_e[i, i:(i+1)] <- c(1,-1)
	}
}
D_e <- t(delta_e)%*%delta_e

knots_e <- seq(from=min(T_m, T_d), to=max(T_m, T_d), length.out=knot_e)
interval_e <- knots_e[2] - knots_e[1]
knots_e <- c(c(min(T_m, T_d)-interval_e*2,min(T_m, T_d)-interval_e), knots_e)
knots_e <- c(knots_e, c(max(T_m, T_d)+interval_e,max(T_m, T_d)+interval_e*2))
B_des_e_m <- splineDesign(knots_e, x=T_m, ord=order)
B_des_e_d <- splineDesign(knots_e, x=T_d, ord=order)
ei_e <- eigen(D_e)
B_des_e_m <- B_des_e_m%*%ei_e$vectors
B_des_e_d <- B_des_e_d%*%ei_e$vectors
D_e <- diag(c(ei_e$values[1:(length(ei_e$values)-2)],0,0))
}else{
if(mod[3]=='l')
{
	B_des_e_m <- matrix(NA, num_m, 2)
	B_des_e_m[,1] <- l_m_1
	B_des_e_m[,2] <- l_m_2
	B_des_e_d <- matrix(NA, num_d, 2)
	B_des_e_d[,1] <- l_d_1
	B_des_e_d[,2] <- l_d_2
	D_e <- matrix(0,2,2)
	knots_e <- c(min(T_m, T_d),max(T_m, T_d))
}else{
	B_des_e_m <- matrix(1, num_m, 1)
	B_des_e_d <- matrix(1, num_d, 1)
	D_e <- matrix(0,1,1)
	knots_e <- c(min(T_m, T_d))
}

}

n_a <- ncol(B_des_a_m)
n_c <- ncol(B_des_c_m)
n_e <- ncol(B_des_e_m)

lower <- 0

#if(is.na(init_var[1]))
#{
	var_b_a <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))
#}else{
#	if((init_var[1]<20)&(init_var[1]>low_var))
#	{var_b_a <- init_var[1]}else{var_b_a <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))}
#}

if(mod[1] %in% c('l','c'))
{
	var_b_a <- lower
	n_a <- ifelse(mod[1]=='l',2,1) 
}

#if(is.na(init_var[2]))
#{
	var_b_c <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))
#}else{
#	if((init_var[2]<20)&(init_var[2]>low_var))
#	{var_b_c <- init_var[2]}else{var_b_c <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))}
#}

if(mod[2] %in% c('l','c'))
{
	var_b_c <- lower
	n_c <- ifelse(mod[2]=='l',2,1) 
}

#if(is.na(init_var[3]))
#{
	var_b_e <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))
#}else{
#	if((init_var[3]<20)&(init_var[3]>low_var))
#	{var_b_e <- init_var[3]}else{var_b_e <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))}
#}

if(mod[3] %in% c('l','c'))
{
	var_b_e <- lower
	n_e <- ifelse(mod[3]=='l',2,1) 
}

beta_a <- runif(n_a,min=-0.2,max=0.2)
beta_c <- runif(n_c,min=-0.2,max=0.2)
beta_e <- runif(n_e,min=-0.2,max=0.2)

lik <- 100000
lik_pre <- 200000

liks <- c()
betas <- matrix(0,0,n_a+n_c+n_e)
vars <- matrix(0,0,3)
hessians <- matrix(0,0,9)

if((mod[1]!='d')&(mod[2]!='d')&(mod[3]!='d'))
{
	low_a <- -15
	upp_a <- 15
	low_c <- -15
	upp_c <- 15
	low_e <- -8
	upp_e <- 10
	result <- optim(c(beta_a,beta_c,beta_e), loglik_AtCtEt_epsp_g, gr_AtCtEt_epsp_g, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, var_b_a=0, var_b_c=0, var_b_e=0, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(rep(low_a,n_a),rep(low_c,n_c),rep(low_e,n_e)), upper = c(rep(upp_a,n_a),rep(upp_c,n_c),rep(upp_e,n_e)), method = "L-BFGS-B", control=list(maxit = 3000))
	beta_a <- result$par[1:n_a]
	beta_c <- result$par[(1+n_a):(n_a+n_c)]
	beta_e <- result$par[(1+n_a+n_c):(n_a+n_c+n_e)]
	lik <- loglik_AtCtEt_epsp(c(0,0,0), pheno_m=pheno_m, pheno_d=pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, beta_a=beta_a, D_a=D_a, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, beta_c=beta_c, D_c=D_c, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, beta_e=beta_e, D_e=D_e)		
	AtCtEtp_model <- list(D_a = D_a, D_c = D_c, D_e = D_e, pheno_m = pheno_m, pheno_d = pheno_d, T_m = T_m, T_d = T_d, knot_a=knots_a, knot_c=knots_c, knot_e=knots_e, beta_a=beta_a, beta_c=beta_c, beta_e=beta_e, con=result$convergence, lik=lik/2, iter=lik/2, var_b_a=lower, var_b_c=lower, var_b_e=lower, mod=mod)

}else{
	while(abs(lik-lik_pre)>eps)
	{
		lik_pre <- lik
		if((mod[1]=='d')&(mod[2]=='d')&(mod[3]=='d'))
		{
			result <- optim(c(beta_a,beta_c,beta_e), loglik_AtCtEt_epsp_g, gr_AtCtEt_epsp_g, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, var_b_a=var_b_a, var_b_c=var_b_c, var_b_e=var_b_e, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(rep((-1)*limit,n_a+n_c),rep((-1)*limit_e,n_e)), upper = rep(limit,n_a+n_c+n_e), method = "L-BFGS-B", control=list(maxit = 3000))
			betas <- rbind(betas, result$par)
			beta_a <- result$par[1:n_a]
			beta_c <- result$par[(1+n_a):(n_a+n_c)]
			beta_e <- result$par[(1+n_a+n_c):(n_a+n_c+n_e)]
			result <- optim(c(var_b_a,var_b_c,var_b_e), loglik_AtCtEt_epsp, gr_AtCtEt_epsp, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, beta_a=beta_a, beta_c=beta_c, beta_e = beta_e, D_a=D_a, D_c=D_c, D_e=D_e, lower = rep(low_var,3), upper = rep(upp_var,3), method = "L-BFGS-B", control=list(maxit = 3000), hessian=TRUE)
      
		}else
		{	
			v_a_t <- var_b_a
			v_c_t <- var_b_c
			v_e_t <- var_b_e
			
			if(mod[1]!='d')
			{v_a_t <- 0}
			if(mod[2]!='d')
			{v_c_t <- 0}
			if(mod[3]!='d')
			{v_e_t <- 0}
			
			low_a <- (-1)*limit
			upp_a <- limit
			low_c <- (-1)*limit
			upp_c <- limit
			low_e <- (-1)*limit_e
			upp_e <- limit_e
			
			result <- optim(c(beta_a,beta_c,beta_e), loglik_AtCtEt_epsp_g, gr_AtCtEt_epsp_g, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, var_b_a=v_a_t, var_b_c=v_c_t, var_b_e=v_e_t, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(rep(low_a,n_a),rep(low_c,n_c),rep(low_e,n_e)), upper = c(rep(upp_a,n_a),rep(upp_c,n_c),rep(upp_e,n_e)), method = "L-BFGS-B", control=list(maxit = 3000))
			beta_a <- result$par[1:n_a]
			beta_c <- result$par[(1+n_a):(n_a+n_c)]
			beta_e <- result$par[(1+n_a+n_c):(n_a+n_c+n_e)]
			
			betas <- rbind(betas,c(beta_a,beta_c,beta_e))

			low_a <- low_c <- low_e <- low_var
			upp_a <- upp_c <- upp_e <- upp_var
			if((mod[1]=='l')|(mod[1]=='c'))
			{low_a <- upp_a <- lower}
			if((mod[2]=='l')|(mod[2]=='c'))
			{low_c <- upp_c <- lower}
			if((mod[3]=='l')|(mod[3]=='c'))
			{low_e <- upp_e <- lower}
		
			result <- optim(c(var_b_a,var_b_c,var_b_e), loglik_AtCtEt_epsp, gr_AtCtEt_epsp, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, beta_a=beta_a, beta_c=beta_c, beta_e=beta_e, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(low_a, low_c, low_e), upper = c(upp_a, upp_c, upp_e), method = "L-BFGS-B", control=list(maxit = 3000), hessian=TRUE)	
		}
		vars <- rbind(vars, result$par)
    hessians <- rbind(hessians, c(result$hessian))

		var_b_a <- result$par[1]
		var_b_c <- result$par[2]
		var_b_e <- result$par[3]
		lik <- result$value
		liks <- c(liks, result$value)
	}

	min_i <- match(min(liks), liks)
	
	if(robust>0)
	{
    for(rob in 1:ceiling(robust))
	  {
      lik <- 100000
      lik_pre <- 200000
      
      liks_r <- c()
      betas_r <- matrix(0,0,n_a+n_c+n_e)
      vars_r <- matrix(0,0,3)
      hessians_r <- matrix(0,0,9)
      
      beta_a <- runif(n_a,min=init_min,max=init_max)
      beta_c <- runif(n_c,min=init_min,max=init_max)
      beta_e <- runif(n_e,min=init_min,max=init_max)
      if(var_b_a!=0)
      {var_b_a <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))}
      if(var_b_c!=0)
      {var_b_c <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))}
      if(var_b_e!=0)
      {var_b_e <- runif(1,min=var_ran_lo*abs(log(mag)),max=var_ran_up*abs(log(mag)))}
      
      while(abs(lik-lik_pre)>eps)
	    {
	    	lik_pre <- lik
	    	if((mod[1]=='d')&(mod[2]=='d')&(mod[3]=='d'))
	    	{
	    		result <- optim(c(beta_a,beta_c,beta_e), loglik_AtCtEt_epsp_g, gr_AtCtEt_epsp_g, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, var_b_a=var_b_a, var_b_c=var_b_c, var_b_e=var_b_e, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(rep((-1)*limit,n_a+n_c),rep((-1)*limit_e,n_e)), upper = rep(limit,n_a+n_c+n_e), method = "L-BFGS-B", control=list(maxit = 3000))
		    	betas_r <- rbind(betas_r, result$par)
	    		beta_a <- result$par[1:n_a]
	    		beta_c <- result$par[(1+n_a):(n_a+n_c)]
	    		beta_e <- result$par[(1+n_a+n_c):(n_a+n_c+n_e)]
	    		result <- optim(c(var_b_a,var_b_c,var_b_e), loglik_AtCtEt_epsp, gr_AtCtEt_epsp, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, beta_a=beta_a, beta_c=beta_c, beta_e = beta_e, D_a=D_a, D_c=D_c, D_e=D_e, lower = rep(low_var,3), upper = rep(upp_var,3), method = "L-BFGS-B", control=list(maxit = 3000), hessian=TRUE)
      
	    	}else
	    	{	
		    	v_a_t <- var_b_a
		    	v_c_t <- var_b_c
		    	v_e_t <- var_b_e
			
		    	if(mod[1]!='d')
		    	{v_a_t <- 0}
		    	if(mod[2]!='d')
		    	{v_c_t <- 0}
		    	if(mod[3]!='d')
		    	{v_e_t <- 0}
			
		    	low_a <- (-1)*limit
			    upp_a <- limit
			    low_c <- (-1)*limit
			    upp_c <- limit
			    low_e <- (-1)*limit_e
			    upp_e <- limit_e
			
		    	result <- optim(c(beta_a,beta_c,beta_e), loglik_AtCtEt_epsp_g, gr_AtCtEt_epsp_g, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, var_b_a=v_a_t, var_b_c=v_c_t, var_b_e=v_e_t, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(rep(low_a,n_a),rep(low_c,n_c),rep(low_e,n_e)), upper = c(rep(upp_a,n_a),rep(upp_c,n_c),rep(upp_e,n_e)), method = "L-BFGS-B", control=list(maxit = 3000))
	    		beta_a <- result$par[1:n_a]
	    		beta_c <- result$par[(1+n_a):(n_a+n_c)]
	    		beta_e <- result$par[(1+n_a+n_c):(n_a+n_c+n_e)]
			
	    		betas_r <- rbind(betas_r,c(beta_a,beta_c,beta_e))

	    		low_a <- low_c <- low_e <- low_var
	    		upp_a <- upp_c <- upp_e <- upp_var
	    		if((mod[1]=='l')|(mod[1]=='c'))
	    		{low_a <- upp_a <- lower}
	    		if((mod[2]=='l')|(mod[2]=='c'))
	    		{low_c <- upp_c <- lower}
	    		if((mod[3]=='l')|(mod[3]=='c'))
	    		{low_e <- upp_e <- lower}
		
		    	result <- optim(c(var_b_a,var_b_c,var_b_e), loglik_AtCtEt_epsp, gr_AtCtEt_epsp, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, beta_a=beta_a, beta_c=beta_c, beta_e=beta_e, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(low_a, low_c, low_e), upper = c(upp_a, upp_c, upp_e), method = "L-BFGS-B", control=list(maxit = 3000), hessian=TRUE)	
	
	    	}
    		vars_r <- rbind(vars_r, result$par)
        hessians_r <- rbind(hessians_r, c(result$hessian))

    		var_b_a <- result$par[1]
    		var_b_c <- result$par[2]
    		var_b_e <- result$par[3]
	    	lik <- result$value
    		liks_r <- c(liks_r, result$value)
	    }

      if(min(liks_r)<min(liks))
      {
        liks <- liks_r
        vars <- vars_r
        betas <- betas_r
        hessians <- hessians_r
      }
    }
	  
	  min_i <- match(min(liks), liks)
	}
	
	if((mod[1]=='d')&(mod[2]=='d')&(mod[3]=='d'))
	{
		result <- optim(betas[min_i,], loglik_AtCtEt_epsp_g, gr_AtCtEt_epsp_g, pheno_m = matrix(pheno_m), pheno_d = matrix(pheno_d), B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, var_b_a=vars[min_i,1], var_b_c=vars[min_i,2], var_b_e=vars[min_i,3], D_a=D_a, D_c=D_c, D_e=D_e, lower = c(rep((-1)*limit,n_a+n_c),rep((-1)*limit_e,n_e)), upper = rep(limit,n_c+n_a+n_e), method = "L-BFGS-B", control=list(maxit = 3000))
	}else{
		v_a_t <- vars[min_i,1]
		v_c_t <- vars[min_i,2]
		v_e_t <- vars[min_i,3]
		
		if(mod[1]!='d')
		{v_a_t <- lower}
		if(mod[2]!='d')
		{v_c_t <- lower}
		if(mod[3]!='d')
		{v_e_t <- lower}
		
		low_a <- (-1)*limit
			upp_a <- limit
			low_c <- (-1)*limit
			upp_c <- limit
			low_e <- (-1)*limit_e
			upp_e <- limit_e
		
		result <- optim(betas[min_i,], loglik_AtCtEt_epsp_g, gr_AtCtEt_epsp_g, pheno_m = pheno_m, pheno_d = pheno_d, B_des_a_m=B_des_a_m, B_des_a_d=B_des_a_d, B_des_c_m=B_des_c_m, B_des_c_d=B_des_c_d, B_des_e_m=B_des_e_m, B_des_e_d=B_des_e_d, var_b_a=v_a_t, var_b_c=v_c_t, var_b_e=v_e_t, D_a=D_a, D_c=D_c, D_e=D_e, lower = c(rep(low_a,n_a),rep(low_c,n_c),rep(low_e,n_e)), upper = c(rep(upp_a,n_a),rep(upp_c,n_c),rep(upp_e,n_e)), method = "L-BFGS-B", control=list(maxit = 3000))
	}
	beta_a <- result$par[1:n_a]
		beta_c <- result$par[(1+n_a):(n_a+n_c)]
		beta_e <- result$par[(1+n_a+n_c):(n_a+n_c+n_e)]

		if(mod[1]=='d')
		{
			D_a <- t(delta_a)%*%delta_a
			beta_a <- ei_a$vectors%*%beta_a
		}
		if(mod[2]=='d')
		{
			D_c <- t(delta_c)%*%delta_c
			beta_c <- ei_c$vectors%*%beta_c
		}
		if(mod[3]=='d')
		{
			D_e <- t(delta_e)%*%delta_e
			beta_e <- ei_e$vectors%*%beta_e
		}
		
		AtCtEtp_model <- list(D_a = D_a, D_c = D_c, D_e=D_e, pheno_m = pheno_m, pheno_d = pheno_d, T_m = T_m, T_d = T_d, knot_a=knots_a, knot_c=knots_c, knot_e=knots_e, beta_a=beta_a, beta_c=beta_c, beta_e=beta_e, con=result$convergence, lik=min(liks)/2, iter=(liks)/2, var_b_a=vars[min_i,1], var_b_c=vars[min_i,2], var_b_e=vars[min_i,3], mod=mod, hessian = matrix(hessians[min_i,],3,3)/2)	
	
}


class(AtCtEtp_model) <- 'AtCtEtp_model'

#print('Estimates of beta_a:')
#print(beta_a)
#print('Estimates of beta_c:')
#print(beta_c)
#print('Estimates of beta_e:')
#print(beta_e)

return(invisible(AtCtEtp_model))

}

Try the ACEt package in your browser

Any scripts or data that you put into this service are public.

ACEt documentation built on May 31, 2023, 9:25 p.m.