11
$\begingroup$

I have a collection of data, which I originally thought was normally distributed. Then I actually looked at it, and realised it wasn't, mostly because the data is skewed, and I also did a shapiro-wilks test.

I'd still like to analyse it using statistical methods, and so I'd like to hypothesis test for skew-normality.

So I'd like to know if there's a way to test for skew normality, and if possible, a library to doing the test for me.

$\endgroup$
8
$\begingroup$

Regarding how to fit data to a skew-normal distribution You could calculate the maximum likelihood estimator from first principles. First note that the probability density function for the skew normal distribution with location parameter $\xi$, scale parameter $\omega$ and shape parameter $\alpha$ is

$$ \frac{2}{\omega} \phi\left(\frac{x-\xi}{\omega}\right) \Phi\left(\alpha \left(\frac{x-\xi}{\omega}\right)\right) $$

where $\phi(\cdot)$ is the standard normal density function and $\Phi(\cdot)$ is the standard normal CDF. Note that this density is a member of the class described in my answer to this question.

The log-likelihood based on a sample of $n$ independent observations from this distribution is:

$$ -n\log(\omega) + \sum_{i=1}^{n} \log \phi\left(\frac{x-\xi}{\omega}\right) + \log \Phi\left(\alpha \left(\frac{x-\xi}{\omega}\right)\right)$$

It's a fact that there is no closed form solution for this MLE. But, it can be solved numerically. For example, in R, you could code up the likelihood function as (note, I've made it less compact/efficient than possible to make it completely transparent how this calculates the likelihood function above):

set.seed(2345)

# generate standard normal data, which is a special case
n = 100 
X = rnorm(n) 

# Calculate (negative) log likelihood for minimization
# P[1] is omega, P[2] is xi and P[3] is alpha
L = function(P)
{

    # positivity constraint on omega
    if( P[1] <= 0 ) return(Inf)

    S = 0
    for(i in 1:n) 
    {
        S = S - log( dnorm( (X[i] - P[2])/P[1] ) ) 
        S = S - log( pnorm( P[3]*(X[i] - P[2])/P[1] ) ) 
    }


    return(S + n*log(P[1]))
}

Now we just numerically minimize this function (i.e. maximize the likelihood). You can do this without having to calculate derivatives by using the Simplex Algorithm, which is the default implementation in the optim() package in R.

Regarding how to test for skewness: We can explicitly test for skew-normal vs. normal (since normal is a submodel) by constraining $\alpha = 0$ and doing a likelihood ratio test.

# log likelihood constraining alpha=0. 
L2 = function(Q) L(c(Q[1],Q[2],0))

# log likelihood from the constrained model
-optim(c(1,1),L2)$value
[1] -202.8816

# log likelihood from the full model
-optim(c(1,1,1),L)$value
[1] -202.0064

# likelihood ratio test statistic
LRT = 2*(202.8816-202.0064)

# p-value under the null distribution (chi square 1)
1-pchisq(LRT,1)
[1] 0.1858265

So we so not reject the null hypothesis that $\alpha=0$ (i.e. no skew).

Here the comparison was simple, since the normal distribution was a submodel. In other, more general cases, you could compare the skew-normal to other reference distributions by comparing, for example, AICs (as done here) if you're using maximum likelihood estimators in all competing fits. For example, you could fit the data by maximum likelihood under a gamma distribution and under the skew normal and see if the added likelihood justifies the added complexity of the skew-normal (3 parameters instead of 2). You could also consider using the one sample Kolmogorov Smirnov test to compare your data with the best fitting estimate from the skew-normal family.

