Simulating data from skewed distribution with zeros

Question

I have some sample data:

   dat <- c(3, 3, 4, 0, 0, 1, 2, 2, 0, 0, 5, 1, 3, 0, 2, 6, 9, 2, 0, 0, 
0, 5, 9, 1, 1, 0, 12, 2, 5, 3, 8, 10, 2, 5, 0, 0, 6, 1, 0, 10, 
5, 0, 2, 0, 1, 1, 1, 9, 2, 4, 0, 5, 2, 0, 0, 0, 8, 1, 1, 7, 0, 
0, 0, 4, 0, 6, 11, 5, 0, 4, 1, 1, 3, 1, 1, 5, 0, 0, 0, 0, 5, 
0, 6, 2, 0, 0, 0, 5, 0, 0, 0, 1, 3, 1, 2, 5, 1, 1, 5, 0, 0, 4, 
1, 10, 2, 0, 2, 5, 0, 0, 0, 2, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 
2, 0, 0, 0, 0, 0, 0, 1, 0, 2, 1, 3, 2, 0, 0, 1, 0, 5, 4, 10, 
1, 0, 1, 1, 1, 3, 3, 0, 1, 0, 0, 0, 1, 8, 0, 2, 1, 2, 2, 5, 1, 
3, 2, 1, 0, 3, 3, 8, 0, 0, 2, 1, 2, 0, 0, 3, 5, 1, 0, 5, 0, 3, 
0, 5, 0, 0, 4, 3, 1, 4, 0, 0, 2, 1, 4, 7, 0, 2, 3, 2, 1, 2, 5, 
2, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 1, 6, 5, 6, 0, 0, 3, 0, 0, 0, 
0, 0, 3, 0, 7, 0, 1, 0, 1, 2, 0, 0, 1, 0, 2, 0, 0, 6, 0, 0)

The histogram looks like this:

The data correspond to some measure of disability and is capped at 12. I'm trying to simulate data for a treated sample (e.g. the above data) and a control sample (e.g. median of the distribution is shifted to the right by 1).

Question Is there a way to simulate the control sample non-parametrically? I want to preserve the zeros, so I cannot simply bootstrap and add 1 to the values. Another option is using some mixture - atom at zero and some distribution $F^+$ on the non-zero part, but it's hard to know what the proportion of zeros and non-zeros will be in the control sample.

Answer 1

Thanks for sharing your dataset dat.

I am assuming you want a fictitious control sample ctrl, of about the same size $n = 246$ as dat, of similar shape, but with a sample median of $2$ instead of $1.$ Also, I'm assuming you'd want to use a two-sample Wilcoxon rank sum test (not a paired Wilcoxon signed rank test) to compare the two samples: actual dat and fictitious ctrl.

Depending on the mechanism for decreasing scores from ctrl to dat you have in mind, there would be many ways to simulate the fictitious sample ctrl. One would be that ctrl would be $0, 1$ or $2$ points higher than dat, but without exceeding $12.$

Importing your dataset dat into R, I get the following summary for dat.

summary(dat)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.000   1.963   3.000  12.000 
n = length(dat);  n
[1] 246
table(dat)
dat
  0   1   2   3   4   5   6   7   8   9  10  11  12 
103  41  31  18   9  21   7   3   4   3   4   1   1 

Then increase the non-zero dat scores by $0,1,$ or $2$ in appropriate proportions to get a preliminary ctrl.1.

ctrl.1 = dat[dat>0]
table(ctrl.1)
ctrl.1
 1  2  3  4  5  6  7  8  9 10 11 12 
41 31 18  9 21  7  3  4  3  4  1  1 

set.seed(2021)
ctrl.1 = ctrl.1 + sample(0:2, 143, rep=T, p=c(1,2,2))
table(ctrl.1)
ctrl.1
 1  2  3  4  5  6  7  8  9 10 11 12 
 7 30 27 20 15 10 17  6  4  2  1  4 

At this point it would be possible to have a few scores above 12. If so, either drop them or maintain the same sample size by changing them to 12.

Now put the $0$'s back, and see what summary we have.

ctrl = c(rep(0,103), ctrl)
summary(ctrl)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   2.000   2.622   4.000  12.000 

The sample median and upper quartile are increased by about $1,$ as desired. Here are boxplots (dat at left) of the two samples. They seem to have similar shapes. Also, the nonoverlapping notches in the sides of the boxes indicate that the two medians may be significantly different at the 5% level.

boxplot(dat, ctrl, notch=T, col="skyblue2", pch=20)

Moreover, a two-sample Wilcoxon RS test can distinguishing between datasets dat and ctrl; the P-value is $0.042 < 0.05 = 5\%.$

wilcox.test(dat, ctrl)

        Wilcoxon rank sum test with continuity correction

data:  dat and ctrl
W = 27182, p-value = 0.04222
alternative hypothesis: 
  true location shift is not equal to 0