1
$\begingroup$

Does anyone more experienced with R's survey package know about an alternative for svychisq for the case of a contingency table with at least one cell with a zero value?

There does not seem to exist an equivalent of fisher.test for complex survey data.

Please find below the link to the database and the code that I have used.

https://drive.google.com/file/d/1kBKqX1hS3RZuMFw1Q7LgI2bpiw_fF7SB/view?usp=sharing

library(survey)

elsi <- read.csv("ELSI-Portugues.csv", header = T)


elsi$religiao<- elsi$s28
elsi$religiao[elsi$s28 == 99] <- NA
elsi$religiao[elsi$s28 == 88] <- NA 
elsi$religiao<-factor(elsi$religiao,levels = c(1,2,3,4,5,6,7,8,9),labels =
                        c("sem religiao",
                          "catolica",
                          "protestante",
                          "evangelica",
                          "espirita/kardecista",
                          "budista",
                          "islamita",
                          "africanas",
                          "outra")) 
elsi$religiao <- droplevels(elsi$religiao)

elsi$age<-rep(NA,length(elsi$id))
elsi$age[elsi$idade>=60 & elsi$idade<70]<- 1
elsi$age[elsi$idade>=70 & elsi$idade<80]<- 2
elsi$age[elsi$idade>=80 & elsi$idade<90]<- 3
elsi$age[elsi$idade>=90]<- 4
elsi$age<-factor(elsi$age, levels = c(1,2,3,4), labels = c("60 a 69","70 a 79","80 a 89","90 ou +"))

svy_elsiTotal <- svydesign(id=~UPA, weights = ~peso_calibrado_n, strata = ~estrato, data = elsi)

svy_elsi <- subset(svy_elsiTotal, idade >= 60)


> svytable(~age+religiao, svy_elsi)
         religiao
age       sem religiao     catolica  protestante   evangelica espirita/kardecista      budista     islamita
  60 a 69  105.6011348 1722.6946604   14.9541722  553.7651576         103.7053472    3.2060205    0.0000000
  70 a 79   52.1208415  841.3150983    7.9219624  286.0446580          35.3520821    1.3206040    0.0000000
  80 a 89   19.1663548  276.7195582    2.9866518   69.7200985          14.6224819    0.0000000    0.0000000
  90 ou +    0.9560516   26.6008151    0.0000000    8.1345739           0.8632814    0.0000000    0.0000000
         religiao
age          africanas        outra
  60 a 69    5.2064232   38.7326433
  70 a 79    1.8645671   17.4661145
  80 a 89    0.0000000    6.2003879
  90 ou +    0.0000000    0.0000000

> svychisq(~age+religiao, svy_elsi)

Error in solve.default(denom, numr) : 
  system is computationally singular: reciprocal condition number = 5.90227e-18

