--- title: "Description of ciuupi" author: "Paul Kabaila, Rheanna Mainzer and Ayesha Perera" date: "13 June 2024" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Description of ciuupi} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Description of the confidence interval that utilizes uncertain prior information (CIUUPI) Suppose that $$y = X \beta + \varepsilon,$$ where $y$ is a random $n$-vector of responses, $X$ is a known $n \times p$ matrix with linearly independent columns, $\beta$ is an unknown parameter $p$-vector, and $\varepsilon \sim N(0, \, \sigma^2 \, I)$, with $\sigma^2$ assumed known. Suppose that the parameter of interest is $\theta = a^{\top} \beta$. Let $\tau = c^{\top} \beta - t = 0$, where $a$ and $c$ are specified linearly independent nonzero $p$-vectors and $t$ is a specified number. Suppose that we have uncertain prior information that $\tau = 0$. Define the scaled expected length of a confidence interval to be $$\frac{\text{expected length of this confidence interval}} {\text{expected length of the standard } 1-\alpha \text{ confidence interval for }\theta}.$$ This package computes a confidence interval, with minimum coverage $1 - \alpha$, for $\theta$ that utilizes the uncertain prior information that $\tau = 0$ in the following sense. This confidence interval has scaled expected length that (a) is substantially less 1 when $\tau=0$ and (b) has maximum value that is not too large. In addition, this confidence interval approaches the standard $1-\alpha$ confidence interval for $\theta$ when the data and the prior information are strongly discordant. Let $\widehat{\beta}$ denote the least squares estimator of $\beta$. Then $\widehat{\theta} = a^{\top} \widehat{\beta}$ and $\widehat{\tau} = c^{\top} \widehat{\beta} - t$ are the least squares estimators of $\theta$ and $\tau$, respectively. Also let $v_{\theta} = \text{Var}(\widehat{\theta})/\sigma^2 = a^{\top} (X^{\top}X)^{-1}a$ and $v_{\tau} = \text{Var}(\widehat{\tau})/\sigma^2 = c^{\top} (X^{\top}X)^{-1}c$. The known correlation between $\widehat{\theta}$ and $\widehat{\tau}$ is $\rho = a^{\top} (X^{\top}X)^{-1}c / (v_{\theta} \, v_{\tau})^{1/2}$. Let $\gamma = \tau / (\sigma \, v_{\tau}^{1/2})$ and $\hat{\gamma} = \hat{\tau}/(\sigma \, v_{\tau}^{1/2})$. The confidence interval for $\theta$ that utilizes uncertain prior information that $\tau=0$ has the form $$ \text{CI}(b,s) =\left[ \hat{\theta} - v_{\theta}^{1/2} \, \sigma \, b(\hat{\gamma}) - v_{\theta}^{1/2} \, \sigma \, s(\hat{\gamma}), \, \hat{\theta} - v_{\theta}^{1/2} \, \sigma \, b(\hat{\gamma}) + v_{\theta}^{1/2} \, \sigma \, s(\hat{\gamma}) \right], $$ where $b: \mathbb{R} \rightarrow \mathbb{R}$ is an odd continuous function and $s: \mathbb{R} \rightarrow [0, \infty)$ is an even continuous function. In addition, $b(x)=0$ and $s(x)= z_{1-\alpha/2}$ for all $|x| \ge d$, where $z_{1-\alpha/2}$ denotes the $1 - \alpha/2$ quantile of the standard normal distribution. For given $d > 0$ and positive integer $q$, let $h = d / q$ and specify the functions $b$ and $s$ by the $(2q - 1)$-vector $$\Big(b(h),...,b((q-1)h), s(0),s(h)...,s((q-1)h)\Big).$$ The values of $b(kh)$ and $s(kh)$ for $k=-q,\dots,q$ are deduced from this vector using the assumptions made about the functions $b$ and $s$. The values of $b(x)$ and $s(x)$ for any $x \in [-d,d]$ are found via cubic spline interpolation using the values of $b(kh)$ and $s(kh)$ for $k=-q,\dots,q$. This is through natural (default) or clamped cubic spline interpolation. Version 1.0 of ciuupi was described by Mainzer and Kabaila (2019), who chose $d = 6$ and $q = 6$. However, it was subsequently found that this choice is no longer adequate when $\alpha$ is very small and/or $|\rho|$ is close to 1. Extensive numerical explorations have been used to find a formula (in terms of $\alpha$ and $|\rho|$) for a 'goldilocks' value of $d$ that is neither too large nor too small. Version 1.2 uses this formula and sets $q = \lceil d/0.75 \rceil$ and $h = d/q$. The confidence interval that utilizes the uncertain prior information (CIUUPI) is found by computing the value of the vector $\big(b(h),...,b((q-1)h), s(0),s(h)...,s((q-1)h)\big)$, such that the confidence interval $\text{CI}(b,s)$ has minimum coverage probability $1 - \alpha$ and the desired scaled expected length properties. The method used is an obvious extension of that described by Mainzer and Kabaila (2019). However, version 1.2 employs a far more efficient algorithm for the computation of $\lambda^*$, which forms part of the computation of the CIUUPI. This confidence interval has the following three practical applications. Firstly, if $\sigma^2$ has been accurately estimated from previous data then it may be treated as being effectively known. Secondly, for sufficiently large $n - p$ ($n-p \ge 30$, say), if we replace the assumed known value of $\sigma^2$ by its usual estimator in the formula for the confidence interval then the resulting interval has, to a very good approximation, the same coverage probability and expected length properties as when $\sigma^2$ is known. Thirdly, some more complicated models can be approximated by the linear regression model with $\sigma^2$ known when certain unknown parameters are replaced by estimates (see e.g. Kabaila and Mainzer, 2019). ## Functions in this package We attach the package `ciuupi` using ```{r setup} library(ciuupi) ``` We illustrate the application of the functions in the package using the real-life example described in Discussion 5.8 on p.3426 of Kabaila and Giri (2009). This example is obtained by extracting a $2 \times 2$ factorial data set from the $2^3$ factorial data set described in Table 7.5 of Box et al. (1963). The resulting design matrix $X$ is specified by the command ```{r} x1 <- c(-1, 1, -1, 1) x2 <- c(-1, -1, 1, 1) X <- cbind(rep(1, 4), x1, x2, x1*x2) ``` A description of the parameter of interest and the uncertain prior information is given in Discussion 5.8 on p.3426 of Kabaila and Giri (2009). The parameter of interest is $\theta = a^{\top} \beta$, where the column vector $a$ is specified by the command ```{r} a <- c(0, 2, 0, -2) ``` The uncertain prior information is that the two-factor interaction term is zero i.e. $\tau = 0$, where $\tau = c^{\top} \beta$ and the column vector $c$ is specified by the command ```{r} c <- c(0, 0, 0, 1) ``` ### `acX_to_rho` To compute the CIUUPI, we need to first compute the known correlation $\rho$ between $\widehat{\theta}$ and $\widehat{\tau}$, given by $\rho = a^{\top} (X^{\top}X)^{-1}c \big / \big(a^{\top} (X^{\top}X)^{-1}a \, c^{\top} (X^{\top}X)^{-1}c \big)^{1/2}$. This is done using the function `acX_to_rho` as follows. ```{r} # Compute rho rho <- acX_to_rho(a, c, X) print(rho) ``` The exact value of $\rho$ is $-1/\sqrt{2}$. ### Compute the functions $b$ and $s$ that specify the CIUUPI ### `bs_ciuupi` The function `bs_ciuupi` is the "engine" of the package. The CIUUPI is required to have minimum coverage probability $1 - \alpha$, for some $\alpha$ specified by the user of the package. The functions $b$ and $s$ of the CIUUPI are determined by $\alpha$ and $\rho$. Suppose that $\alpha = 0.05$. The functions $b$ and $s$ can be specified either by natural cubic spline interpolation (natural=1) or by clamped cubic spline interpolation (natural=0). By default, natural=1. For this default, the list that specifies the functions $b$ and $s$ of the CIUUPI is computed using the following command bs.list.example <- bs_ciuupi(alpha, rho) This command takes about 5 minutes to run and so `bs.list.example` has been stored as data for quick access. We can specify the functions $b$ and $s$ by clamped cubic spline interpolation as follows. The list that specifies the functions $b$ and $s$ of the CIUUPI is computed using the following command bs_ciuupi(alpha, rho, natural = 0) This command takes about 45 minutes to run and leads to almost identical functions $b$ and $s$ of the CIUUPI as when bs_ciuupi(alpha, rho) is used. The output of `bs_ciuupi` is a list that includes the element bsvec, which is the vector $$\Big(b(h),...,b((q-1)h), s(0),s(h)...,s((q-1)h)\Big)$$ that, using cubic spline interpolation, determines the functions $b$ and $s$ that specify the CIUUPI for all possible values of $\sigma$ and observed response vector $y$. We stress that bsvec does NOT depend on either $\sigma$ or the observed response vector $y$. The graphs of the functions $b$ and $s$ are plotted using the functions `plot_b` and `plot_s`, respectively. ### Plot the graphs of the functions $b$ and $s$ ### `plot_b` and `plot_s` We plot the graph of the odd function $b$ of the CIUUPI using ```{r, fig.align='center', fig.width=5, fig.height=4} plot_b(bs.list.example) ``` We plot of the graph of the even function $s$ of the CIUUPI using ```{r, fig.align='center', fig.width=5, fig.height=4} plot_s(bs.list.example) ``` ### The performance of the CIUUPI Define the scaled expected length of the CIUUPI to be the expected length of the CIUUPI divided by the expected length of the standard $1 - \alpha$ confidence interval for $\theta$, given by $$ \left[ \hat{\theta} - z_{1-\alpha/2} \, v_{\theta}^{1/2} \, \sigma, \, \hat{\theta} + z_{1-\alpha/2} \, v_{\theta}^{1/2} \, \sigma \right]. $$ The performance of the CIUUPI is assessed by its coverage probability and its squared scaled expected length. ### `plot_cp` and `plot_squared_sel` The coverage probability of the CIUUPI is an even function of the unknown parameter $\gamma$. We plot the graph of the coverage probability minus $1 - \alpha$, as a function of $|\gamma|$, using ```{r, fig.align="center", fig.width=5, fig.height=4} plot_cp(bs.list.example) ``` The scaled expected length of the CIUUPI is an even function of the unknown parameter $\gamma$. The squared scaled expected length may be interpreted as the following measure of efficiency of the standard $1-\alpha$ confidence interval for $\theta$ relative to the CIUUPI with minimum coverage $1-\alpha$. For a given large sample size used to find the the standard $1-\alpha$ confidence interval for $\theta$, it is the sample size required by the CIUUPI with minimum coverage $1-\alpha$ divided by the sample size used to find the standard $1-\alpha$ confidence interval for $\theta$ to achieve the same expected length, for the specified value of $\gamma$. We plot the graph of the squared scaled expected length, as a function of $|\gamma|$, using ```{r, fig.align="center", fig.width=5, fig.height=4} plot_squared_sel(bs.list.example) ``` ### Computation of the observed value of the CIUUPI Now we introduce the additional information provided by $\sigma^2$ and the observed response vector $y$. This permits us to compute the observed value of the CIUUPI and also to compute the standard $1 - \alpha$ confidence interval for $\theta$. ### `ciuupi_observed_value` and `ci_standard` For this example, the observed response vector $y$ is (87.2, 88.4, 86.7, 89.2) and the known value of $\sigma$ is 0.8. We specify these values using the following commands ```{r} y <- c(87.2, 88.4, 86.7, 89.2) sig <- 0.8 ``` The uncertain prior information is that $\tau = 0$, where $\tau = c^{\top} \beta - t$, with $t = 0$. The observed value of the CIUUPI is computed and printed using the following commands ```{r} alpha <- 0.05 t <- 0 ciuupi_observed_value(a, c, X, alpha, bs.list.example, t, y, sig = sig) ``` For comparison, the standard $1-\alpha$ confidence interval for $\theta$, with $\alpha = 0.05$, is computed and printed using the following commands ```{r} alpha <- 0.05 ci_standard(a, X, y, alpha, sig) ``` ## References Box, G.E.P., Connor, L.R., Cousins, W.R., Davies, O.L., Hinsworth, F.R., Sillitto, G.P. (1963) The Design and Analysis of Industrial Experiments, 2nd edition, reprinted. Oliver and Boyd, London. Kabaila, P. and Mainzer, R. (2019). An alternative to confidence intervals constructed after a Hausman pretest in panel data. Mainzer, R. and Kabaila, P. (2019). ciuupi: An R package for computing confidence intervals that utilize uncertain prior information.