[R] Survey-weighted ordered logistic regression

[Thomas Lumley]

On Tue, 18 Jul 2006, Debarchana Ghosh wrote:

> Hi,
>
> I am trying to fit a model with an ordered response variable (3 levels) and
> 13 predictor variables. The sample has complex survey design and I've used
> 'svydesign' command from the survey package to specify the sampling design.
> After reading the manual of 'svyglm' command, I've found that you can fit a
> logistic regression (binary response variable) by specifying the
> family=binomial in svyglm function. However I'm unable to fit an ordered
> logistic model in 'svyglm function'.
>

You can do this most easily using withReplicates().  If you have a survey 
object created with svydesign, use as.svrepdesign() to turn it into a 
replicate-weights object and then do eg

library(MASS) 
ologit <- withReplicates( repwtdesign,
                 quote(coef(polr(y~x1+x2+x3+x4, weights=.weights))))

You will get a lot of (harmless) warnings about "non-integer #successes in 
a binomial glm!".  This will not give you standard errors for the 
intercept terms - if you want them you can do

ologit.int <- withReplicates( repwtdesign,
                 quote(polr(y~x1+x2+x3+x4, weights=.weights)$zeta))


You could also use svymle() and get Taylor linearisation standard errors, 
but that would require writing out the loglikelihood and its gradient, 
which is tedious (although you could borrow most of it from MASS::polr)


 	-thomas

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle