5
$\begingroup$

I am trying to model data about altruistic behavior in a simple lab experiment. I have one value for each participant in the sample (N=479), describing how altruistic that person was. As you can see in the graph below, the data are positively skewed, with most people having an altruism score of 0. However, there are also a few people with a negative altruism score - these "spiteful" people are willing to pay a cost to hurt others.

enter image description here

For simulation purposes, I am trying to fit this data to a parametric distribution. However, most distributions that can be used to model positively-skewed data, like gamma or lognormal, cannot take negative values.

One solution to this problem would be to transform my data by adding a constant so that all negative values become positive, and then fit the data to, e.g., a gamma distribution. However, this does not feel entirely satisfactory: while a gamma distribution has a strict lower bound (at 0), the variable I am trying to model doesn't (there is no hard theoretical limit on how spiteful someone could be).

Is there a probability distribution that allows for skew and also can take negative values?

edit: here is the data

  0.00  0.00  0.15  0.74  0.15  0.00  0.35  0.00  0.55  0.00  0.15  0.15 1.55  0.15  0.55  0.00  0.15  0.00 0.15 -0.45  0.15  0.75  0.00  0.35  0.00  0.00  0.00  0.00  0.15  0.55  0.35  0.00  0.35  0.00  0.55  0.00 0.15  1.55  0.35  1.55  0.00  0.15  0.00  0.15  1.35  0.15  0.15  0.00  0.00 -0.26  0.00  0.00  0.55  0.15 0.55  0.15  0.00 -0.45  1.35  0.55  0.00  0.35  0.15  0.35  0.95  0.15  0.00  0.00  0.15  0.00  0.15  0.00 0.00  0.00  0.75  0.00  0.00  0.00  0.00  0.15  0.55  0.74 -0.26  1.55  0.15 -0.26  0.00  0.00  0.15  0.00  0.00  0.15  0.15  0.15  0.15  0.00  0.15  0.00 -0.26  0.35  0.00  0.35  1.55  0.00  1.55  0.15  0.00  0.00  0.00  0.00  0.55  0.00  0.15 -0.45  0.94 -0.26  0.15  0.15  0.00  0.00  0.00 -0.45  0.00  0.00  0.00  0.35  1.55 -0.26  0.55  0.55  0.35  0.35  0.35  0.00  0.00  0.55  0.74  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.15  0.15  0.55  0.15  0.15  0.00  0.00  0.00  0.75  0.15  0.55  0.00  0.00  0.15  0.15  0.55  0.35  0.00  0.00  0.00  1.35  0.00  0.00  0.00  0.75  0.15  0.00 -0.26  0.00 -0.45  0.15  0.15  0.00  0.00  0.35  0.00  0.55  0.35  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.94  1.55  0.15  0.00  0.55  0.00  0.00  0.00  1.55  0.15  0.15 -0.26  0.35  0.00  0.15  0.15  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  0.15  0.00  0.00  0.94  0.00  0.35  0.00  0.00  0.55  0.55  0.15  0.35  0.00  0.00 -0.26  0.00  0.15  0.00  0.15  0.00  0.35  0.00  0.15  0.00  0.15  1.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15 0.00 -0.45  0.00  0.00  0.74  0.74  0.00  0.35  0.15  0.55  0.00  0.15  0.15  0.15  0.55  1.55  0.00  0.00  0.00 -0.45  0.15  1.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  0.55  0.94  0.00  0.15  0.00  0.00  0.35  0.94  0.75  0.00  0.00  0.15  0.00  0.15  0.55  0.35  0.00  0.00  0.35  1.55 -0.26  0.00  0.00  0.21  0.00  0.00  0.15  0.00  1.15 -0.45  0.00  0.55  0.00  0.35  0.35  0.35  0.15  0.15  0.15  0.00  0.95 -0.26  0.94  0.00  1.35  0.15  0.35  0.00  0.55  0.65  0.00  1.15 -0.45  0.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  -0.26  0.15  0.35  0.35  0.00  0.94  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.94 -0.45  0.00  0.00  0.00 0.35  0.15  0.00  0.00  0.15  1.35  0.00  0.15  0.00  0.15  0.35  0.94  0.35  0.00  0.00  0.55  0.55  0.00  0.00  0.00  0.55  1.35  0.00  0.55  0.55  0.00  0.15  0.35  0.15  0.00  0.00  0.00  1.35  0.15  0.00  0.35  0.35  0.00  0.15  1.55  0.55  0.00  0.35  0.15  0.15  0.15  0.15  0.15  0.94  0.35  0.95  0.15  0.00  0.00  0.00  0.15  0.00  0.55  0.15  0.00  0.15  1.55  0.00  0.75 -0.45  0.00  0.00  0.35  0.00  0.94  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  1.15  0.00  0.15  0.15 -0.26  0.35  0.35  0.35  0.15  0.15 -0.26  0.35 -0.45 0.00  0.55  0.00 -0.26  0.15  0.55 -0.26  0.00  0.15  0.15  0.94  0.94  0.15  0.00  0.00  0.00  0.35  0.55  0.15  0.00  0.00  0.00  0.15  0.00  0.00  0.00  0.00  0.55  0.00
