16
$\begingroup$

A random population sample was surveyed. They were asked if they eat vegetarian diet. If they answered yes, they were also asked to specify how long they’ve been eating vegetarian diet without interruption. I want to use this data to calculate average length of adherence to vegetarianism. In other words, when someone becomes vegetarian, I want to know long on average they stay vegetarian. Let’s assume that:

  • All respondents gave correct and accurate responses
  • World is stable: popularity of vegetarianism is not changing, average length of adherence is not changing either.

My reasoning so far

I found it helpful to analyse a toy model of the world, where at the beginning of every year two people become vegetarians. Every time, one of them stays vegetarian for 1 year and another for 3 years. Obviously, the average length of adherence in this world is (1 + 3) / 2 = 2 years. Here is a graph that illustrates the example. Each rectangle represents a period of vegetarianism:

an illustration

Let’s say we take a survey in the middle of year 4 (red line). We get the following data:

a table

We’d get the same data if we took the survey at any year, starting year 3. If we just average the responses we get:

(2* 0.5 + 1.5 + 2.5)/4 = 1.25

We underestimate because we assume that everyone stopped being vegetarians right after survey, which is obviously incorrect. To obtain an estimate that is closer to the real average times that these participants would remain vegetarian, we can assume that on average, they reported a time about halfway through their period of vegetarianism and multiply reported durations by 2. In a large survey drawing randomly from the population (like the one I’m analysing), I think this is a realistic assumption. At least it’d give a correct expected value. However, if doubling is the only thing we do, we get average of 2.5, which is an overestimate. This is because the longer person stays vegetarian, the more likely (s)he is to be in the sample of current vegetarians.

I then thought that the probability that someone is in the the sample of current vegetarians is proportional to their length of vegetarianism. To account for this bias, I tried to divide number of current vegetarians by their predicted length of adherence:

yet another table

However, this gives an incorrect average as well:

(2*1 + ⅓ * 3 + ⅕ * 5)/(2 + ⅓ + ⅕) = 4 / 2.533333 = 1.579 years

It’d give the correct estimate if number of vegetarians were divided by their correct lengths of adherence:

(1 + ⅓ * (1 + 3 + 5))/(1 + ⅓ * 3) = 2 years

But it doesn't work if I use predicted lengths of adherence and they are all I have in reality. I don’t know what else to try. I read a bit about survival analysis but I’m not sure how to apply it in this case. Ideally, I would also like to be able to calculate a 90% confidence interval. Any tips would be greatly appreciated.

EDIT: It may be possible that the question above has no answer. But there was also another study that asked a random sample of people if they are/were vegetarian and how many times they've been vegetarian in the past. I also know age of everyone in both studies and some other things. Maybe this information can be used in conjunction to the survey of current vegetarians to get the mean somehow. In reality, the study I talked about is just one piece of the puzzle, but a very important one and I want to get more out of it.

$\endgroup$
  • 1
    $\begingroup$ That's not an option atm. This data definitely provides some evidence for length of adherence, I just don't know how to use it. $\endgroup$ – Saulius Šimčikas Aug 19 '17 at 19:57
  • 1
    $\begingroup$ At least one of your images appears to have disappeared (403 error when I use the URL). $\endgroup$ – barrycarter Aug 19 '17 at 21:02
  • 2
    $\begingroup$ @kjetilbhalvorsen For the problem in doesn't matter if vegetarians keep being vegetarians for life. At some point, they will stop being vegetarians, either by eating meat or by dying. $\endgroup$ – Pere Aug 19 '17 at 22:25
  • 4
    $\begingroup$ @kjetil Your "real vegetarians" comment sounds somewhat like a No True Scotsman. The ordinary definition of a vegetarian says nothing about what might happen in the future, nor about why someone is vegetarian, but only about their behavior at the time the attribute is being considered. If someone is vegetarian now, they're vegetarian now, for whatever reason they happen to be one. I don't think our personal feelings about the idea of eating meat or the reasons why we might feel was we do is on topic here; they belong somewhere else. $\endgroup$ – Glen_b -Reinstate Monica Aug 20 '17 at 0:11
  • 2
    $\begingroup$ Since people who are vegetarians for longer are more likely to be selected to appear in your sample, this means that the probability density function of your sample data is proportional to one minus the cumulative distribution function of the adherence lengths. To make an example out of your example, the distribution of lengths is [0, 0.5, 0, 0.5] (50% last for 1 year, 50% for 3 years), giving a CDF of [0, 0.5, 0.5, 1], with one minus that being [1, 0.5, 0.5, 0] which is proportional to the [2, 1, 1, 0] counts of your sample. $\endgroup$ – PhiNotPi Aug 20 '17 at 2:36
10
$\begingroup$

Let $f_X(x)$ denote the pdf of adherence length $X$ of vegetarianism in the population. Our objective is to estimate $EX=\int_0^\infty xf_X(x)dx$.

Assuming that the probability of being included in the survey (the event $S$) is proportional to $X$, the pdf of adherence length $X$ among those included in the survey is $$ f_{X|S}(x) = \frac{xf_X(x)}{\int x f_X(x) dx}=\frac{xf_X(x)}{EX}. $$ At the time of being included in the survey, only a time $Z$ has passed. Conditional on $X$ (and $S$), the reported time being a vegetarian is uniform with pdf $$ f_{Z|X=x}(z) = \frac1x, 0<z<x. $$ Hence, using the law of total probability, the overall distribution of time $Z$ passed as vegetarian among those included in the survey becomes \begin{align} f_Z(z) &= \int_z^\infty f_{Z|X=x}(z)f_{X|S}(x)dx \\&= \int_z^\infty \frac1x \frac{xf_X(x)}{EX}dx \\&= \frac{1-F_X(z)}{EX}, \end{align} where $F_X(z)$ is the cdf of $X$. Since $X$ is a positive variable $F_X(0)=P(X\le 0)=0$ and so $f_Z(0)=1/EX$.

