diff --git a/.Rbuildignore b/.Rbuildignore index 0f76123..b3c2b29 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -4,3 +4,4 @@ ^docs$ ^pkgdown$ ^scratch$ +^codecov\.yml$ diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml new file mode 100644 index 0000000..8ee5250 --- /dev/null +++ b/.github/workflows/test-coverage.yaml @@ -0,0 +1,50 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master, test_codecov] + pull_request: + branches: [main, master, test_codecov] + +name: test-coverage + +jobs: + test-coverage: + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + + steps: + - uses: actions/checkout@v3 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::covr + needs: coverage + + - name: Test coverage + run: | + covr::codecov( + quiet = FALSE, + clean = FALSE, + install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") + ) + shell: Rscript {0} + + - name: Show testthat output + if: always() + run: | + ## -------------------------------------------------------------------- + find ${{ runner.temp }}/package -name 'testthat.Rout*' -exec cat '{}' \; || true + shell: bash + + - name: Upload test results + if: failure() + uses: actions/upload-artifact@v3 + with: + name: coverage-test-failures + path: ${{ runner.temp }}/package diff --git a/DESCRIPTION b/DESCRIPTION index b29603a..3932e74 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: dsem Type: Package Title: Fit Dynamic Structural Equation Models -Version: 1.0.1 +Version: 1.0.2 Date: 2024-01-17 Authors@R: c(person(given = "James", diff --git a/NEWS.md b/NEWS.md index 38844c4..0894130 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +# dsem 1.0.2 + +* Eliminate `eval` usage +* Add codecov Action and badge + # dsem 1.0.1 * Fix bug in `simulate.dsem` to keep up with changing interface in `dsem` diff --git a/R/dsem.R b/R/dsem.R index d8b3d83..4f6ed94 100644 --- a/R/dsem.R +++ b/R/dsem.R @@ -62,6 +62,8 @@ #' \item{tmb_inputs}{The list of inputs passed to \code{\link[TMB]{MakeADFun}}} #' \item{opt}{The output from \code{\link[stats]{nlminb}}} #' \item{sdrep}{The output from \code{\link[TMB]{sdreport}}} +#' \item{interal}{Objects useful for package function, i.e., all arguments +#' passed during the call} #' } #' #' @references @@ -134,7 +136,7 @@ function( sem, # Data = list( "RAM" = as.matrix(na.omit(ram[,1:4])), "RAMstart" = as.numeric(ram[,5]), - "familycode_j" = sapply(family, FUN=switch, "fixed"=0, "normal"=1 ), + "familycode_j" = sapply(family, FUN=switch, "fixed"=0, "normal"=1, "gamma"=4 ), "y_tj" = tsdata ) # Construct parameters @@ -183,11 +185,19 @@ function( sem, } obj = MakeADFun( data=Data, parameters=Params, random=Random, map=Map, DLL="dsem" ) if(control$quiet==FALSE) list_parameters(obj) + internal = list( + sem = sem, + tsdata = tsdata, + family = family, + estimate_delta0 = estimate_delta0, + control = control + ) out = list( "obj"=obj, "ram"=ram, "sem_full"=out$model, "tmb_inputs"=list("data"=Data, "parameters"=Params, "random"=Random, "map"=Map), - "call" = match.call() ) + #"call" = match.call(), + "internal" = internal ) # Export stuff if( control$run_model==FALSE ){ @@ -226,7 +236,9 @@ function( sem, if( isTRUE(control$getsd) ){ if( isTRUE(control$verbose) ) message("Running sdreport") Hess_fixed = optimHess( par=out$opt$par, fn=obj$fn, gr=obj$gr ) - out$sdrep = sdreport( obj, hessian.fixed=Hess_fixed ) + out$sdrep = sdreport( obj, + hessian.fixed = Hess_fixed, + getJointPrecision = control$getJointPrecision ) }else{ out$sdrep = NULL } @@ -263,6 +275,8 @@ function( sem, #' by hand (only helpful for advanced users to change starting values or restart at intended values) #' @param map list of fixed and mirrored parameters, constructed by \code{dsem} by default but available #' to override this default and then pass to \code{\link[TMB]{MakeADFun}} +#' @param getJointPrecision whether to get the joint precision matrix. Passed +#' to \code{\link[TMB]{sdreport}}. #' #' @return #' An S3 object of class "dsem_control" that specifies detailed model settings, @@ -280,7 +294,8 @@ function( nlminb_loops = 1, run_model = TRUE, use_REML = TRUE, parameters = NULL, - map = NULL ){ + map = NULL, + getJointPrecision = FALSE ){ # Return structure( list( @@ -294,7 +309,8 @@ function( nlminb_loops = 1, run_model = run_model, use_REML = use_REML, parameters = parameters, - map = map + map = map, + getJointPrecision = getJointPrecision ), class = "dsem_control" ) } @@ -466,7 +482,7 @@ function( object, # pull out objects for easy use obj = object$obj parfull = obj$env$parList() - tsdata = eval(object$call$tsdata) + tsdata = object$internal$tsdata # Extract parameters, and add noise as desired par_zr = outer( obj$env$last.par.best, rep(1,nsim) ) @@ -492,17 +508,26 @@ function( object, Q_kk = newrep$Q_kk tmp = rmvnorm_prec( newrep$delta_k + as.vector(newrep$xhat_tj), Q_kk, nsim=1 ) # Modify call - newcall = object$call + #newcall = object$call # Get control - newcall$control = eval(newcall$control) - newcall$control$parameters = newparfull - newcall$control$parameters$x_tj[] = tmp + #newcall$control = eval(newcall$control) + #newcall$control$parameters = newparfull + #newcall$control$parameters$x_tj[] = tmp # Rebuild model with new GMRF values - newcall$control$run_model = FALSE - newfit = eval(newcall) + #newcall$control$run_model = FALSE + #newfit = eval(newcall) + control = object$internal$control + control$parameters = newparfull + control$parameters$x_tj[] = tmp + control$run_model = FALSE + newfit = dsem( sem = object$internal$sem, + tsdata = object$internal$tsdata, + family = object$internal$family, + estimate_delta0 = object$internal$estimate_delta0, + control = control ) out[[r]] = newfit$obj$simulate()$y_tj }else{ - out[[r]] = obj$simulate( par_zr[,r] ) + out[[r]] = obj$simulate( par_zr[,r] )$y_tj } colnames(out[[r]]) = colnames(tsdata) tsp(out[[r]]) = attr(tsdata,"tsp") @@ -579,43 +604,43 @@ function( object, # https://stats.stackexchange.com/questions/1432/what-do-the-residuals-in-a-logistic-regression-mean # Normal deviance residuals - if( FALSE ){ - x = rnorm(10) - y = x + rnorm(10) - Glm = glm( y ~ 1 + x, family="gaussian") - mu = predict(Glm,type="response") - r1 = y - mu - r1 - resid(Glm) - } - # Poisson deviance residuals - if( FALSE ){ - x = rnorm(10) - y = rpois(10, exp(x)) - Glm = glm( y ~ 1 + x, family="poisson") - mu = predict(Glm,type="response") - # https://stats.stackexchange.com/questions/398098/formula-for-deviance-residuals-for-poisson-model-with-identity-link-function - r1 = sign(y - mu) * sqrt(2*(y*log((y+1e-10)/mu) - (y-mu))) - r1 - resid(Glm) - } - # Binomial deviance residuals - if( FALSE ){ - p = 0.5 - y = rbinom(10, prob=p, size=1) - Glm = glm( y ~ 1, family="binomial") - mu = predict(Glm, type="response") - r1 = sign(y - mu) * sqrt(-2*(((1-y)*log(1-mu) + y*log(mu)))) - r1 - resid(Glm) - } - # Gamma deviance residuals - if( FALSE ){ - mu = 1 - cv = 0.8 - y = rgamma( n=10, shape=1/cv^2, scale=mu*cv^2 ) - Glm = glm( y ~ 1, family=Gamma(link='log')) - mu = predict(Glm, type="response") - r1 = sign(y - mu) * sqrt(2 * ( (y-mu)/mu - log(y/mu) )) - r1 - resid(Glm) - } + #if( FALSE ){ + # x = rnorm(10) + # y = x + rnorm(10) + # Glm = glm( y ~ 1 + x, family="gaussian") + # mu = predict(Glm,type="response") + # r1 = y - mu + # r1 - resid(Glm) + #} + ## Poisson deviance residuals + #if( FALSE ){ + # x = rnorm(10) + # y = rpois(10, exp(x)) + # Glm = glm( y ~ 1 + x, family="poisson") + # mu = predict(Glm,type="response") + # # https://stats.stackexchange.com/questions/398098/formula-for-deviance-residuals-for-poisson-model-with-identity-link-function + # r1 = sign(y - mu) * sqrt(2*(y*log((y+1e-10)/mu) - (y-mu))) + # r1 - resid(Glm) + #} + ## Binomial deviance residuals + #if( FALSE ){ + # p = 0.5 + # y = rbinom(10, prob=p, size=1) + # Glm = glm( y ~ 1, family="binomial") + # mu = predict(Glm, type="response") + # r1 = sign(y - mu) * sqrt(-2*(((1-y)*log(1-mu) + y*log(mu)))) + # r1 - resid(Glm) + #} + ## Gamma deviance residuals + #if( FALSE ){ + # mu = 1 + # cv = 0.8 + # y = rgamma( n=10, shape=1/cv^2, scale=mu*cv^2 ) + # Glm = glm( y ~ 1, family=Gamma(link='log')) + # mu = predict(Glm, type="response") + # r1 = sign(y - mu) * sqrt(2 * ( (y-mu)/mu - log(y/mu) )) + # r1 - resid(Glm) + #} # Poisson: sign(y - mu) * sqrt(2*(ifelse(y==0, 0, y*log(y/mu)) - (y-mu))) # Binomial: -2 * ((1-y)*log(1-mu) + y*log(mu)) @@ -701,11 +726,18 @@ function( object, if(type=="link") out = parfull$x_tj if(type=="response") out = report$mu_tj }else{ - newcall = object$call - newcall$tsdata = newdata + #newcall = object$call + #newcall$tsdata = newdata # Rebuild model with new data - newcall$run_model = FALSE - newfit = eval(newcall) + #newcall$run_model = FALSE + #newfit = eval(newcall) + control = object$internal$control + control$run_model = FALSE + newfit = dsem( sem = object$internal$sem, + tsdata = newdata, + family = object$internal$family, + estimate_delta0 = object$internal$estimate_delta0, + control = object$internal$control ) # Optimize random effects given original MLE and newdata newfit$obj$fn( object$opt$par ) # Return predictor @@ -803,7 +835,7 @@ function( object, model = model[model[,2]==0,c(1,3,4)] out = sem( model, S = Sprime, - N = nrow(eval(object$call$tsdata)) ) + N = nrow(object$internal$tsdata) ) # pass out return(out) diff --git a/R/make_dsem_ram.R b/R/make_dsem_ram.R index ee05363..6b3b44f 100644 --- a/R/make_dsem_ram.R +++ b/R/make_dsem_ram.R @@ -199,21 +199,21 @@ function( sem, # MATH CHECK IN ROXYGEN DOCS ABOVE if( FALSE ){ - rho = 0.8 - sigma = 0.5 - Rho = Gamma = matrix(0, nrow=4, ncol=4) - Rho[cbind(2:4,1:3)] = rho - Gamma = I = diag(4) - diag(Gamma)[] = sigma - # DSEM covariance - solve(I-Rho) %*% Gamma %*% t(Gamma) %*% t(solve(I-Rho)) - # Stated covariance - sigma^2 * rbind( - c(1, rho, rho^2, rho^3), - c(rho, 1+rho^2, rho*(1+rho^2), rho^2*(1+rho^2) ), - c(rho^2, rho*(1+rho^2), 1+rho^2+rho^4, rho*(1+rho^2+rho^4) ), - c(rho^3, rho^2*(1+rho^2), rho*(1+rho^2+rho^4), 1+rho^2+rho^4+rho^6 ) - ) + #rho = 0.8 + #sigma = 0.5 + #Rho = Gamma = matrix(0, nrow=4, ncol=4) + #Rho[cbind(2:4,1:3)] = rho + #Gamma = I = diag(4) + #diag(Gamma)[] = sigma + ## DSEM covariance + #solve(I-Rho) %*% Gamma %*% t(Gamma) %*% t(solve(I-Rho)) + ## Stated covariance + #sigma^2 * rbind( + # c(1, rho, rho^2, rho^3), + # c(rho, 1+rho^2, rho*(1+rho^2), rho^2*(1+rho^2) ), + # c(rho^2, rho*(1+rho^2), 1+rho^2+rho^4, rho*(1+rho^2+rho^4) ), + # c(rho^3, rho^2*(1+rho^2), rho*(1+rho^2+rho^4), 1+rho^2+rho^4+rho^6 ) + #) } ####### Error checks diff --git a/R/zzz.R b/R/zzz.R index 9ec976c..9b2a058 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -4,7 +4,7 @@ # See: https://mail.google.com/mail/u/0/#inbox/QgrcJHsNjpqZNGJhNMVHcfGFDLLMfrvqqHl .onLoad <- function(libname, pkgname) { #checkDepPackageVersion(dep_pkg="TMB", this_pkg="dsem") - #checkDepPackageVersion(dep_pkg="Matrix", this_pkg="TMB") + checkDepPackageVersion(dep_pkg="Matrix", this_pkg="TMB") } .onUnload <- function(libpath) { diff --git a/README.md b/README.md index 54dc8ce..7a65cdb 100644 --- a/README.md +++ b/README.md @@ -26,3 +26,5 @@ Please see package vignettes for more details regarding syntax and features. [![](https://cranlogs.r-pkg.org/badges/dsem)](https://cran.r-project.org/package=dsem) [![DOI](https://zenodo.org/badge/656795688.svg)](https://zenodo.org/doi/10.5281/zenodo.10304770) + +[![Codecov test coverage](https://codecov.io/gh/James-Thorson-NOAA/dsem/branch/test_codecov/graph/badge.svg)](https://app.codecov.io/gh/James-Thorson-NOAA/dsem?branch=test_codecov) diff --git a/codecov.yml b/codecov.yml new file mode 100644 index 0000000..04c5585 --- /dev/null +++ b/codecov.yml @@ -0,0 +1,14 @@ +comment: false + +coverage: + status: + project: + default: + target: auto + threshold: 1% + informational: true + patch: + default: + target: auto + threshold: 1% + informational: true diff --git a/man/dsem.Rd b/man/dsem.Rd index 8bcef87..e5e38a1 100644 --- a/man/dsem.Rd +++ b/man/dsem.Rd @@ -40,6 +40,8 @@ An object (list) of class `dsem`. Elements include: \item{tmb_inputs}{The list of inputs passed to \code{\link[TMB]{MakeADFun}}} \item{opt}{The output from \code{\link[stats]{nlminb}}} \item{sdrep}{The output from \code{\link[TMB]{sdreport}}} +\item{interal}{Objects useful for package function, i.e., all arguments + passed during the call} } } \description{ diff --git a/tests/testthat/test-platform.R b/tests/testthat/test-platform.R index 915cf34..827399f 100644 --- a/tests/testthat/test-platform.R +++ b/tests/testthat/test-platform.R @@ -1,8 +1,4 @@ - -context("Testing cross platform and R version compatibility") - -# Eastern Bering Sea pollcok test_that("dsem example is working ", { #skip_on_ci() sem = " @@ -47,11 +43,37 @@ test_that("dsem example is working ", { Z = ts( cbind(Data, "neg_Gov_wage"=-1*Data[,'Gov_wage']) ) # Fit model - fit = dsem( sem=sem, tsdata=Z ) + fit = dsem( sem=sem, + tsdata=Z, + control = dsem_control(getJointPrecision=TRUE) ) # Check objective function expect_equal( as.numeric(fit$opt$obj), 587.4755, tolerance=1e-2 ) # Convert and plot using phylopath as_fitted_DAG(fit) + + # Various other utilities + plot(fit) + vcov(fit, which="fixed") + vcov(fit, which="random") + vcov(fit, which="both") + print(fit) + logLik(fit) + as_sem(fit) + predict(fit, type="link") + predict(fit, type="response") + predict(fit, type="link", newdata=Z) + simulate(fit, variance = "none") + simulate(fit, variance = "random") + simulate(fit, variance = "both") + simulate(fit, resimulate_gmrf=TRUE) + + # Refit with measurement errors + fit1 = dsem( sem=sem, + tsdata=Z, + family = c("normal","gamma",rep("fixed",ncol(Z)-2)), + control = dsem_control(getsd=FALSE, newton_loops=0) ) + residuals(fit1, type="deviance") + residuals(fit1, type="response") })