$\endgroup$
  • 1
    $\begingroup$ +1, I think this is a clear, thorough & constructive answer. I have 1 nitpick / concern in the final paragraph re the use of the AIC. An issue I have w/ the various information criteria is that they assume all parameters contribute equally to the ability of the model to fit the data. When evaluating different multiple regression models, I think this is fine; however if examining different types of distributions, it is not clear to me a-priori that all parameters afford equal flexibility. Thus, I am uncomfortable w/ that. What is your position on this issue? $\endgroup$ – gung May 28 '12 at 18:04
  • $\begingroup$ +1 I am just slightly worried about some issues with the Azzalini skew normal such as: (1) the Fisher information matrix of $\alpha$, the skewness parameter, is singular at $\alpha=0$ which indicates inferential problems at this point, specially when using the likelihood ratio statistic; (2) The profile likelihood of $\alpha$ is typically very flat; (3) It has two inflection points and the joint MLE of $(\mu,\sigma,\alpha)$ does not exists for some data sets. $\endgroup$ – user10525 May 28 '12 at 18:31
  • $\begingroup$ @gung, this is a good point. I was using AIC as an example more than anything - something else could be used - but I've seen people use AIC to compare models with different error distributions, which is effectively making this same assumption that all parameters are "created equal". Have you seen any literature on this subject? I'd be interested. $\endgroup$ – Macro May 28 '12 at 18:34
  • $\begingroup$ @Procrastinator, when $\alpha=0$, you have the regular normal distribution as a sub-model. It is straightforward to make this comparison without inverting the fisher information. Re: (2), yes this is true, it's also true that the level of skewness achievable with this distribution is pretty modest (maxes out at around $\pm .9$ which happens as $\alpha$ diverges to $\pm \infty$), so anytime the data set displays that level of skew, the likelihood of $\alpha$ will be pretty flat (i.e. slowly increasing as $\alpha$ diverges). I suppose that last remark was also related to (3) $\endgroup$ – Macro May 28 '12 at 18:39
  • 1
    $\begingroup$ @Macro Thanks for this. That also depends on how you measure skewness, moment-based measures are not very used nowadays because they only exist for ligth tailed distributions. The problem is also that the profile likelihood of $\alpha$ has two inflection points (as shown in my first link). In Azzalini's website, he also mentions that the MLE does not exist for some data sets that he also characterise. This is a delicate point about this distribution which has produced lots of criticisms. Thanks for the discussion. $\endgroup$ – user10525 May 28 '12 at 18:45
5
$\begingroup$

I am a statistician who has been working in this profession for over 30 years and before reading this post I had never heard of the skew normal distribution. If you have highly skewed data why do specifically want to look at skew normal as opposed to lognormal or gamma? Anytime you have a parametric family of distributions such as the gamma, lognormal or skew normal you can apply a goodness of fit test such as chi-square or Kolmogorov-Smirnov.

$\endgroup$
  • 5
    $\begingroup$ The Azzalini skew normal is a popular distribution proposed in 1985. It has support on the whole real line. $\endgroup$ – user10525 May 28 '12 at 11:40
  • 1
    $\begingroup$ @Procrastinator I know that now and probably should have heard of it before. But I guess my point is that since I hadn't heard of it perhaps it is a little more obscure than other skewed distributions. I get the point about the support over the whole real line as opposed to my examples which have support only on [0.∞) or {a,∞) if a shift parameter is added. These distributions are all just approximations to describe how the data are distributed. Do we really know that all negative values are possible? In practical cases the data probably have lower and upper bounds. $\endgroup$ – Michael Chernick May 28 '12 at 12:07
  • 1
    $\begingroup$ @Procrastinator That has nothing to do with my comment. I am saying that real data are often really bounded even when they can be well approximated with unbounded distributions. $\endgroup$ – Michael Chernick May 28 '12 at 13:47
  • 1
    $\begingroup$ @Procrastinator Not quite. If you are given a finite set of number you can't tell from the data whether or not they come from a discrete or continuous distribution. The same is true for boundedness. I am saying that apart from the data you would know strictly on the basis of what you are measuring whether or not it is bounded and continuous or discrete. If for example you are measuring weight of a person you know the weight is greater than 0 and bounded above by physical limitations say 5000 pounds. $\endgroup$ – Michael Chernick May 28 '12 at 14:21
  • 1
    $\begingroup$ Also even though weight measurements can only be determined to a certain number of decimal places it is reasonable to treat weight as continuous. Now if you are going to flip a coin 10 times you know that the number of heads you will get must be an integer between 0 and 10 (so discrete and bounded). My point is that bounds on the distribution are usually very apparent. It is not as clear when deciding between continuous and discrete distributions. $\endgroup$ – Michael Chernick May 28 '12 at 14:26
