Search code examples
rsurvey

R survey package function svyciprop with "likelihood" method


I'm trying to calculate confidence intervals of proportions using the R survey package function svyciprop with the "likelihood" method.

Here's some sample code:

df <- data.frame(id =  c(1, 1, 1, 2, 2, 2), var = c("a", "b", "a", "b", "a", "b"))
survey_design <- svydesign(id = ~id, data = df)
svyciprop(~I(var == "a"), survey_design, method = "likelihood")

This generates the error message:

Error in seq.int(xmin, xmax, length.out = n) : 'from' must be finite

I can find nothing in the package documentation that explains how to get this to work.

Many thanks!


Solution

  • the documentation for svyciprop is found by typing ?svyciprop or by googling svyciprop but the documentation won't cover something as specific as your error.

    since all R code is available for users to read, you can debug the function that you are using. debug survey:::svyciprop which leads you to survey:::confint.svyglm which leads you to MASS:::confint.glm which leads you to MASS:::confint.profile.glm and so on. the internet has lots of explanations of how to use the debug function in R. there are a lot of moving parts here

    you're getting some Inf values from the glm objects deep in this calculation that's causing it to break. it's probably related to your example data set being a bit too perfect (and unrealistic). ;)

    if i throw out a single observation from your df then it works.

    library(survey)
    df <- data.frame(id =  c(1, 1, 2, 2, 2), var = c("a", "b", "b", "a", "b"))
    survey_design <- svydesign(id = ~id, data = df)
    svyciprop(~I(var == "a"), survey_design, method = "likelihood")
    

    the example data sets found at the bottom of ?svyciprop also work.

    your particular issue is that the confidence interval you're requesting is impossible.

    # watch how the confidence interval tends toward zero and one as you widen it.
    svyciprop(~I(var == "a"), survey_design, method = "likelihood",level=0.8)
    svyciprop(~I(var == "a"), survey_design, method = "likelihood",level=0.9)
    svyciprop(~I(var == "a"), survey_design, method = "likelihood",level=0.93)
    svyciprop(~I(var == "a"), survey_design, method = "likelihood",level=0.95) # this is the default