Correct se for mean proportion in a survey with varying number of trials

Question

I am analysing fishing surveys using the survey package. My PSU is the person and I want to obtain estimates of the proportion of successful fishing trips. Each person has a variable number of fishing trips (trials) and successes.

I can use svyratio() to calc p.hat = mean(successes)/mean(trials) or I can calculate p.i for each person as #succesess/#trials and use svymean() to get p.hat = mean(p.i)

For both approaches the SE calculation uses the #Persons (replicates/PSU’s), not the total number of trials. I’m wondering if this matters?

Of the two approaches, I would prefer to use the latter approach (where p.hat = mean(p.i)), as this is what has been used on a different (non-survey) data set. However, do I need to calculate an alternative se to that obtained via svymean()? And if so, is it just:

p <- svymean(~Pi, s)[[1]]

se.p <- sqrt(p *(1- p)/N))

? Thank you

Example data/code:

d <- data.frame(cbind(NTrips=round(runif(50, min=1, max=50),0), p=runif(50, min=0, max=1), 
weights=rnorm(50,1000,500)))
d <- cbind(Person=paste("P",seq(1:50),sep="_"),d)
d$NSuccess <- d$NTrips*d$p
d$PropST <-d$NSuccess/d$NTrips


library(survey)
s <- svydesign(~1,data=d, weights=~weights)
N <- nrow(d)  # total number of replicates/PSUs

(pr <- svyratio(~NSuccess,~NTrips,s))
(se.pr <- sqrt(pr[[1]] *(1- pr[[1]])/N))
(pm <- (svymean(~PropST, s)))
(se.pm <- sqrt(pm[[1]] *(1- pm[[1]])/N))

Answer 1

Choosing which population characteristic to estimate

These two statistics you describe are not estimates of the same thing. Rather, they are estimates for two different population characteristics. Which population characteristic you care about will depend on the substantive research question you're trying to answer.

The population ratio $\frac{N_{successes}}{N_{trials}}$ is a essentially a summary of the average trial. This answers the question: "if a trial happens, how likely is it to be successful?"

The population mean $\sum_{i \in U} n_{successes,i} / n_{trials,i}$ is essentially a summary of the average person. This answers the question: "for the average person in the population, how likely is it that their average trial is successful?"

You can imagine these two measures being very different if there is a correlation between individuals' success rate and their number of trials. For example, there might be a few prolific fishers who make lots of fishing trips (i.e., "trials") and have lots of successes on those trips.

Standard errors

In both cases, the estimated standard errors from the survey package are reflecting the fact that you only sampled a portion of the persons in the population. Because you're ultimately interested in trips but your sampling units are persons, you effectively have a cluster sample. That's fine. As long as you make sure to take this clustering into account in your analysis, you should be accounting for most or all of your sampling variance.

In general, it is not a good idea to use the following code to estimate your standard errors. That formula is incoherent, because it's combining estimated proportions among trials with the denominator N of persons; in other words, it's confusing sampling units (people) and population units (trials).

(se.pr <- sqrt(pr[[1]] *(1- pr[[1]])/N))

Whichever population characteristic you decide you're interested in, for the standard error estimates it's best to just use the methods from the survey package, which are tailored to the sampling design.

(pr <- svyratio(~NSuccess,~NTrips,s))
(se.pr <- SE(pr))
(pm <- (svymean(~PropST, s)))
(se.pm <- SE(pr))

Recommended reading

Sharon Lohr's classic "Sampling: Design and Analysis" provides an excellent introduction to working with survey data with features such as stratification or clustering. It also has a nice discussion of population ratios vs. sample means and why you would estimate one or the other.