Improve this question
$\endgroup$
6
  • $\begingroup$ Of course there are. Any positive valued continuous function defined on an interval where the lower endpoint is negative is not symmetric and integrates to 1 will work. One such example is the skew normal distribution which you can look up on Wikipedia. $\endgroup$
    – Michael R. Chernick
    Apr 5, 2018 at 22:30
  • $\begingroup$ Do you mind sharing the data? $\endgroup$
    – Georg M. Goerg
    Apr 5, 2018 at 22:58
  • $\begingroup$ With a sample size of 479, why do you feel the need to fit a parametric probability distribution? Why not just get a nonparametric density estimate (such as what one gets with the density function in R)? This is not to suggest that you don't have a need to fit a parametric probability distribution but only that you might want to express why you think that is necessary (assuming that it is necessary). $\endgroup$
    – JimB
    Apr 5, 2018 at 23:23
  • $\begingroup$ I have two reasons in mind: 1)Using a parametric distribution makes it easier to write code that generates data from this distribution (at least, using what I know in R). 2)Assuming I publish the outcome of this research, it makes it easier for readers of the paper to reproduce my results, just by using the parameters of the distribution given in the paper. $\endgroup$
    – Choirofangels
    Apr 6, 2018 at 0:05
  • $\begingroup$ If the parametric distribution is a good fit, then for whatever it's worth I think your reasons are good. If you can't find one with a good fit, then consider stats.stackexchange.com/questions/82797/… to eliminate reason (1). $\endgroup$
    – JimB
    Apr 6, 2018 at 4:25

3 Answers 3

Reset to default
2
$\begingroup$

