[R] svyglm and sandwich estimator of variance

2008-12-19 Thread Roberta Pereira Niquini
Hi, 

I would like to estimate coefficients using poisson regression and then get 
standard errors that are adjusted for heteroskedasticity, using a complex 
sample survey data. Then I will calculate prevalence ratio and confidence 
intervals.
Can sandwich estimator of variance be used when observations aren’t 
independent? In my case, observations are independent across groups 
(clusters), but not necessarily within groups.  Can I calculate the standard 
errors with robust variance, in complex sample survey data using R? 

Outputs:

design_tarv<-svydesign(ids=~X2, strata=~X3, data=banco, weights=~X4)

banco.glm7 <- svyglm(y ~x1, data = banco,  family = poisson (link= "log"), 
design= design_tarv)
summary(banco.glm7)

Call:
svyglm(y ~ x1, data = banco, family = poisson(link = "log"), 
design = design_tarv)

Survey design:
svydesign(ids = ~X2, strata = ~X3, data = banco, 
weights = ~X4)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.918930.04696 -19.570  < 2e-16 ***
x1  0.197100.06568   3.001  0.00603 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 0.5722583)

Number of Fisher Scoring iterations: 5


library(sandwich)

vcovHC(banco.glm7) 
 (Intercept)x1
(Intercept)  4.806945e-13  -4.771409e-13
x1 -4.771409e-137.127168e-13

sqrt(diag(vcovHC(banco.glm7, type="HC0"))) 
(Intercept)   x1 
6.923295e-078.426314e-07

# I think this result isn’t correct, because standard errors are so small.


Thank you for the help, 
Roberta Niquini.







--
ENSP - Fiocruz

__
R-help@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.


[R] prevalence ratio and confidence intervals

2008-12-12 Thread Roberta Pereira Niquini
Hi everybody,

I would like to estimate prevalence ratio and confidence intervals. 

I tried to do a log-binomial regression, but there was a failure of 
convergence.
Now, I would like to learn how to do a poisson regression with robust 
variance.
I am trying to estimate coefficients with poisson regression and then get 
standard errors that are adjusted for heteroskedasticity. 

glm22<- svyglm(y~x1+x2+x3+offset(log(x4)), data = banco,  family = poisson, 
design= design_tarv)

# Y has a binomial distribution (0/1)
# X1, X2, X3 e X4 are categorical variables.
#I am using the library(survey) because it is an analysis of Complex Sample 
Survey Data .

summary(glm22)

Call:
svyglm(y~x1+x2+x3+ offset(log(x4)),data = banco, family = poisson, design = 
design_tarv)

Survey design:
svydesign(ids = ~conglomerado, strata = ~estrato, data = banco, 
weights = ~peso)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.612240.07223 -77.699  < 2e-16 ***
x1   0.338470.07428   4.557 0.000155 ***
x2   0.177450.07059   2.514 0.019765 *  
x3   0.335080.09447   3.547 0.001808 ** 
x4   0.243820.08808   2.768 0.011217 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 0.7535822)

Number of Fisher Scoring iterations: 5

# Using family=quasipoisson, I found the same values.


library(sandwich)

vcovHAC(glm22)

 (Intercept)x1   x2   x3x4
(Intercept)1.060857e-12-1.306035e-13-5.139155e-13 -9.788354e-13 -3.428080e-13
x1 -1.306035e-13  7.237868e-13   -3.263182e-13  -1.620593e-13  1.704392e-13
x2 -5.139155e-13  -3.263182e-13  1.250564e-12   7.207572e-13   -9.350062e-13
x3 -9.788354e-13  -1.620593e-13  7.207572e-13   1.707176e-12   -2.244859e-13
x4 -3.428080e-13   1.704392e-13   -9.350062e-13  -2.244859e-13   2.031640e-12

sqrt(diag(vcovHAC(glm22)))

 (Intercept)   x1x2x3 x4
1.029979e-06 8.507566e-07 1.118286e-06 1.306589e-06 1.425356e-06 


I think these standards errors are very small. 

Is this the correct form to do poisson regression with robust variance?

Thank you for the help, 
Roberta.

__
R-help@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.