Sorry, forgot to "reply all"...
Laust wrote:> Dear list,
>
> I am trying to set up a propensity-weighted regression using the
> survey package. Most of my population is sampled with a sampling
> probability of one (that is, I have the full population). However, for
> a subset of the data I have only a 50% sample of the full population.
> In previous work on the data, I analyzed these data using SAS and
> STATA. In those packages I used a propensity weight of 1/[sampling
> probability] in various generalized linear regression-procedures, but
> I am having trouble setting this up. I bet the solution is simple, but
> I?m a R newbie. Code to illustrate my problem below.
Hi Laust,
You probably need the package author to explain fully, but as far as I
can see, the crux is that a dispersion parameter is being used, based on
Pearson residuals, even in the Poisson case (i.e. you effectively get
the same result as with quasipoisson()).
I don't know what the rationale is for this, but it is clear that with
your data, an estimated dispersion parameter is going to be large. E.g.
the data has both 0 cases in 750000 person-years and 1000 cases in 5000
person-years for Denmark, and in your model they are supposed to have
the same Poisson rate.
summary.svyglm starts off with
est.disp <- TRUE
and AFAICS there is no way it can get set to FALSE. Knowing Thomas,
there is probably a perfectly good reason not to just set the dispersion
to 1, but I don't get it either...
>
> Thanks
> Laust
>
> # loading survey
> library(survey)
>
> # creating data
> listc <-
c("Denmark","Finland","Norway","Sweden","Denmark","Finland","Norway","Sweden")
> listw <- c(1,2,1,1,1,1,1,1)
> listd <- c(0,0,0,0,1000,1000,1000,2000)
> listt <- c(750000,500000,900000,1900000,5000,5000,5000,10000)
> list.cwdt <- c(listc, listw, listd, listt)
> country <-
data.frame(country=listc,weight=listw,deaths=listd,yrs_at_risk=listt)
>
> # running a frequency weighted regression to get the correct point
> estimates for comparison
> glm <- glm(deaths ~ country + offset(log(yrs_at_risk)),
> weights=weight, data=country, family=poisson())
> summary(glm)
> regTermTest(glm, ~ country)
>
> # running survey weighted regression
> svy <- svydesign(~0,,data=country, weight=~weight)
> svyglm <- svyglm(deaths ~ country + offset(log(yrs_at_risk)),
> design=svy, data=country, family=poisson())
> summary(svyglm)
> # point estimates are correct, but standard error is way too large
> regTermTest(svyglm, ~ country)
> # test indicates no country differences
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
O__ ---- Peter Dalgaard ?ster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907