In my (limited) experience survey weighting for multilevel models is a can of worms (I mostly follow Andrew Gelman's opinion on this, e.g this paper or this blog post about it; here's an (unanswered) CrossValidated question on the topic). In other words, you might have to solve a statistical problem (what is the right way to handle this?) before you worry about the computational problem.
I do think Stata and/or SAS offer such capabilities (see CrossValidated link above), but I can't tell you much about them.
However, if you don't want/need random effects, i.e. you just want to handle a problem with a large number of fixed effects, you should be able to use a sparse model matrix as in this question, as in the example below.
It will still be up to you to figure out how to deal with the weights. I don't believe that the weights=
argument in lm
, glm
, lme4
, etc. handle survey weights. I believe the survey
package is very sophisticated, but I don't know if its surveyglm()
function can handle/has an option for sparse model matrices ...
library("Matrix")
library("MatrixModels")
set.seed(101)
ngrps <- 1e3
nobs <- 10
ntot <- nobs*ngrps
d <- data.frame(y=rnorm(ntot),
x=rnorm(ntot),
f=gl(ngrps,nobs))
object.size(X1 <- model.matrix(~x+f,data=d)) ## 80 MB
object.size(X2 <- sparse.model.matrix(~x+f,data=d)) ## 709K
system.time(m1 <- MatrixModels:::lm.fit.sparse(X2,d$y)) ## 0.004 seconds
system.time(m2 <- lm.fit(X1,d$y)) ## 9 seconds