This suggests estimating $EX$ by perhaps first estimating $f_Z(z)$ non-parametrically from the observed data $z_1,z_2,\dots,z_n$. One option is kernel density estimation, using Silverman's reflection method around $z=0$ since the domain of $f_Z(z)$ has a lower bound at $z=0$. This method applied to simulated data is shown as the red curve in the figure below. Having obtained an estimate $\hat f_Z(0)$ of $f_Z(z)$ at $z=0$, an estimate of $EX$ is then given by $\widehat{EX}=1/\hat f_Z(0)$.

enter image description here

This non-parametric method is not ideal however since it does not exploit the fact that $f_Z(z)$ is a non-increasing function. Also, if $f_X(0)=F_X'(0)>0$, $f_Z(0)$ may become severely underestimated and $EX$ overestimated. Finding an estimate of $EX$ in such situations without making more assumptions seems difficult, essentially because short adherence times present in this situation hardly show up in the observed data as a result of the biased sampling.

Alternatively, one could make some distributional assumptions about $f_X(x)$ and fit a parametric model by maximising the likelihood $$ L(\theta)=\prod_{i=1}^n \frac{1-F_X(z_i;\theta)}{EX(\theta)} $$ numerically (blue curve in above figure).

R code simulating data and implementing both methods:

# Simulate lognormal duration length in population
set.seed(1)
n <- 1e+4
x <- rlnorm(n,mean=2,sd=.2)
# Biased sampling
x.given.S <- sample(x, size=n/10, prob=x, replace=TRUE)
# Duration at time of sampling
z <- runif(length(x.given.S),min=0, max=x.given.S)
hist(z,prob=TRUE,main="")

# Compute kernel density estimate with reflection around z=0
to <- max(x) + 3
fhat <- density(z,from = -to, to=to)
m <- length(fhat$y)
fhat$y <- fhat$y[(m/2+1):m] + fhat$y[(m/2):1]
fhat$x <- fhat$x[(m/2+1):m]
lines(fhat,col="red")
# Estimate of EX
1/fhat$y[1]
# True value (mean of above lognormal)
exp(2+.2^2/2)

# Maximum likelihood
nll <- function(theta, z) {
  - sum(plnorm(z, theta[1], theta[2], log.p=TRUE, lower.tail = FALSE)) + length(z)*(theta[1] + theta[2]^2/2)
}
fit <- optim(c(0,1),nll,z=z)
fit$par
EXhat <- exp(fit$par[1]+fit$par[2]^2/2) # MLE of EX
EXhat
curve(plnorm(z, fit$par[1], fit$par[2], lower.tail=FALSE)/EXhat, xname="z", col="blue",add=TRUE)
$\endgroup$
  • 1
    $\begingroup$ Hey, thank you very much for answering, I haven't yet took the time to understand everything, just wanted to add that I do know general distribution from that another study. (the only problem with the other study it that it made people choose between options for how long they've been vegetarian and one of the options was "More than 10 years" and the average depends almost entirely on how much longer than 10 years people remain vegetarian) $\endgroup$ – Saulius Šimčikas Aug 20 '17 at 12:47
  • $\begingroup$ Ok, I hope there are no major flaws in my reasoning. I see that @PhiNotPi arrive at the same pdf in his comment to the OP. $\endgroup$ – Jarle Tufto Aug 20 '17 at 13:06
  • $\begingroup$ @Saulius If your have access to the second right censored data set and the underlying distributions indeed can be assumed to be identical, then the ideal solution would be to combine the likelihood for that data set (which is straightforward to write down if it is just some right censored sampling) and then maximise the joint likelihood. $\endgroup$ – Jarle Tufto Aug 20 '17 at 13:29
  • $\begingroup$ that one is not right censored: imgur.com/U8ofZ3A I now realise that I had to mention this at the start but I thought that my problem had some more straightforward solution... $\endgroup$ – Saulius Šimčikas Aug 20 '17 at 13:39
  • $\begingroup$ @Saulius Those data are interval censored. Again, it is straightforward to compute the likelihood. $\endgroup$ – Jarle Tufto Aug 20 '17 at 13:44
0
$\begingroup$

(I've dithered over adding this, as it appears @JarleTufto has already given a nice mathematical approach; However I'm not clever enough to understand his answer, and now I'm curious if it is exactly the same approach, or if the approach I describe below ever has its uses.)

What I would do is guess an average length, and guess a few distributions around it, and then, for each, make a simulation of my population, and sample it regularly.

You said to assume the total population of vegetarians is not changing, so each time my model has someone stopping, a brand new vegetarian is created. We need to run the model for a number of simulated years to make sure it has settled down, before we can start to sample. After that I think you can take samples every simulated month (*) until you have enough to form your 90% confidence interval.

*: or whatever resolution works with your data. If people gave their answer to the nearest year, sampling every 6 months is good enough.

Out of all your guesses, you choose the mean and distribution which (averaged over all the samples you took) gives you the closest result to what your real-life survey gave.

I would iterate my guesses a few times, to narrow in on the best match.

The best distribution may not be single-peaked. The ex-vegetarians I personally can think of stopped because of major lifestyle changes (typically marrying/living with a non-vegetarian, or moving country, or falling seriously ill and a doctor suggesting it might be diet); on the other side is the power of habit: the longer you've been a vegetarian the more likely you are to carry on being one. If your data had asked age and relationship status, we could throw that in the above simulation too.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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