Description Usage Arguments Value Author(s) References Examples
The csvy function performs designbased domain mean estimation with monotonicity and blockmonotone shape constraints.
For example, in a one dimensional situation, we assume that \bar{y}_{U_t} are nondecreasing over T domains. If this monotonicity is not used in estimation, the population domain means can be estimated by the HorvitzThompson estimator or the Hajek estimator. To use the monotonicity information, this csvy function starts from the Hajek estimates \bar{y}_{S_t} = (∑_{k\in S_t}y_k/π_k)/N_t and the isotonic estimator (\hat{θ}_1,…,\hat{θ}_T)^T minimizes the weighted sum of squared deviations from the sample domain means over the set of ordered vectors; that is, \bold{\hat{θ}} is the minimizer of (\tilde{\bold{y}}_{S}  \bold{θ})^T \bold{W}_s (\tilde{\bold{y}}_{S}  \bold{θ}) subject to \bold{Aθ} ≥q \bold{0}, where \bold{W}_S is the diagonal matrix with elements \hat{N}_1/\hat{N},…,\hat{N}_D/\hat{N}, and \hat{N} = ∑_{t=1}^T \hat{N}_t and \bold{A} is a m\times T constraint matrix imposing the monotonicity constraint.
Domains can also be formed from multiple covariates. In that case, a grid will be used to represent the domains. For example, if there are two predictors x_1 and x_2, and x_1 has values on D_1 domains: 1,…,D_1, x_2 has values on D_2 domains: 1,…,D_2, then the domains formed by x_1 and x_2 will be a D_1\times D_2 by 2 grid.
To get 100(1α)\% approximate confidence intervals or surfaces for the domain means, we apply the method in Meyer, M. C. (2018). \hat{p}_J is the estimated probability that the projection of y_s onto \cal C lands on \cal F_J, and the \hat{p}_J values are obtained by simulating many normal random vectors with estimated domain means and covariance matrix I, where I is a M \times M matrix, and recording the resulting sets J.
The user needs to provide a survey design, which is specified by the svydesign function in the survey package, and also a data frame containing the response, predictor(s), domain variable, sampling weights, etc. So far, only stratified sampling design with simple random sampling without replacement (STSI) is considered in the examples in this package.
Note that when there is any empty domain, the user must specify the total number of domains in the nD argument.
1 2 
formula 
A formula object which gives a symbolic description of the model to be fitted. It has the form "response ~ predictor". The response is a vector of length n. For now, the response can only be gaussian. A predictor can be a nonparametrically modelled variable with a monotonicity or convexity restriction, or a combination of both. In terms of a nonparametrically modelled predictor, the user is supposed to indicate the relationship between the domain mean and a predictor x in the following way: Assume that μ is the vector of domain means and x is a predictor:

