## File Name: srm.R
## File Version: 0.625
srm <- function(model.syntax = NULL,
data = NULL, group.var = NULL,
rrgroup_name = NULL, person_names = c("Actor","Partner"),
fixed.groups=FALSE, var_positive=-1, optimizer="srm",
maxiter=300, conv_dev=1e-8, conv_par=1e-6,
do_line_search=TRUE, line_search_iter_max=6,
verbose=TRUE, use_rcpp=TRUE, shortcut=TRUE, use_woodbury=TRUE)
{
CALL <- match.call()
v1 <- s1 <- Sys.time()
method <- "ml"
#method <- "uls"
inv_type <- "pinv"
sd_noise <- 0
# use_woodbury <- FALSE
z0 <- Sys.time()
## Step 0: Determine the number of groups
if (is.null(group.var)) {
ngroups = 1L
} else {
ngroups = length(unique(data[,c(group.var)]))
}
## Step 1: use model syntax to generate a parameter table
parm.table <- SRM_PARTABLE_MAKE(model.syntax = model.syntax, ngroups = ngroups,
data_colnames = colnames(data), method=method )
rrvar_names <- attr(parm.table, "rrvar_names")
personcov_names <- attr(parm.table, "personcov_names")
dyadcov_names <- attr(parm.table, "dyadcov_names")
# cat(" *** make partable ") ; z1 <- Sys.time(); print(z1-z0) ; z0 <- z1
#-- add more information to parameter table
if (sd_noise>0){
parm.table <- SRM_PARTABLE_ADD_NOISE(parm.table=parm.table, sd=sd_noise)
}
res <- SRM_PARTABLE_EXTEND( parm.table=parm.table, var_positive=var_positive,
optimizer=optimizer, method=method )
parm.table <- res$parm.table
parm_table_free <- res$parm_table_free
lower <- res$lower
upper <- res$upper
NOP <- res$NOP
npar <- res$npar
symm_matrices <- res$symm_matrices
optimizer <- res$optimizer
# cat(" *** refine partable ") ; z1 <- Sys.time(); print(z1-z0) ; z0 <- z1
# create lists of model parameters
parm_list <- SRM_CREATE_PARAMETER_LISTS(parm.table = parm.table, ngroups = ngroups )
parm_list <- SRM_INCLUDE_PARAMETERS_PARM_LIST( parm.table=parm.table,
parm_list=parm_list, symm_matrices=symm_matrices, include_fixed=TRUE )
# cat(" *** create parameter lists") ; z1 <- Sys.time(); print(z1-z0) ; z0 <- z1
## Step 2: make list with data
data_list <- SRM_MAKE_DATA_LIST(srm_data = data, person_names = person_names ,
rrgroup_name = rrgroup_name, group.var = group.var,
fixed.groups = fixed.groups, rrvar_names=rrvar_names,
personcov_names=personcov_names, dyadcov_names=dyadcov_names,
use_rcpp=TRUE, do_checks=FALSE)
rrgroup_name <- attr(data_list, "rrgroup_name")
# cat(" *** make data list") ; z1 <- Sys.time(); print(z1-z0) ; z0 <- z1
# make design matrices
data_list <- SRM_CREATE_DESIGN_MATRICES( data_list = data_list,
ngroups = ngroups, rrgroup_name=rrgroup_name, use_rcpp=use_rcpp)
# cat(" *** create design matrices") ; z1 <- Sys.time(); print(z1-z0) ; z0 <- z1
## Step 3: estimate the model
x0 <- parm.table[ attr(parm.table, "parm_extract"), "est"]
parm_list <- SRM_INCLUDE_PARAMETERS_PARM_LIST( parm.table=parm.table,
parm_list=parm_list, symm_matrices=symm_matrices, include_fixed=TRUE )
#- define optimization function for stats::nlminb()
srm_optim_function <- function(x, parm.table, parm_list, symm_matrices,
npar, data_list, NOP, parm_table_free, only_value=FALSE, only_gradient=FALSE)
{
compute_gradient <- ! only_value
res <- SRM_OPTIMIZE(x=x, parm.table=parm.table, parm_list=parm_list,
symm_matrices=symm_matrices, npar=npar, data_list=data_list,
NOP=NOP, parm_table_free=parm_table_free, use_rcpp=use_rcpp,
shortcut=shortcut, compute_gradient=compute_gradient,
compute_hessian=FALSE, method=method, inv_type=inv_type,
use_woodbury=use_woodbury)
if (only_value){
val <- -res$ll_new
} else {
if (only_gradient){
val <- -res$grad
} else {
val <- list( objective=-res$ll_new, gradient=-res$grad)
}
}
return(val)
}
#- optimization function for external optimizers
srm_objective <- function(x, parm.table, parm_list, symm_matrices,
npar, data_list, NOP, parm_table_free, method)
{
res <- srm_optim_function( x=x, parm.table=parm.table, parm_list=parm_list,
symm_matrices=symm_matrices, npar=npar, data_list=data_list, NOP=NOP,
parm_table_free=parm_table_free, only_value=TRUE, only_gradient=FALSE)
return(res)
}
srm_gradient <- function(x, parm.table, parm_list, symm_matrices,
npar, data_list, NOP, parm_table_free, method)
{
res <- srm_optim_function( x=x, parm.table=parm.table, parm_list=parm_list,
symm_matrices=symm_matrices, npar=npar, data_list=data_list, NOP=NOP,
parm_table_free=parm_table_free, only_value=FALSE, only_gradient=TRUE)
return(res)
}
# collect optimization function arguments
args_opt0 <- list( parm.table=parm.table, parm_list=parm_list,
symm_matrices=symm_matrices, npar=npar, data_list=data_list,
NOP=NOP, parm_table_free=parm_table_free)
args_fs <- args_opt0
args_fs$use_rcpp <- use_rcpp
args_fs$shortcut <- shortcut
args_fs$compute_hessian <- TRUE
args_fs$use_woodbury <- use_woodbury
s2 <- Sys.time()
time_pre <- s2-s1
# args_fs$method <- method
# args_opt0$method <- method
#--- optimization
if (optimizer=="srm"){
#*** internal optimizer: optimization of the log-likelihood with Fisher scoring
res <- SRM_OPTIMIZER_INTERNAL(x0=x0, fn=SRM_OPTIMIZE, args_fs=args_fs,
conv_dev=conv_dev, conv_par=conv_par, maxiter=maxiter,
do_line_search=do_line_search, line_search_iter_max=line_search_iter_max,
verbose=verbose)
} else {
#*** external optimizer: stats::nlminb
res <- SRM_OPTIMIZER_EXTERNAL( optimizer=optimizer, x0=x0, fn=srm_objective,
gr=srm_gradient, args_opt0=args_opt0, maxiter=maxiter, lower=lower,
upper=upper, conv_dev=conv_dev, conv_par=conv_par, verbose=verbose )
}
res_opt <- res
time_opt <- res$time_diff
coef <- x0 <- res$par
names(coef) <- attr(parm.table, "par_names")
#- compute Hessian matrix
s1 <- Sys.time()
args_fs$x <- x0
res <- do.call(what=SRM_OPTIMIZE, args=args_fs)
parm.table <- as.data.frame(res$parm.table)
loglike <- res$ll_new
dev <- -2*loglike
grad <- res$grad
grad_maxabs <- max(abs(grad)) / abs(loglike)
infomat <- res$expected_infomat
vcov <- SRM_GINV(x=-infomat)
rownames(vcov) <- colnames(vcov) <- names(coef)
parm_list <- res$parm_list
# compute standard errors
se <- sqrt(diag(vcov))
parm.table$se <- se[ parm.table$index ]
#- informations about input data
res <- SRM_INFO_INPUT_DATA(data_list=data_list)
nrr <- res$nrr
npersons <- res$npersons
ndyads <- res$ndyads
# covariance matrices at levels
sigma <- list()
for (level in c("U","D")){
for (gg in 1:ngroups){
sigma[[level]][[gg]] <- SRM_COMPUTE_SEM_SIGMA( parm_list=parm_list,
level=level, group=gg )
}
}
s2 <- Sys.time()
time <- list(time_pre=time_pre, time_opt=time_opt, time_post=s2-s1)
#--- output
res <- list( coef=coef, vcov=vcov, grad=grad, se=se, loglike=loglike, dev=dev,
res_opt=res_opt, parm.table=parm.table, parm_list=parm_list,
data_list=data_list, ngroups=ngroups, nrr=nrr, npersons=npersons,
ndyads=ndyads, grad_maxabs=grad_maxabs, sigma=sigma,
CALL=CALL, time=time, time_start=v1 )
class(res) <- 'srm'
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.