RJ 2025-4 issue

Ref #115

Author: mitchelloharawild

_articles/RJ-2025-032/RJ-2025-032.R

@@ -0,0 +1,146 @@
+# Generated by `rjournal_pdf_article()` using `knitr::purl()`: do not edit by hand
+# Please edit RJ-2025-032.Rmd to modify this file
+
+## ----setup, include=FALSE-----------------------------------------------------
+knitr::opts_chunk$set(
+  echo = TRUE,       # show code by default
+  warning = FALSE,   # suppress warnings
+  message = FALSE,   # suppress messages
+  fig.align = "center",  # center figures
+  fig.width = 6,     # default figure width in inches
+  fig.height = 4,  # default figure height in inches
+  dpi = 300,         # high-res figures for PDF
+  out.width = "100%",
+  cache = TRUE
+)
+options(csurvey.multicore = FALSE)
+library(Matrix)
+library(data.table)
+library(coneproj)
+#library(foreign)
+library(tidyverse)
+library(csurvey)
+library(MASS)
+library(survey)
+
+
+## -----------------------------------------------------------------------------
+library(csurvey)
+data(nhdat2, package = "csurvey")
+dstrat <- svydesign(ids = ~id,  strata = ~str, data = nhdat2,  weight = ~wt)
+
+
+## -----------------------------------------------------------------------------
+ans <- csvy(chol ~ incr(age), design = dstrat, n.mix = 100)
+
+
+## -----------------------------------------------------------------------------
+cat("CIC (constrained):", ans$CIC, "\n")
+
+
+## -----------------------------------------------------------------------------
+cat("CIC (unconstrained):", ans$CIC.un, "\n")
+
+
+## -----------------------------------------------------------------------------
+cat(svycontrast(ans, list(avg = c(rep(-1, 13)/13, rep(1, 12)/12))), "\n")
+
+
+## ----nh1big, fig.cap="Estimates of average cholesterol level for 25 ages, with 95% confidence intervals, for a stratified sample in the R dataset `nhdat2`, $n=1933$.", fig.align='center', echo=FALSE----
+knitr::include_graphics("figures/nhanes1.png")
+
+
+## -----------------------------------------------------------------------------
+set.seed(1)
+ans <- csvy(chol ~ incr(age)*incr(wcat)*icat, design = dstrat)
+
+
+## -----------------------------------------------------------------------------
+domains <- data.frame(age = c(24, 35), wcat = c(2, 4), icat = c(2, 3))
+pans <- predict(ans, newdata = domains, se.fit = TRUE)
+cat("Predicted values, confidence intervals and standard errors for specified domains:\n")
+print (pans)
+
+
+## ----nh2, fig.cap="Constrained estimates of population domain means for 400 domains in a 25x4x4 grid. The increasing population domain estimates for the 25 ages are shown within the waist size and income categories. The blue bands indicate 95% confidence intervals for the population domain means, with two specific domains, namely, (age, waist, income) = (24, 2, 2) and (35, 4, 3) marked in red. Empty domains are marked with a red 'x' sign.", fig.align='center', echo=FALSE----
+knitr::include_graphics("figures/nhanes_grid3.png")
+
+
+## ----nh2un, fig.cap="Unconstrained estimates of population domain means for 400 domains in a 25x4x4 grid. The population domain estimates for the 25 ages are shown within the waist size and income categories. The green bands indicate 95% confidence intervals for the population domain means. Empty domains are marked with a red 'x' sign.", fig.align='center', echo=FALSE----
+knitr::include_graphics("figures/nhanes_grid3_un.png")
+
+
+## ----surface9, fig.cap="Estimates of average log(salary) by field of study and year of degree, for observations where highest degree is a Bachelor's, for each of the nine regions.", fig.align='center', echo=FALSE, out.width="60%"----
+knitr::include_graphics("figures/new_surfaces9.png")
+
+
+## ----NEreg, fig.cap="Estimates of average log(salary) for the 75 domains in each of three regions. The blue dots represent the constrained domain mean estimates, while the grey dots represent the unconstrained domain mean estimates. The blue band is the 95% confidence interval for the domains, using the constraints; the grey band is the 95% unconstrained domain mean confidence interval.", fig.align='center', echo=FALSE, out.width="100%"----
+knitr::include_graphics("figures/newplot4.png")
+
+
+## ----test, fig.cap="Estimates of average log(salary) by father's education level, for each of five regions and four fields, for subjects whose degree was attained in 2016-2017. The solid blue lines connect the estimates where the average salary is constrained to be increasing in father's education, and the solid red lines connect unconstrained estimates of average salary.", fig.align='center', echo=FALSE, out.width="100%"----
+knitr::include_graphics("figures/daded.png")
+
+
+## ----comppv, echo=FALSE, results='asis'---------------------------------------
+library(knitr)
+library(kableExtra)
+
+years <- c("2008-09","2010-11","2012-13","2014-15","2016-17","2018-19")
+vals <- matrix(
+  c(".008", "n/a",
+    "<.001", ".018",
+    "<.001", "<.001",
+    "<.001", "n/a",
+    ".003", ".417",
+    "<.001", "n/a"),
+  nrow = 1, byrow = TRUE
+)
+df <- as.data.frame(vals, stringsAsFactors = FALSE)
+colnames(df) <- rep(c("one", "two"), length(years))
+
+kable(df, booktabs = TRUE,
+      caption = "One-sided and two-sided $p$-values for the test of the null hypothesis that salary is constant in father's education level. The two-sided test results in n/a when the grid has at least one empty domain.",
+      escape = TRUE) %>%
+  add_header_above(setNames(rep(2, length(years)), years)) %>%
+  kable_styling(latex_options = c("hold_position"))
+
+
+## -----------------------------------------------------------------------------
+load("./nscg19_2.rda")
+data <- nscg2 |>
+  dplyr::filter(hd_year %in% c(2008, 2009))
+
+rds <- svrepdesign(data = data, repweights = dplyr::select(data, "RW0001":"RW0320"), weights = ~w,
+                  combined.weights = TRUE, mse = TRUE, type = "other",
+                  scale = 1, rscale = 0.05)
+
+set.seed(1)
+ans <- csvy(logSalary ~ incr(daded) * field * region, design = rds, test = TRUE)
+
+
+## ----eval=T-------------------------------------------------------------------
+summary(ans)
+
+
+## -----------------------------------------------------------------------------
+data(nhdat, package = "csurvey")
+dstrat <- svydesign(ids = ~ id,  strata = ~ str, data = nhdat,  weight = ~ wt)
+set.seed(1)
+ans <- csvy(chol ~ incr(age) * incr(wcat) * gender, design = dstrat,
+            family = binomial(link = "logit"), test = TRUE)
+
+
+## -----------------------------------------------------------------------------
+summary(ans)
+
+
+## ----eval=FALSE, echo=FALSE---------------------------------------------------
+# ctl <- list(angle = 0, x1size = 2, x2size = 2, x1lab = "waist", x2_labels = c("male", "female"),
+#   subtitle.size=6)
+# plot(ans, x1 = "wcat", x2 = "gender", type="both", control = ctl)
+
+
+## ----nhanesbin, fig.cap="Estimates of probability of high cholesterol level for each combination of age, waist and gender. The blue dots represent the constrained domain mean estimates, while the green dots represent the unconstrained domain mean estimates. The blue band is the 95% confidence interval for the domains, using the constraints; the green band is the 95% unconstrained domain mean confidence interval.", fig.align='center', echo=FALSE----
+knitr::include_graphics("figures/nhanes_bin.png")
+