Publish 2025-4 issue

_articles/RJ-2014-031/RJ-2014-031.Rmd

@@ -18,7 +18,7 @@ abstract: Assessing the assumption of multivariate normality is required by many
   robust Mahalanobis distances. Moreover, this package offers functions to check the
   univariate normality of marginal distributions through both tests and plots. Furthermore,
   especially for non-R users, we provide a user-friendly web application of the package.
-  This application is available at <http://www.biosoft.hacettepe.edu.tr/MVN/>.
+  This application is available at <https://biosoft.shinyapps.io/mvn-shiny-app/> (link updated 27 Jan 2026).
 author:
 - name: Selcuk Korkmaz
   affiliation: Department of Biostatistics, Hacettepe University
@@ -190,7 +190,7 @@ $\mu$ and variance $\sigma^2$ as given below (updated 2025-08-29: $\mu$ and $\si
 
 $$\begin{aligned}
   \mu       = & \, 1 - a^{-p/2}\left(1 + \frac{p\,\beta^{2}}{a} + \frac{p(p+2)\,\beta^{4}}{2a^{2}}\right)\\
-  \sigma^2  = & \, 2(1+4\beta^{2})^{-p/2} + 2a^{-p}\!\left(1 + \frac{2p\beta^{4}}{a^{2}} + \frac{3p(p+2)\beta^{8}}{4a^{4}}\right)\\ 
+  \sigma^2  = & \, 2(1+4\beta^{2})^{-p/2} + 2a^{-p}\!\left(1 + \frac{2p\beta^{4}}{a^{2}} + \frac{3p(p+2)\beta^{8}}{4a^{4}}\right)\\
               & \, - 4w_{\beta}^{-p/2}\!\left(1 + \frac{3p\beta^{4}}{2w_{\beta}} + \frac{p(p+2)\beta^{8}}{2w_{\beta}^{2}}\right)
 \end{aligned}$$
 
@@ -226,15 +226,15 @@ obtained from the normality transformation proposed by
 $$\label{eq:xandwjs}
 \begin{array}{l l l l}
   \text{if } \: 4\leq n \leq 11; & \qquad x = n \quad & \text{and} & \quad w_j = -\text{log}\left[\gamma - \text{log}\left(1- W_j\right)\right] \\
- \text{if } \: 12\leq n \leq 2000; & \qquad x = \text{log} (n) \quad & \text{and} & \quad w_j = \text{log}\left(1- W_j\right) 
+ \text{if } \: 12\leq n \leq 2000; & \qquad x = \text{log} (n) \quad & \text{and} & \quad w_j = \text{log}\left(1- W_j\right)
 \end{array}  (\#eq:xandwjs)  $$
 
 As seen from equation \@ref(eq:xandwjs), $x$ and $w_j$-s change with the
 sample size ($n$). By using equation \@ref(eq:xandwjs), transformed
 values of each random variable can be obtained from equation
 \@ref(eq:zjs).
 
-$$\label{eq:zjs}  
+$$\label{eq:zjs}
   Z_j = \frac{w_j - \mu}{\sigma}  (\#eq:zjs)  $$
 
 where $\gamma$, $\mu$ and $\sigma$ are derived from the polynomial
@@ -243,10 +243,10 @@ coefficients are provided by @royston1992approx for different sample
 sizes.
 
 $$\begin{aligned}
- \label{eq:polygms} 
+ \label{eq:polygms}
   \gamma =  & \: a_{0\gamma} + a_{1\gamma}x + a_{2\gamma}x^2 + \cdots + a_{d\gamma}x^{d} \nonumber \\
   \mu =     & \: a_{0\mu} + a_{1\mu}x + a_{2\mu}x^2 + \cdots + a_{d\mu}x^{d} \\
-  \text{log}(\sigma) =  & \: a_{0\sigma} + a_{1\sigma}x + a_{2\sigma}x^2 + \cdots + a_{d\sigma}x^{d} \nonumber  
+  \text{log}(\sigma) =  & \: a_{0\sigma} + a_{1\sigma}x + a_{2\sigma}x^2 + \cdots + a_{d\sigma}x^{d} \nonumber
 \end{aligned}  (\#eq:polygms)  $$
 
 The Royston's test statistic for multivariate normality is as follows:
@@ -259,9 +259,9 @@ is the cumulative distribution function for the standard normal
 distribution such that,
 
 $$\begin{aligned}
-  \label{eq:edf}  
+  \label{eq:edf}
   e       &= p / [1 + (p - 1)\bar{c}] \nonumber \\
-  \psi_j   &= \left\{\Phi^{-1}\left[\Phi(-Z_j)/2\right]\right\}^2, \qquad j=1,2,\ldots,p.  
+  \psi_j   &= \left\{\Phi^{-1}\left[\Phi(-Z_j)/2\right]\right\}^2, \qquad j=1,2,\ldots,p.
 \end{aligned}  (\#eq:edf)  $$
 
 As seen from equation \@ref(eq:edf), another extra term $\bar{c}$ has to
@@ -314,7 +314,7 @@ Similarly, the *Iris* data can be loaded from the R database by using
 the following R code:
 
 ``` r
-# load Iris data 
+# load Iris data
 data(iris)
 ```
 
@@ -346,23 +346,23 @@ result
 ```
 
 ``` r
-   Mardia's Multivariate Normality Test 
---------------------------------------- 
-   data : setosa 
+   Mardia's Multivariate Normality Test
+---------------------------------------
+   data : setosa
 
-   g1p            : 3.079721 
-   chi.skew       : 25.66434 
-   p.value.skew   : 0.1771859 
+   g1p            : 3.079721
+   chi.skew       : 25.66434
+   p.value.skew   : 0.1771859
 
-   g2p            : 26.53766 
-   z.kurtosis     : 1.294992 
-   p.value.kurt   : 0.1953229 
+   g2p            : 26.53766
+   z.kurtosis     : 1.294992
+   p.value.kurt   : 0.1953229
 
-   chi.small.skew : 27.85973 
-   p.value.small  : 0.1127617 
+   chi.small.skew : 27.85973
+   p.value.small  : 0.1127617
 
-   Result         : Data are multivariate normal. 
---------------------------------------- 
+   Result         : Data are multivariate normal.
+---------------------------------------
 ```
 
 Here:
@@ -397,15 +397,15 @@ result
 ```
 
 ``` r
-  Henze-Zirkler's Multivariate Normality Test 
---------------------------------------------- 
-  data : setosa 
+  Henze-Zirkler's Multivariate Normality Test
+---------------------------------------------
+  data : setosa
 
-  HZ      : 0.9488453 
-  p-value : 0.04995356 
+  HZ      : 0.9488453
+  p-value : 0.04995356
 
-  Result  : Data are not multivariate normal. 
---------------------------------------------- 
+  Result  : Data are not multivariate normal.
+---------------------------------------------
 ```
 
 Here, `HZ` is the value of the Henze-Zirkler's test statistic at
@@ -429,15 +429,15 @@ result
 ```
 
 ``` r
-  Royston's Multivariate Normality Test 
---------------------------------------------- 
-  data : setosa 
+  Royston's Multivariate Normality Test
+---------------------------------------------
+  data : setosa
 
-  H       : 31.51803 
-  p-value : 2.187653e-06 
+  H       : 31.51803
+  p-value : 2.187653e-06
 
-  Result  : Data are not multivariate normal. 
---------------------------------------------- 
+  Result  : Data are not multivariate normal.
+---------------------------------------------
 ```
 
 Here, `H` is the value of the Royston's test statistic at significance
@@ -519,7 +519,7 @@ following code chunk is used to perform the Shapiro-Wilk's normality
 test on each variable:
 
 ``` r
-uniNorm(setosa, type = "SW") 
+uniNorm(setosa, type = "SW")
 ```
 
 ``` r
@@ -542,7 +542,7 @@ are given in Table [1](#tbl:setosa).
   ------------------------------------------
   Test              Test Statistic   p-value
   --------------- ---------------- ---------
-  Mardia                           
+  Mardia
 
   Skewness                  11.249     0.338
 
@@ -626,7 +626,7 @@ statistical significance of bivariate normal distribution of the
   ------------------------------------------
   Test              Test Statistic   p-value
   --------------- ---------------- ---------
-  Mardia                           
+  Mardia
 
   Skewness                   0.760     0.944
 
@@ -692,11 +692,11 @@ outliers. For this example, we will use another subset of the *Iris*
 data, which is *versicolor* flowers, with the first three variables.
 
 ``` r
-versicolor <- iris[51:100, 1:3] 
+versicolor <- iris[51:100, 1:3]
 # Mahalanobis distance
-result <- mvOutlier(versicolor, qqplot = TRUE, method = "quan") 
+result <- mvOutlier(versicolor, qqplot = TRUE, method = "quan")
 # Adjusted Mahalanobis distance
-result <- mvOutlier(versicolor, qqplot = TRUE, method = "adj.quan") 
+result <- mvOutlier(versicolor, qqplot = TRUE, method = "adj.quan")
 ```
 
 <figure id="fig:mvOutlier">
@@ -815,4 +815,4 @@ Research Fund of Marmara University \[FEN-C-DRP-120613-0273\].
 
 [^1]: <http://www.rstudio.com/shiny/>
 
-[^2]: <http://www.biosoft.hacettepe.edu.tr/MVN>
+[^2]: <https://biosoft.shinyapps.io/mvn-shiny-app/> (link updated 27 Jan 2026)

_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")
+