1
$\begingroup$

So my solution in the end was to download the fGarch package, and snormFit provided by fGarch to get MLEs for the parameters to a Skewed-Normal.

Then I plugged those parameters, with the dsnorm function provided by fGarch, in to a Kolmogorov-Smirnov test.

$\endgroup$
  • $\begingroup$ How can you combine MLE and the Kolmogorov-Smirnov test, which is a nonparametric test? $\endgroup$ – user10525 May 28 '12 at 13:17
  • $\begingroup$ I'd like to point out that I have no idea what I'm doing, and just wandering along blindly. I assumed that KS worked like chi-squared, and it looked at the differences between what sample data I have, and what the distribution itself should be. R's ks.test accepts firstly the sample data, and then a distribution, along with the parameters to that distribution. I used the MLEs as the parameters. Also, my assumption/remberance of how Chi-squared works could also have been wrong... $\endgroup$ – Squidly May 28 '12 at 13:27
  • $\begingroup$ I think you have to be careful with this approach. First, you have to answer: do you want to conduct a goodness of fit test for skew normality of your data? or, do you want to check if the data is normal or skew normal, i.e. $H_0: \lambda=0$? $\endgroup$ – user10525 May 28 '12 at 13:29
  • $\begingroup$ The second one $ H_0 = \lambda = 0 $ I was assuming that if the fit was poor enough, I'd get a tiny P value, and I could happily reject $H_0$. Which I did. $\endgroup$ – Squidly May 28 '12 at 13:32
  • 1
    $\begingroup$ @Procrastinator There are many goodness of fit tests based on the empirical cdf. Kolmogorov Smirnov is one. These test can be used to compare the empirical cdf to any specific distribution (and with adjustments when unknown parameters are estimated before testing. You are absolutely right that rejecting a distribution for normality for example does not tell how the distributions differ. But MrBones if he wants to do this formally he can test for significant skewness and then do a KS or chi square test for skew normal. Adjustments can be made for multiple testing. Why be smug with the OP? $\endgroup$ – Michael Chernick May 28 '12 at 14:55
0
$\begingroup$

Check out http://www.egyankosh.ac.in/bitstream/123456789/25807/1/Unit6.pdf and http://en.wikipedia.org/wiki/Skewness

You could use the Karl Pearson test for skewness. The ratio of the third moment to the cube of standard deviation is called the coefficient of skewness. Symmetrical distributions would have skewness = 0

$\endgroup$
  • 2
    $\begingroup$ I'm not wanting to work out the skew of my data. I know it's skewed. I want to know if my data follows a skew-normal distribution. $\endgroup$ – Squidly May 28 '12 at 11:02
  • 2
    $\begingroup$ Which software are you using, R has a package 'sn' (skew normal) that contains functions that calculate the ML estimates. I am not sure of the exact functions though - check out this site azzalini.stat.unipd.it/SN for details on the R package $\endgroup$ – NaN May 28 '12 at 11:21
0
$\begingroup$

in SPSS you can get a estimate of the skewness (by going to analyze and then descriptives and then mark skewness) then you get a score of skewness and S.E (standard error) of skewness. Divide the skewness by its S.E and if your score is between +-1.96 its normally skewd. If its not skewd then there are many non-parametric tests out there! Good luck and all the best!

$\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.