Not a complete answer, but people will see why this is not a comment. I get 473 values here and suggest that the ambition to fit a smooth distribution here needs to surmount marked granularity in the data. What protocol lies behind this? The small moral that binning can conceal fine structure as well as noise should need little emphasis. The square root scale (compare J.W. Tukey's "rootograms") is chosen as a matter of convenience. 

      score |      freq.           %      cum. %
------------+-----------------------------------
       -.45 |         14        2.96        2.96
       -.26 |         16        3.38        6.34
          0 |        186       39.32       45.67
        .15 |        119       25.16       70.82
        .21 |          1        0.21       71.04
        .35 |         46        9.73       80.76
        .55 |         36        7.61       88.37
        .65 |          1        0.21       88.58
        .74 |          7        1.48       90.06
        .75 |          6        1.27       91.33
        .94 |         13        2.75       94.08
        .95 |          3        0.63       94.71
       1.15 |          5        1.06       95.77
       1.35 |          7        1.48       97.25
       1.55 |         13        2.75      100.00
------------+-----------------------------------
      total |        473      100.00

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ The binning comes from the nature of the experiment (getting less coarse data would have required too many experimental trials). I understand that this level of granularity is less than ideal, but how problematic would you say it is (minor inconvenience or fatal problem)? $\endgroup$
    – Choirofangels
    Apr 6, 2018 at 0:58
  • $\begingroup$ I am more puzzled by the small gaps (e.g. 0.74, 0.75; 0.94, 0.95) than by the large ones! Also, any story on why no-one at around $-$0.15 given that 0 is a strong mode? $\endgroup$
    – Nick Cox
    Apr 6, 2018 at 1:07
  • $\begingroup$ Were there limits on possible scores? My guess is a rescaled beta distribution is a fair starting point. $\endgroup$
    – Nick Cox
    Apr 6, 2018 at 1:40
  • $\begingroup$ There were limits on possible scores: a participant couldn't have a score lower than -.45 or higher than 1.55. However, these bounds are fairly arbitrary: one could have designed a similar experiment with more trials that could have enabled higher and lower scores. In other words, -.45 does not represent "maximum possible spitefulness" but rather the maximum amount of spitefulness that could be detected in this particular experiment. Therefore, ideally I would be looking for a distribution that has no bounds either. $\endgroup$
    – Choirofangels
    Apr 6, 2018 at 4:01
  • 1
    $\begingroup$ As the limits can be attained, I think beta is off the table. I don't have further suggestions as I can't tell between respecting your data, respecting your protocol and respecting your desire to infer the distribution of a hypothetical latent variable. How can anyone guess whether the asymmetry is a fact of psychology or imparted by your procedure? $\endgroup$
    – Nick Cox
    Apr 6, 2018 at 8:22
2
$\begingroup$

Just some ideas, but yours seems to be a difficult case for distribution fitting. I will disregard the bounds for now, and show some example code in R.

text  <-  "0.00  0.00  0.15  0.74  0.15  0.00  0.35  0.00  0.55  0.00  0.15  0.15 1.55  0.15  0.55  0.00  0.15  0.00 0.15 -0.45  0.15  0.75  0.00  0.35  0.00  0.00  0.00  0.00  0.15  0.55  0.35  0.00  0.35  0.00  0.55  0.00 0.15  1.55  0.35  1.55  0.00  0.15  0.00  0.15  1.35  0.15  0.15  0.00  0.00 -0.26  0.00  0.00  0.55  0.15 0.55  0.15  0.00 -0.45  1.35  0.55  0.00  0.35  0.15  0.35  0.95  0.15  0.00  0.00  0.15  0.00  0.15  0.00 0.00  0.00  0.75  0.00  0.00  0.00  0.00  0.15  0.55  0.74 -0.26  1.55  0.15 -0.26  0.00  0.00  0.15  0.00  0.00  0.15  0.15  0.15  0.15  0.00  0.15  0.00 -0.26  0.35  0.00  0.35  1.55  0.00  1.55  0.15  0.00  0.00  0.00  0.00  0.55  0.00  0.15 -0.45  0.94 -0.26  0.15  0.15  0.00  0.00  0.00 -0.45  0.00  0.00  0.00  0.35  1.55 -0.26  0.55  0.55  0.35  0.35  0.35  0.00  0.00  0.55  0.74  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.15  0.15  0.55  0.15  0.15  0.00  0.00  0.00  0.75  0.15  0.55  0.00  0.00  0.15  0.15  0.55  0.35  0.00  0.00  0.00  1.35  0.00  0.00  0.00  0.75  0.15  0.00 -0.26  0.00 -0.45  0.15  0.15  0.00  0.00  0.35  0.00  0.55  0.35  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.94  1.55  0.15  0.00  0.55  0.00  0.00  0.00  1.55  0.15  0.15 -0.26  0.35  0.00  0.15  0.15  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  0.15  0.00  0.00  0.94  0.00  0.35  0.00  0.00  0.55  0.55  0.15  0.35  0.00  0.00 -0.26  0.00  0.15  0.00  0.15  0.00  0.35  0.00  0.15  0.00  0.15  1.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15 0.00 -0.45  0.00  0.00  0.74  0.74  0.00  0.35  0.15  0.55  0.00  0.15  0.15  0.15  0.55  1.55  0.00  0.00  0.00 -0.45  0.15  1.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  0.55  0.94  0.00  0.15  0.00  0.00  0.35  0.94  0.75  0.00  0.00  0.15  0.00  0.15  0.55  0.35  0.00  0.00  0.35  1.55 -0.26  0.00  0.00  0.21  0.00  0.00  0.15  0.00  1.15 -0.45  0.00  0.55  0.00  0.35  0.35  0.35  0.15  0.15  0.15  0.00  0.95 -0.26  0.94  0.00  1.35  0.15  0.35  0.00  0.55  0.65  0.00  1.15 -0.45  0.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  -0.26  0.15  0.35  0.35  0.00  0.94  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.94 -0.45  0.00  0.00  0.00 0.35  0.15  0.00  0.00  0.15  1.35  0.00  0.15  0.00  0.15  0.35  0.94  0.35  0.00  0.00  0.55  0.55  0.00  0.00  0.00  0.55  1.35  0.00  0.55  0.55  0.00  0.15  0.35  0.15  0.00  0.00  0.00  1.35  0.15  0.00  0.35  0.35  0.00  0.15  1.55  0.55  0.00  0.35  0.15  0.15  0.15  0.15  0.15  0.94  0.35  0.95  0.15  0.00  0.00  0.00  0.15  0.00  0.55  0.15  0.00  0.15  1.55  0.00  0.75 -0.45  0.00  0.00  0.35  0.00  0.94  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  1.15  0.00  0.15  0.15 -0.26  0.35  0.35  0.35  0.15  0.15 -0.26  0.35 -0.45 0.00  0.55  0.00 -0.26  0.15  0.55 -0.26  0.00  0.15  0.15  0.94  0.94  0.15  0.00  0.00  0.00  0.35  0.55  0.15  0.00  0.00  0.00  0.15  0.00  0.00  0.00  0.00  0.55  0.00"
dat  <-  scan(textConnection(text))

One idea is to try some Box-Cox transform, but the standard formulation cannot be used as there are some negative values.

library(MASS) # We will use this later
library(car) # We use boxCox from car
boxCox(dat ~ 1, family="yjPower",param="gamma")  # this produces a plot, indicating lambda around -0.5  (not shown)
 boxCox(dat ~ 1, family="bcnPower")  # This gives somewhat conflicting advice, lambda close to zero but positive, again plot not shown

I will leave that kind of method only as an idea here. Now, for a different idea, we can try to fit a skew-normal distribution (for this data I believe more in this last idea than the former one):

library(sn)
skew_mod <-  selm(dat ~ 1) # selm is "skew-elliptic lm"
summary(skew_mod)
Call: selm(formula = dat ~ 1)
Number of observations: 473 
Family: SN 
Estimation method: MLE
Log-likelihood: -159.4612 
Parameter type: CP 

CP residuals:
     Min       1Q   Median       3Q      Max 
-0.69865 -0.24865 -0.09865  0.10135  1.30135 

Regression coefficients
     estimate  std.err  z-ratio Pr{>|z|}
mean  0.24865  0.01651 15.06170        0

Parameters of the SEC random component
       estimate std.err
s.d.     0.3603   0.013
gamma1   0.7498   0.035

and finally we show a histogram with the estimated skew-normal density overlaid:

hist(dat,prob=TRUE,nclass="scott") # "scott" is from MASS
plot(function(x) dsn(x, dp=skew_mod@param$dp), from=-0.5, to=1.75, col="red", add=TRUE)

The probability mass outside the bounds is so low that that should not cause problems. For sampling one could simply reject simulated values outside bounds. But the fit is not too good: the skewness is captured, but the histogram has much more mass at the center. Maybe a skew-t distribution could give a better fit. I will leave that for the OP, but the fitting can be done with the software I have illustrated here. The other answer by Nick Cox address the obvious granularity in the data, I have ignored that aspect, but before doing anything you should ask yourself why is that. If it is important to maintain that granularity in simulations, what I have done here is not enough.

histogram with skew-normal density overlaid

Improve this answer
$\endgroup$
1
$\begingroup$

Thanks to everyone who helped! I have finally found that a skewed laplace distribution provides an excellent fit.

#package rmutil has the skewed laplace distribution
library(rmutil)
#package fitdistrplus is used to fit the data to a distribution
library(fitdistrplus)
#estimate which parameters of a skewed laplace distribution fit the data best
fit.sklaplace <- fitdist(AltruismScore, "skewlaplace", start=list(0,.2,.5))
#this yields the following parameters: m=0.00, s=.23, f=.63.

#plot the distribution thus obtained over the data
range <- seq(-1,2.5,by=.01)
histogram<-hist(AltruismScore, prob=TRUE, xlab="AltruismScore", main="Histogram of AltruismScore", col="lightgreen", xaxt="n")
histogram
lines(range, dskewlaplace(range, 00, .23, .63), col="darkblue", lwd=2)
axis(side=1,at=histogram$mids)

enter image description here

edit: histograms with better bins enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ Your histogram doesn't do you any favours here. Your data show a mode at 0 and that's not quite obvious in the histogram. Evidently your bins are of the form $a < x \le b$ so that when $b = 0$ the zeros are included in the bin $(-0.2, 0]$. I'd try bin limits that go $-0.7(0.2)1.7$ or $-0.65(0.1)1.65$ so that the bin that includes $0$ has $0$ as its midpoint. On the larger issue, I remain sceptical that an unbounded distribution is of that much use or interest for a bounded response. $\endgroup$
    – Nick Cox
    Apr 23, 2018 at 1:27
  • $\begingroup$ I have uploaded a new histogram above following your suggestions. The weird thing is that the axis tick at -.4 is left out for some reason. Code: histogram<-hist(AltruismScore, prob=TRUE, xlab="AltruismScore", main="Histogram of AltruismScore", col="lightgreen", xaxt="n", breaks=seq(-.7,1.7,.2)) histogram lines(range, dskewlaplace(range, 00, .23, .63), col="darkblue", lwd=2) axis(side=1,at=round(histogram$mids, 2)) $\endgroup$
    – Choirofangels
    Apr 23, 2018 at 3:01
  • $\begingroup$ I see a tick at $-0.6$ but no text -0.6. I think R does that when the axis gets too crowded. $\endgroup$
    – Nick Cox
    Apr 23, 2018 at 7:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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