Skip to content

Commit

Permalink
Merge pull request #27 from sfcheung/devel
Browse files Browse the repository at this point in the history 
Update to v0.0.0.4
  • Loading branch information
@sfcheung
sfcheung authored Nov 7, 2020
2 parents dbf080d + 4f7a413 commit 5218ab8
Show file tree
Hide file tree
Showing 19 changed files with 738 additions and 57 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: semlbci
Title: Likelihood Based Confidence Interval for SEM
Version: 0.0.0.3
Version: 0.0.0.4
Authors@R:
person(given = "Shu Fai",
family = "Cheung",
Expand Down
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,9 @@ export(ci_bound_i)
export(ci_i)
export(find_free)
export(find_variance_in_free)
export(get_lhs_op_rhs)
export(semlbci)
export(set_constraint)
export(set_start)
export(syntax_to_i)
export(update_model)
8 changes: 8 additions & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,11 @@
# semlbci 0.0.0.4

* (Work in progress .... )

* Can find the LBCI for the free and user-defined parameters in the standardized solution.

* Starting values based on Wald CIs are uesd also for standardized solution.

# semlbci 0.0.0.3

* Can find the LBCI for free parameters and user-defined parameters.
Expand Down
184 changes: 158 additions & 26 deletions R/ci_bound_i.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#' value but seems to be necessary when we need to find the LBCI for user-defined parameter.
#' @param wald_ci_start If TRUE and there are no equality constraints in the model, Wald confidence limit will
#' be used as the starting value for the target parameter if it is not a userd-defined paramter.
#' @param standardized If TRUE, the LBCI is for the standardized estimate.
#' @param opts Options to be passed to \code{nloptr}
#' @param ... Optional arguments. Not used.
#'
Expand All @@ -48,6 +49,7 @@ ci_bound_i <- function(i = NULL,
perturbation_factor = .9,
lb_var = -Inf,
wald_ci_start = TRUE,
standardized = FALSE,
opts = list(
"algorithm" = "NLOPT_LD_SLSQP",
"xtol_rel" = 1.0e-10,
Expand All @@ -60,47 +62,177 @@ ci_bound_i <- function(i = NULL,
# Check if the parameter is a user-defined parameter
p_table <- lavaan::parameterTable(sem_out)
i_op <- p_table[i, "op"]
if (i_op == ":=") {
# Get the name of the defined parameter
i_name <- p_table[i, "label"]
if (standardized && (i_op %in% c("=~", "~", "~~", ":="))) {
p_std <- lavaan::standardizedSolution(sem_out,
type = "std.all",
se = FALSE,
zstat = FALSE,
pvalue = FALSE,
ci = FALSE,
remove.eq = FALSE,
remove.ineq = FALSE,
remove.def = FALSE,
output = "data.frame")
p_std$id <- seq_len(nrow(p_std))
i_lor <- get_lhs_op_rhs(i, sem_out)
i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
start0 <- lavaan::parameterTable(sem_out)
# The function to be minimized.
lbci_b_f <- function(param, sem_out, debug, lav_warn) {
# The following block, from set_constraint, is not necessary.
# def.function is defined by the model and does not change with
# the parameter values
# start0 <- lavaan::parameterTable(sem_out)
# start0[find_free(sem_out), "est"] <- param
# if (lav_warn) {
# fit2 <- lavaan::update(sem_out, start = start0, do.fit = FALSE)
# } else {
# suppressWarnings(fit2 <- lavaan::update(sem_out, start = start0, do.fit = FALSE))
# }
k * sem_out@Model@def.function(param)[i_name]
start1 <- start0
start1[start1$free > 0, "est"] <- param
sem_out2 <- sem_out
sem_out2@ParTable <- as.list(start1)
sem_model <- sem_out2@Model
sem_model <- update_model(sem_model,
start1[start1$free > 0, "est"] )
sem_out2@Model <- sem_model
std0 <- lavaan::standardizedSolution(sem_out2,
type = "std.all",
se = FALSE,
zstat = FALSE,
pvalue = FALSE,
ci = FALSE,
remove.eq = FALSE,
remove.ineq = FALSE,
remove.def = FALSE,
output = "data.frame")
k * std0[i_std, "est.std"]
}
# The gradient of the function to be minimized
lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
numDeriv::grad(lbci_b_f, param, sem_out = sem_out,
debug = debug, lav_warn = lav_warn)
debug = debug, lav_warn = lav_warn)
}
} else {
# The function to be minimized.
lbci_b_f <- function(param, sem_out, debug, lav_warn) {
k * param[i]
}
if (i_op == ":=") {
if (standardized) {
# # To be deleted
# # Copied from the block for standardized free parameters.
# # If works, should merge the two blocks together
# p_std <- lavaan::standardizedSolution(sem_out,
# type = "std.all",
# se = FALSE,
# zstat = FALSE,
# pvalue = FALSE,
# ci = FALSE,
# remove.eq = FALSE,
# remove.ineq = FALSE,
# remove.def = FALSE,
# output = "data.frame")
# p_std$id <- seq_len(nrow(p_std))
# i_lor <- get_lhs_op_rhs(i, sem_out)
# i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
# start0 <- lavaan::parameterTable(sem_out)
# # The function to be minimized.
# lbci_b_f <- function(param, sem_out, debug, lav_warn) {
# start1 <- start0
# start1[start1$free > 0, "est"] <- param
# sem_out2 <- sem_out
# sem_out2@ParTable <- as.list(start1)
# sem_model <- sem_out2@Model
# sem_model <- update_model(sem_model,
# start1[start1$free > 0, "est"] )
# sem_out2@Model <- sem_model
# std0 <- lavaan::standardizedSolution(sem_out2,
# type = "std.all",
# se = FALSE,
# zstat = FALSE,
# pvalue = FALSE,
# ci = FALSE,
# remove.eq = FALSE,
# remove.ineq = FALSE,
# remove.def = FALSE,
# output = "data.frame")
# k * std0[i_std, "est.std"]
# }
# # The gradient of the function to be minimized
# lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
# numDeriv::grad(lbci_b_f, param, sem_out = sem_out,
# debug = debug, lav_warn = lav_warn)
# }
} else {
# Get the name of the defined parameter
i_name <- p_table[i, "label"]
# The function to be minimized.
lbci_b_f <- function(param, sem_out, debug, lav_warn) {
k * sem_out@Model@def.function(param)[i_name]
}
# The gradient of the function to be minimized
lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
numDeriv::grad(lbci_b_f, param, sem_out = sem_out,
debug = debug, lav_warn = lav_warn)
}
}
# The gradient of the function to be minimized
grad_c <- rep(0, npar)
grad_c[i] <- k
lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
grad_c
} else {
if (standardized) {
# # To be deleted:
# p_std <- lavaan::standardizedSolution(sem_out,
# type = "std.all",
# se = FALSE,
# zstat = FALSE,
# pvalue = FALSE,
# ci = FALSE,
# remove.eq = FALSE,
# remove.ineq = FALSE,
# remove.def = FALSE,
# output = "data.frame")
# p_std$id <- seq_len(nrow(p_std))
# i_lor <- get_lhs_op_rhs(i, sem_out)
# i_std <- merge(p_std, i_lor, by = c("lhs", "op", "rhs"))$id
# start0 <- lavaan::parameterTable(sem_out)
# # The function to be minimized.
# lbci_b_f <- function(param, sem_out, debug, lav_warn) {
# start1 <- start0
# start1[start1$free > 0, "est"] <- param
# sem_out2 <- sem_out
# sem_out2@ParTable <- as.list(start1)
# sem_model <- sem_out2@Model
# sem_model <- update_model(sem_model,
# start1[start1$free > 0, "est"] )
# sem_out2@Model <- sem_model
# std0 <- lavaan::standardizedSolution(sem_out2,
# type = "std.all",
# se = FALSE,
# zstat = FALSE,
# pvalue = FALSE,
# ci = FALSE,
# remove.eq = FALSE,
# remove.ineq = FALSE,
# remove.def = FALSE,
# output = "data.frame")
# k * std0[i_std, "est.std"]
# }
# # The gradient of the function to be minimized
# lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
# numDeriv::grad(lbci_b_f, param, sem_out = sem_out,
# debug = debug, lav_warn = lav_warn)
# }
} else {
# The function to be minimized.
lbci_b_f <- function(param, sem_out, debug, lav_warn) {
k * param[i]
}
# The gradient of the function to be minimized
grad_c <- rep(0, npar)
grad_c[i] <- k
lbci_b_grad <- function(param, sem_out, debug, lav_warn) {
grad_c
}
}
}
if (wald_ci_start & !sem_out@Model@eq.constraints) {
if (i_op == ":=") {
xstart <- set_start(i, sem_out, which)
xstart <- xstart[xstart$free > 0, "est"]
} else {
xstart <- set_start(i, sem_out, which)
xstart <- xstart[xstart$free > 0, "est"]
if (standardized) {
xstart <- set_start(i, sem_out, which, standardized)
xstart <- xstart[xstart$free > 0, "est"]
} else {
xstart <- set_start(i, sem_out, which)
xstart <- xstart[xstart$free > 0, "est"]
}
}
} else {
xstart <- perturbation_factor * lavaan::coef(sem_out)
Expand Down
35 changes: 35 additions & 0 deletions R/get_lhs_op_rhs.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
#'@title Get lhs, op, and rhs of a parameter from an SEM output
#'
#'@description Get lhs, op, and rhs of a parameter from an SEM output
#'
#'@details
#'
#' Currently supports \code{lavaan} output only.
#'
#'@return
#' A one row data frame of lhs, op, and rhs of the parameter in the parameter table.
#'
#' @param i The position of the target parameters as in the parameter table of lavaan.
#' @param sem_out The SEM output. Currently \code{lavaan} output only.
#'
#'@examples
#' library(lavaan)
#' data(cfa_two_factors)
#' mod <-
#' "
#' f1 =~ x1 + x2 + a*x3
#' f2 =~ x4 + a*x5 + equal('f1=~x2')*x6
#' f1 ~~ 0*f2
#' asq := a^2
#' "
#' fit <- sem(mod, cfa_two_factors)
#'@export

get_lhs_op_rhs <- function(i, sem_out) {
if (!inherits(sem_out, "lavaan")) {
stop("sem_out is not a supported object.")
}
ptable <- lavaan::parameterTable(sem_out)
# Do not check for
as.data.frame(as.list(ptable[i, c("lhs", "op", "rhs")]))
}
Loading

0 comments on commit 5218ab8

Please sign in to comment.