```
$\endgroup$
4
  • $\begingroup$ I'm presuming you mean a cell with a zero? Is this a structural zero (there's no way this cell could have an observation) or a sampling zero (presumably there are observations here but the sample wasn't large enough to pick up any observations? $\endgroup$
    – num_39
    Apr 12, 2022 at 18:09
  • $\begingroup$ @num_39, thanks for your prompt response. I mean cells with sampling zeros. I have just found this article published in Biometrics in 2015 (ncbi.nlm.nih.gov/pmc/articles/PMC4567525) by Stuart Lipsitz and colleagues, who developed an alternative to the chi-square test for contingency tables with cells with zeros. In the online supplement to that article the authors show the code for a SAS macro that performs that method. Unfortunately, I was not able to find any R package implementing it. Does anyone know of a package implementing that method? $\endgroup$
    – user156625
    Apr 12, 2022 at 19:53
  • $\begingroup$ I didn't look at this article closely, but svychisq also implements two Wald-type tests. Would either of these be sufficient for your purposes? Thomas Lumley is also quite active here so maybe he'll have a suggestion. $\endgroup$
    – num_39
    Apr 12, 2022 at 21:44
  • $\begingroup$ @num_39, I had already checked all the statistics options in svychisq (“F”, “Chisq”, “Wald”, “adjWald”, “lincom”, and “saddlepoint”). Unfortunately, none of them works. The error message that I keep receiving is this one :"Error in solve.default(denom, numr) : system is computationally singular: reciprocal condition number = 5.70074e-18". I checked and this error only happens when there are cells in the contingency table with zeros. If you happen to know Thomas Lumley, please let him know about my request for help. I searched his book on Complex Surveys and could not find any answer. $\endgroup$
    – user156625
    Apr 12, 2022 at 22:45

3 Answers 3

4
$\begingroup$

Ok, with the new question there is less hope. The problem is the very sparse data. In the 90+ age group there are five religion categories that have only 2 observations between them. Worse, the only buddista are in stratum 4 and there are no africanas in stratum 1. (And there are no islamita, but we can delete that category)

The problem isn't individual zero cells, which would make a full loglinear interaction model have poorly-defined variances. The problem is that even the main-effects model, describing the null hypothesis of independence has a singular variance-covariance matrix -- it's an 11-dimensional matrix but it only has rank 3.

I think you need to collapse some of your categories to get a less sparse table and give up on testing for the whole table. None of the usual approaches will work here and I don't know of any way to get a formal test.

I will note, though, that if you rescale the weights to sum to the sample size, the Pearson $X^2$ statistic computed on the weighted table is 23.7, which is smaller than its degrees of freedom. You'd expect a naive $X^2$ statistic like this to overstate the evidence against independence, so I think it's reasonable to say there isn't any real evidence against independence -- though there might be if you collapsed the categories a bit.

$\endgroup$
3
  • $\begingroup$ many thanks for your response! I wouldn't know how to rescale the weights to sum to the sample size. Yet, I found a possible solution for that problem, that did not involve collapsing some categories. That tentative solution involves running a simple logistic regression with svyglm and then using the regTermTest function to evaluate hypothesis of independence. In the example above, that approach did not incur in the singular matrix error. Please let me know if you believe this is an appropriate approach. If you do, I can post it as a complementary answer. $\endgroup$
    – user156625
    Apr 19, 2022 at 23:39
  • $\begingroup$ weights<-weights/mean(weights) rescales them to sum to the sample size $\endgroup$ Apr 20, 2022 at 22:15
  • $\begingroup$ Dr. Lumley, many thanks for your continued support. I have followed your advice and rescaled the weights from the survey design object and rerun the analyses. Unfortunately, svychisq yielded the same error due to the singular variance-covariance matrix problem that you had previously explained. I will describe in a tentative answer below how I used regTermTest to test the hypothesis of independence even without collapsing categories. $\endgroup$
    – user156625
    Apr 24, 2022 at 10:11
3
$\begingroup$

Ok, so the problem needs a bit more specification, because svychisq already works with zero cells

> svytable(~sch.wide+comp.imp, dclus1)
        comp.imp
sch.wide        No       Yes
     No   778.4809    0.0000
     Yes  913.8689 4501.6505
> svychisq(~sch.wide+comp.imp, dclus1)

    Pearson's X^2: Rao & Scott adjustment

data:  svychisq(~sch.wide + comp.imp, dclus1)
F = 236.89, ndf = 1, ddf = 14, p-value = 3.618e-10

> svychisq(~sch.wide+comp.imp, dclus1,statistic="Wald")

    Design-based Wald test of association

data:  svychisq(~sch.wide + comp.imp, dclus1, statistic = "Wald")
F = 18.253, ndf = 1, ddf = 14, p-value = 0.0007739

> svychisq(~sch.wide+comp.imp, dclus1,statistic="adjWald")

    Design-based Wald test of association

data:  svychisq(~sch.wide + comp.imp, dclus1, statistic = "adjWald")
F = 18.253, ndf = 1, ddf = 14, p-value = 0.0007739

> svychisq(~sch.wide+comp.imp, dclus1,statistic="Chisq")

    Pearson's X^2: Rao & Scott adjustment

data:  svychisq(~sch.wide + comp.imp, dclus1, statistic = "Chisq")
X-squared = 69.975, df = 1, p-value < 2.2e-16

I checked with Stata as well, and with an older version of survey in R (3.36, on R 3.6.0) and the same is true there so it isn't just an issue with recent R changes. The tests based on Rao-Scott corrections to the Pearson score statistic or on the Wald statistic in a linear model will work in the presence of zeros. The ones based on likelihood ratio tests won't, but svychisq doesn't do them.

The score test of Lipsitz et al, once I had implemented it (closer to) correctly, turns out to be very close to the adjusted Wald test that's already implemented. Which is as it should be.

Can you show what you actually did and what results you got?

$\endgroup$
2
  • $\begingroup$ (The somewhat related problem I was recalling that isn't solved is confidence intervals for a proportion $p$ when $\hat p=0$) $\endgroup$ Apr 16, 2022 at 4:18
  • $\begingroup$ many thanks for your response! I have edited my question by providing the database, the code that I have used and the error message. Please let me know if you need more details. $\endgroup$
    – user156625
    Apr 16, 2022 at 21:25
0
$\begingroup$

Here goes my tentative solution to this problem. Instead of using svychisq, I used svyloglin followed by regTermTest . I had thought that svyloglin would not work for the same reason why svychisq had not worked, however it looks like it worked.

Please, criticize it if you have reason to believe that it is wrong by any means.

Link to the database and the code that I have used:

https://drive.google.com/file/d/1kBKqX1hS3RZuMFw1Q7LgI2bpiw_fF7SB/view?usp=sharing

library(survey)

elsi <- read.csv("ELSI-Portugues.csv", header = T) # importing the dataset


elsi$religiao<- elsi$s28
elsi$religiao[elsi$s28 == 99] <- NA
elsi$religiao[elsi$s28 == 88] <- NA 
elsi$religiao<-factor(elsi$religiao,levels = c(1,2,3,4,5,6,7,8,9),labels =
                        c("sem religiao",
                          "catolica",
                          "protestante",
                          "evangelica",
                          "espirita/kardecista",
                          "budista",
                          "islamita",
                          "africanas",
                          "outra"))  # formating the “religiao” variable
elsi$religiao <- droplevels(elsi$religiao) # dropping the “islamita” empty category

# Creating a variable for different age strata
elsi$age<-rep(NA,length(elsi$id))
elsi$age[elsi$idade>=60 & elsi$idade<70]<- 1
elsi$age[elsi$idade>=70 & elsi$idade<80]<- 2
elsi$age[elsi$idade>=80 & elsi$idade<90]<- 3
elsi$age[elsi$idade>=90]<- 4
elsi$age<-factor(elsi$age, levels = c(1,2,3,4), labels = c("60 a 69","70 a 79","80 a 89","90 ou +"))

# Updating the survey design object 
svy_elsiTotal <- svydesign(id=~UPA, weights = ~peso_calibrado_n, strata = ~estrato, data = elsi)

# Restricting the survey design object to people aged 60 years and older
svy_elsi <- subset(svy_elsiTotal, idade >= 60)

# Running svychisq, which doesn’t work.

> svychisq(~age+religiao, svy_elsi)

Error in solve.default(denom, numr) : 
  system is computationally singular: reciprocal condition number = 5.90227e-18

# Running a loglinear  model

modrelig<-svyloglin(~age+religiao,svy_elsi)
regTermTest(modrelig, "religiao")

Wald test for religiao
 in svyloglin(~age + religiao, svy_elsi)
F =  242.6663  on  7  and  255  df: p= < 2.22e-16

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.