data 
A data frame, list or environment containing the variables in the model. It must be the same as the data frame used in the survey design. 
design 
A survey design, which must be specified by the svydesign routine in the survey package. 
nD 
The total number of domains. 
family 
A parameter indicating the error distribution and link function to be used in the model. It can be a character string naming a family function or the result of a call to a family function. This is borrowed from the glm routine in the stats package. For now, the only family is gaussian. 
amat 
A k \times M matrix imposing shape constraints in each dimension, where M is the total number of domains. If the user doesn't provide the constraint matrix, a subroutine in the csurvey package will create a constraint matrix according to shape constraints specified in the formula. The default is amat = NULL. 
level 
Confidence level of the approximate confidence surfaces. The default is 0.95. 
n.mix 
The number of simulations used to get the approximate confidence intervals or surfaces. If n.mix = 0, no simulation will be done and the face of the final projection will be used to compute the covariance matrix of the constrained estimate. The default is n.mix = 100L. 
test 
A logical scalar. If test == TRUE, then the pvalue for the test H_0:θ is in V versus H_1:θ is in C is returned. C is the constraint cone of the form \{β: Aβ ≥ 0\}, and V is the null space of A. The default is test = TRUE. 
The output is a list of values used for estimation, inference and visualization.
design 
The survey design used in the model. 
muhat 
Estimated shapeconstrained domain means. 
muhat.un 
Estimated unconstrained domain means. 
lwr 
Approximate lower confidence band or surface for the shapeconstrained domain mean estimate. 
upp 
Approximate upper confidence band or surface for the shapeconstrained domain mean estimate. 
lwru 
Approximate lower confidence band or surface for the unconstrained domain mean estimate. 
uppu 
Approximate upper confidence band or surface for the unconstrained domain mean estimate. 
amat 
The k \times M constraint matrix imposing shape constraints in each dimension, where M is the total number of domains. 
grid 
A M \times p grid, where p is the total number of predictors or dimensions. 
nd 
A vector of sample sizes in all domains. 
Ds 
A vector of the number of domains in each dimension. 
cov.c 
Constrained covariance estimate of domain means. 
cov.un 
Unconstrained covariance estimate of domain means. 
cic 
The cone information criterion proposed in Meyer(2013a). It uses the "null expected degrees of freedom" as a measure of the complexity of the model. See Meyer(2013a) for further details of cic. 
Xiyue Liao
Xu, X. and Meyer, M. C. (2021) Onesided testing of population domain means in surveys.
Oliva, C., Meyer, M. C., and Opsomer, J.D. (2020) Estimation and inference of domain means subject to qualitative constraints. Survey Methodology
Meyer, M. C. (2018) A Framework for Estimation and Inference in Generalized Additive Models with Shape and Order Restrictions. Statistical Science 33(4) 595–614.
Wu, J., Opsomer, J.D., and Meyer, M. C. (2016) Survey estimation of domain means that respect natural orderings. Canadian Journal of Statistics 44(4) 431–444.
Meyer, M. C. (2013a) Semiparametric additive constrained regression. Journal of Nonparametric Statistics 25(3), 715.
Lumley, T. (2004) Analysis of complex survey samples. Journal of Statistical Software 9(1) 1–19.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138  data(api)
mcat=apipop$meals
for(i in 1:10){mcat[trunc(apipop$meals/10)+1==i]=i}
mcat[mcat==100]=10
D1=10
gcat=apipop$col.grad
for(i in 1:10){gcat[trunc(apipop$col.grad/10)+1==i]=i}
gcat[gcat >= 5]=4
D2=4
nsp=c(200,200,200)*1 ## sample sizes per stratum
es=sample(apipop$snum[apipop$stype=='E'&!is.na(apipop$avg.ed)&!is.na(apipop$api00)],nsp[1])
ms=sample(apipop$snum[apipop$stype=='M'&!is.na(apipop$avg.ed)&!is.na(apipop$api00)],nsp[2])
hs=sample(apipop$snum[apipop$stype=='H'&!is.na(apipop$avg.ed)&!is.na(apipop$api00)],nsp[3])
sid=c(es,ms,hs)
pw=1:6194*0+4421/nsp[1]
pw[apipop$stype=='M']=1018/nsp[2]
pw[apipop$stype=='H']=755/nsp[3]
fpc=1:6194*0+4421
fpc[apipop$stype=='M']=1018
fpc[apipop$stype=='H']=755
strsamp=cbind(apipop,mcat,gcat,pw,fpc)[sid,]
dstrat<svydesign(ids=~snum, strata=~stype, fpc=~fpc, data=strsamp, weight=~pw)
rds=as.svrepdesign(dstrat, type="JKn")
# Example 1: monotonic in one dimension
ansc1 = csvy(api00~decr(mcat),data=strsamp,design=dstrat, nD=D1)
# checked estimated domain means
# ansc1$muhat
# checked sorted estimated domain means
# ansc1$muhat.s
# Example 2: monotonic in two dimensions
ansc2 = csvy(api00~incr(gcat)*decr(mcat),data=strsamp,design=dstrat, nD=(D1*D2))
plotpersp(ansc2, ci="up", th=140)
plotpersp(ansc2, th=140)
# Example 3: monotonic in three dimensions
D1 = 5
D2 = 5
D3 = 6
Ds = c(D1, D2, D3)
M = cumprod(Ds)[3]
x1vec = 1:D1
x2vec = 1:D2
x3vec = 1:D3
grid = expand.grid(x1vec, x2vec, x3vec)
N = M*100*4
Ns = rep(N/M, M)
mu.f = function(x) {
mus = x[1]^(0.25)+4*exp(0.5+2*x[2])/(1+exp(0.5+2*x[2]))+sqrt(1/4+x[3])
mus = as.numeric(mus$Var1)
return (mus)
}
mus = mu.f(grid)
H = 4
nh = c(180,360,360,540)
n = sum(nh)
Nh = rep(N/H, H)
#generate population
y = NULL
z = NULL
set.seed(1)
for(i in 1:M){
Ni = Ns[i]
mui = mus[i]
ei = rnorm(Ni, 0, sd=1)
yi = mui + ei
y = c(y, yi)
zi = i/M + rnorm(Ni, mean=0, sd=1)
z = c(z, zi)
}
x1 = rep(grid[,1], times=Ns)
x2 = rep(grid[,2], times=Ns)
x3 = rep(grid[,3], times=Ns)
domain = rep(1:M, times=Ns)
cts = quantile(z, probs=seq(0,1,length=5))
strata = 1:N*0
strata[z >= cts[1] & z < cts[2]] = 1
strata[z >= cts[2] & z < cts[3]] = 2
strata[z >= cts[3] & z < cts[4]] = 3
strata[z >= cts[4] & z <= cts[5]] = 4
freq = rep(N/(length(cts)1), n)
w0 = Nh/nh
w = 1:N*0
w[strata == 1] = w0[1]
w[strata == 2] = w0[2]
w[strata == 3] = w0[3]
w[strata == 4] = w0[4]
pop = data.frame(y = y, x1 = x1, x2 = x2, x3 = x3, domain = domain, strata = strata, w=w)
ssid = stratsample(pop$strata, c("1"=nh[1], "2"=nh[2], "3"=nh[3], "4"=nh[4]))
sample.stsi = pop[ssid, ,drop=FALSE]
ds = svydesign(id=~1, strata =~strata, fpc=~freq, weights=~w, data=sample.stsi)
#domain means are increasing w.r.t x1, x2 and block monotonic in x3
ord = c(1,1,2,2,3,3)
ans = csvy(y~incr(x1)*incr(x2)*block.Ord(x3,order=ord), data=sample.stsi, design=ds, n.mix=0)
#3D plot of estimated domain means: x1 and x2
plotpersp(ans)
#3D plot of estimated domain means: x3 and x2
plotpersp(ans, x3, x2)
#3D plot of estimated domain means: x3 and x2 for each domain of x1
plotpersp(ans, x3, x2, categ="x1")
#3D plot of estimated domain means: x3 and x2 for each domain of x1
plotpersp(ans, x3, x2, categ="x1", NCOL = 3)
# Example 4: unconstrained in one dimension
#no constraint on x1
ans = csvy(y~x1*incr(x2)*incr(x3), data=sample.stsi, design=ds, n.mix=0)
#3D plot of estimated domain means: x1 and x2
plotpersp(ans)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.