Nothing
#' Fitting a 1-D SPM model with constant parameters
#'
#' @description This function implements a analytical solution to estimate the parameters in the continuous SPM model by assuming all the parameters are constants.
#'
#' @param spm_data A dataset for the SPM model. See the STPM package for more details about the format.
#' @param a The initial value for the paramter \eqn{a}. The initial value will be predicted if not specified.
#' @param b The initial value for the paramter \eqn{b}. The initial value will be predicted if not specified.
#' @param q The initial value for the paramter \eqn{q}. The initial value will be predicted if not specified.
#' @param f The initial value for the paramter \eqn{f}. The initial value will be predicted if not specified.
#' @param f1 The initial value for the paramter \eqn{f_1}. The initial value will be predicted if not specified.
#' @param mu0 The initial value for the paramter \eqn{\mu_0} in the baseline hazard. The initial value will be predicted if not specified.
#' @param theta The initial value for the paramter \eqn{\theta} in the baseline hazard. The initial value will be predicted if not specified.
#' @param lower A vector of the lower bound of the parameters.
#' @param upper A vector of the upper bound of the parameters.
#' @param control A list of the control parameters for the optimization paramters.
#' @param global A logical variable indicating whether the MLSL (TRUE) or the L-BFGS (FALSE) algorithm is used for the optimization.
#' @param verbose A logical variable indicating whether initial information is printed.
#' @param ahessian A logical variable indicating whether the approximate (FALSE) or analytical (TRUE) Hessian is returned.
#' @return est The estimates of the parameters.
#' @return hessian The Hessian matrix of the estimates.
#' @return lik The minus log-likelihood.
#' @return con A number indicating the convergence. See the 'nloptr' package for more details.
#' @return message Extra message about the convergence. See the 'nloptr' package for more details.
#' @references He, L., Zhbannikov, I., Arbeev, K. G., Yashin, A. I., and Kulminski, A.M., 2017. Genetic stochastic process model for detecting pleiotropic and interaction effects with longitudinal data.
#' @export
#' @examples {
#' library(stpm)
#' dat <- simdata_cont(N=500)
#' colnames(dat) <- c("id", "xi", "t1", "t2", "y", "y.next")
#' res <- spm_con_1d(as.data.frame(dat), a=-0.05, b=2, q=1e-8, f=80, f1=90, mu0=1e-3, theta=0.08)
#'}
spm_con_1d <- function(spm_data, a = NA, b = NA, q = NA, f = NA, f1 = NA, mu0 = NA, theta = NA, lower = c(), upper = c(), control = list(xtol_rel = 1e-6), global = FALSE, verbose = TRUE, ahessian = FALSE)
{
#a <- -0.1
#b <- 1.2
#q <- 0.001
#f <- 25
#f1 <- 26
#mu0 <- 0.01
#theta <- 0
data <- spm_data
data[which(is.na(data$y.next)),'y.next'] <- data[which(is.na(data$y.next)),'y']
data$r0 <- 0
aggdata <-aggregate(data[,c('id','t1','t2','xi')], by=list(data$id), FUN='max')
tau <- aggdata$t2
n_j <- as.data.frame(table(data$id))[,2]
m0 <- data$y
r0 <- data$r0
yij <- data$y.next
tij <- data$t2
t0 <- data$t1
delta <- aggdata$xi
name_par <- c('a','b','q','f','f1','mu0','theta')
if(is.na(f))
{f <- mean(m0, na.rm=TRUE)}
if(is.na(f1))
{f1 <- mean(m0, na.rm=TRUE)}
if(is.na(a))
{a <- (-0.1)*abs(max(m0, na.rm=TRUE))/max(t0,na.rm=TRUE)}
if(is.na(mu0))
{mu0 <- sum(delta, na.rm=TRUE)/sum(tij-t0, na.rm=TRUE)}
if(is.na(theta))
{theta <- 0.00001}
if(is.na(b))
{b <- sqrt(abs(mean(yij-m0, na.rm=TRUE)))/abs(mean(tij-t0,na.rm=TRUE))}
if(is.na(q))
{q <- mu0/(0.25*f^2)}
if(a>=0)
{stop('The value for the argument \'a\' must be negative.')}
if(b<0)
{stop('The value for the argument \'b\' must be positive.')}
if(q<0)
{stop('The value for the argument \'q\' must be positive.')}
#if(mu0<0)
#{stop('The value for the argument \'mu0\' must be positive.')}
param <- c(a,b,q,f,f1,mu0,theta)
if(length(lower)==0)
{lower <- c(a*5, b*0.1, q*0.01, f*0.5, f1*0.5, mu0*0.1, theta)}
if(length(upper)==0)
{upper <- c(a*0.05, b*5, q*10, f*2, f1*2, mu0*10, theta)}
if(verbose == TRUE)
{
print('Initial values:')
print(param)
print('Lower bounds:')
print(lower)
print('Upper bounds:')
print(upper)
}
# re <- lik_con_1d(param, m0,r0,tau,yij,delta,tij, n_j, t0)
# re <- c(re, gr_con_1d(param, m0,r0,tau,yij,delta,tij, n_j, t0))
if(global == FALSE)
{
re <- lbfgs( x0=param,fn=lik_con_1d,gr=gr_con_1d,m0=m0,r0=r0,tau=tau,yij=yij,delta=delta,tij=tij, n_j=n_j, t0=t0,lower=lower,upper=upper, control = control)
}else{
re <- mlsl( x0=param,fn=lik_con_1d,gr=gr_con_1d,m0=m0,r0=r0,tau=tau,yij=yij,delta=delta,tij=tij, n_j=n_j, t0=t0,lower=lower,upper=upper, control = control)
}
hessian <- optim(re$par, lik_con_1d, gr_con_1d, m0,r0,tau,yij,delta,tij, n_j, t0, method = "L-BFGS-B", lower = re$par, upper = re$par, hessian=TRUE)
a_hes <- hessian$hessian
if(ahessian == TRUE)
{
a_hes <- hes_loglik(re$par,m0,r0,tau,yij,delta,tij, n_j,t0)
}
colnames(a_hes) <- name_par
rownames(a_hes) <- name_par
#par_re <- matrix(re$par, 1, 7)
#colnames(par_re) <- c('a','b','q','f','f1','mu0','theta')
##### Calculate p-values for estimates ####
coef <- re$par
names(coef) <- name_par
fi <- solve(a_hes) #### Fischer Information Matrix. Note we're not using the negative of the hessian here
stderr <- sqrt(diag(fi))
zscore <- coef/stderr
pvalue <- 2*(1 - pnorm(abs(zscore)))
results.par_re <- cbind(coef,stderr,zscore,pvalue)
colnames(results.par_re) <- c("Coeff.", "Std. Err.", "z", "p value")
list(est=results.par_re,hessian=a_hes, lik=re$value, con=re$convergence, message = re$message)
}
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.