svycontrast

I got asked for more detail about the svycontrast() function, so I thought I’d post it here too. The function is related to the CONTRASTS you get in SAS, but focused on estimation rather than testing.

The input to svycontrast() is a $$p$$-vector of estimates $$\hat\theta$$ (which I’ll consider as a column vector) and an estimated $$p\times p$$ covariance matrix $$\hat\Xi$$

There are two main cases:

Linear

Given a $$p$$-vector of coefficients $$b$$, the function computes $$b^T\hat\theta$$ and $$b^T\hat\Xi b$$. Given a list of $$k$$ $$p$$-vectors of coefficients the function pastes them into a $$p\times k$$ matrix $$B$$ and computes $$B^T\theta$$ and $$B^T\hat\Xi B$$.

For example, from the help page

library(survey)
data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
a <- svymean(~api00+enroll+api99, dclus1)
b<-svycontrast(a, list(avg=c(0.5,0,0.5), diff=c(1,0,-1)))
b
##      contrast      SE
## avg   625.574 23.8362
## diff   37.191  3.0852

and doing it step by step

coef(a)
##    api00   enroll    api99
## 644.1694 549.7158 606.9781
crossprod(c(0.5,0,0.5),coef(a))
##          [,1]
## [1,] 625.5738
crossprod(c(1,0,-1),coef(a))
##          [,1]
## [1,] 37.19126
sqrt(crossprod(c(0.5,0,0.5), vcov(a))%*%c(0.5,0,0.5))
##          [,1]
## [1,] 23.83622
sqrt(crossprod(c(1,0,-1), vcov(a))%*%c(1,0,-1))
##          [,1]
## [1,] 3.085197

These might be rounded differently when they print, but are the same,eg

crossprod(c(0.5,0,0.5),coef(a))-coef(b)[1]
##      [,1]
## [1,]    0
sqrt(crossprod(c(1,0,-1), vcov(a))%*%c(1,0,-1)) - SE(b)[2]
##      [,1]
## [1,]    0

And the covariance term:

crossprod(c(0.5,0,0.5), vcov(a))%*%c(1,0,-1)
##           [,1]
## [1,] -16.30776
vcov(b)[1,2]
## [1] -16.30776

You can also use names to indicate coefficients, so that this is the same

svycontrast(a, list(avg=c(api00=0.5,api99=0.5), diff=c(api00=1,api99=-1)))
##      contrast      SE
## avg   625.574 23.8362
## diff   37.191  3.0852

Non-linear

Given a quoted expression in which the free variables are the names of the coefficients, svycontrast() treats it as a function $$f()$$ and computes $$f(\hat \theta)$$ and $$f'(\hat\theta)^T\hat\Xi f'(\hat\theta)$$, using deriv() to do the symbolic differentiation.

As a trivial case, you can, of course, do linear combinations this way and get the same as above

svycontrast(a, list(quote(api00/2+api99/2), quote(api00-api99)))
##        nlcon      SE
## [1,] 625.574 23.8362
## [2,]  37.191  3.0852

Less trivially: geometric means, where $$f(\theta)=\exp\theta$$ and so $$f'(\hat\theta)$$ is a diagonal $$2\times 2$$ matrix with diagonal entries $$\exp(\hat\theta)$$

meanlogs <- svymean(~log(api00)+log(api99), dclus1)
meanlogs
##              mean     SE
## log(api00) 6.4541 0.0378
## log(api99) 6.3905 0.0411
geomeans<-svycontrast(meanlogs,
list(api00=quote(exp(log(api00))), api99=quote(exp(log(api99)))))
exp(coef(meanlogs))-coef(geomeans)
## log(api00) log(api99)
##          0          0
B<-diag(exp(coef(meanlogs)))
crossprod(B, vcov(meanlogs))%*%B - vcov(geomeans)
##       api00 api99
## api00     0     0
## api99     0     0

Even less trivially, svykappa() does a bit of quasiquotation as well

svycontrast(probs, list(kappa = bquote((.(obs) - .(expect))/(1 - .(expect)))))

where obs are the names of diagonal entries of a square contingency table and expect are the products that define the expected values under independency. For a $$2\times 2$$ table the expression comes out as

(.a_..A_ + .b_..B_ - (.a_ * .A_ + .b_ * .B_))/(1 - (.a_ * .A_ +
.b_ * .B_))

and svycontrast() uses the deriv() function to differentiate that expression with respect to all its entries.