Estimation problem: 'division by 0' about samplics HOT 7 CLOSED

Hi,

I updated samplics and modified your code to work. Basically, I added the parameter singletons which indicate the option for dealing with the singletons. The options are guided by an ENUM class called SinglePSUEst; see parameter single_psu. The Enum takes these four options: error, skip, certainty, and combine.

class SinglePSUEst(Enum):
    """Estimation options for strata with singleton PSU"""

    error = "Raise Error when one PSU in a stratum"
    skip = "Set variance to zero and skip stratum with one PSU"
    certainty = "Use SSUs or lowest units to estimate the variance"
    combine = "Combine the strata with the singleton psu to another stratum"

The way to use the four options are the following:

• single_psu=SinglePSUEst.error: Let it crash. This good for the first run to identify your singletons
• single_psu=SinglePSUEst.skip: skip the singletons and set the variance to 0.
• single_psu=SinglePSUEst.certainty: treat the singletons as certainties and use SSU (if provided) or the individual records to estimate the variance
• single_psu=SinglePSUEst.combine: combine the singletons strata to other strata. You will have to specify strata_comb. strata_comb is a dictionary to map the old strata to the new strata i.e. {old_stratum1: new_stratum1, old_stratum2: new_straum2, ...}

Given the above, here is your modified code (I used "skip" but your should be able to use the other options)

import pandas as pd
from samplics.estimation import ReplicateEstimator, TaylorEstimator
from samplics.utils import SinglePSUEst

IPM3 = pd.read_csv("IPM3.txt")


def Mean(IPM, var, singletons):  # IPM debe contener los pesos
    IPM = IPM.reset_index()
    mean = TaylorEstimator("mean")
    mean.estimate(
        y=IPM[var],
        samp_weight=IPM["factor07"],
        stratum=IPM["estrato"],
        psu=IPM["conglome"],
        single_psu=singletons,
    )
    return pd.Series([mean.point_est, mean.stderror], index=["Mean", "SE"])



# Run 1
IPM3.groupby("dpto").apply(Mean, "mieperho", SinglePSUEst.skip)

# Run 2: Error when we only use cases that have values different from 0 in "IPM_censored"
IPM3.loc[IPM3["IPM_censored"] != 0].groupby("dpto").apply(Mean, "mieperho", SinglePSUEst.skip)

With regards

from samplics.

In addition, I run the estimation in STATA, with svyset command and default settings and works fine.

from samplics.

Hi

Thank you for reporting this bug. My suspicion is that the number of PSUs in the stratum is one.
It is possible to share a minimum reproducible example that is code+data?
You can send an email if that is easier [email protected].

Best regards

from samplics.

Try with this file:
IPM3 - copia.txt

using this code:

import samplics
from samplics.estimation import TaylorEstimator, ReplicateEstimator
import pandas as pd

IPM3 = pd.read_csv('file.txt')


def Mean(IPM, var): # IPM debe contener los pesos
    IPM = IPM.reset_index()
    mean = TaylorEstimator("mean")
    mean.estimate(
        y=IPM[var],
        samp_weight=IPM["factor07"],
        stratum=IPM["estrato"],
        psu=IPM["conglome"])
    return pd.Series([mean.point_est, mean.stderror], index=['Mean','SE'])

# no Error
IPM3.groupby('dpto').apply(Mean, "mieperho")

# Error when we only use cases that have values different from 0 in "IPM_censored"
IPM3.loc[IPM3["IPM_censored"]!=0].groupby('dpto').apply(Mean, "mieperho")

from samplics.

the problem may be in

when dividing by n-1

from samplics.

Yes, if you have only one PSU then variance cannot be calculated for that stratum that is n-1=0.
So the approach can be either drop that stratum or use a lower level to compute the variance.
I will look at solutions from other software and make a decision for samplics.

from samplics.

Thanks! this may help:
https://www.stata.com/manuals13/svysvyset.pdf (at the end of page 4)

from samplics.