I am running a linear mixed model with 4 fixed factors and 1 random factor. The response variable is %growth and it has negative values (some of my animals shrunk). The problem I'm having is the residuals are not normally distributed according to the shapiro wilk test. The histogram seems like a normal distribution, but the normal quantile plot has an s-curve. In looking this up, I've found this curve labeled fat tail, short tail, and long tail.

I have tried transforms (log, sqrt, cube root) but nothing seems to make it normal (log doesn't work at all bc of the negative values). My next step would be to try a generalize linear mixed model, but I use JMP and it doesn't offer that functionality yet. Do I have any options aside from learning R?

My n is 304, and I have heard that normality matters less as sample size increases, but I don't think 304 is very large.

Edit: I have heard that % growth is not a great response variable for ANOVA, and that perhaps I should do an ANCOVA with final length as the response variable, and initial length as a covariate. However, this also gives me an s-curve and non-normal residuals.

residual histogram and normal quantile plot

  • $\begingroup$ 304 is not huge, but it is large enough for small violations from normality to hopefully be unproblematic. Sadly, "hopefully" is not such a scientifically useful term. I'd be surprised if the problem disappears with a mixed model. The model of final length as a response variable seems risky to me, as you will likely need to interact all your other variables with initial length. Can JMP compute bootstrapped p-values? That would seem to me to be the nice solution here. $\endgroup$
    – Tim
    Feb 28, 2017 at 1:04
  • $\begingroup$ Have you tried using log(final length) as response variable and log(initial length) as a covariate or offset? $\endgroup$
    – EdM
    Feb 28, 2017 at 1:24
  • $\begingroup$ @Tim I will look into bootstrapping today and get back to you. You are right about the interactions with initial length - I tried it and it caused a loss of degrees of freedom for all my fixed effects. $\endgroup$ Feb 28, 2017 at 13:47
  • $\begingroup$ @EdM All of the interaction terms would likely cause a loss of DF, but I will try the log x log ANCOVA $\endgroup$ Feb 28, 2017 at 13:47
  • $\begingroup$ @Tim I'm not sure how I would do bootstrapping in JMP (I've never done bootstrapping at all). My search resulted in a bootstrap function that is only available in JMP Pro which I do not have, and it is performed on one-way analysis. I think even if I had the enhanced program, that option would not work as I have 4 factors plus the random factor. $\endgroup$ Feb 28, 2017 at 16:21

3 Answers 3


Growth rates must be distributed as some variation of the Cauchy distribution. I have written a series of papers on this. The Cauchy distribution has no mean so it has no variance or covariance. You can find my author page at https://papers.ssrn.com/sol3/cf_dev/AbsByAuth.cfm?per_id=1541471

Start with the paper titled "The distribution of returns," and then switch to the paper on Bayesian methods. Generally speaking, there is no admissible non-Bayesian solution though in specific cases there is a maximum likelihood solution that can be used if a null hypothesis method is required. The Bayesian likelihood function is always minimally sufficient.

You can communicate with me through the address on the author page. Because there is no variance or covariance, ANOVA and ANCOVA are impossible.

EDIT With regard to the comments:

1)If I say it's Cauchy, and therefore rule out ANOVA and ANCOVA, what are my options?

Bayesian regression is still available. Your likelihood function would be $$\frac{1}{\pi}\frac{\sigma}{\sigma^2+(y-\beta_0-\beta_1x_1-b_2x_2\dots\beta_nx_n)^2}$$

The meaning is very different from OLS, though. OLS is from a convergent process, such as water going down a drain. This can be conceived as a double pendulum problem, and as such has limited predictive capapcity. For example, if you had $y|x$ you could think of $x$ as the upper pendulum and $y$ as the lower pendulum attached to the top pendulum. So, while $y$ is affected by the movement of $x$ it does not mean they are even moving in the same direction with a positive correlation. A pendulum swinging to the left could cause the other pendulum from momentum to swing to the right.

There is a tight linkage between the double pendulum problem, which is the first real observed problem in chaos theory, and regression in this case.

The proper interpretation is that, for example, if $y=1.1x$ then 50% of the time $y$ will be greater than $1.1x$ and fifty percent of the time it will be less than $1.1x$. You may be able to make stronger statements if there are other properties in your system, such as non-negativity.

2) I've read that Cauchy is problematic if residuals are almost normal.

This does not matter. You can find, or even construct, cases where the Cauchy distribution is indistinguishable from a normal distribution. Generally speaking there is no admissible non-Bayesian solution for most standard problems. This is a problem for someone trained only in Frequentist methods, but is not a problem per se. If a null hypothesis is required by the nature of the problem, then the only close solution would be quantile regression or Theil's regression. The problem with either is that, in the above equation, $x_1$ and $x_2$ are not independent, but they are also not correlated.

The question is not Gaussian versus Cauchy by some empirical test, but which should I have from theory. A certain percentage of the time, data drawn from a pure normal distribution will falsify a test of normality through chance alone. While it is sometimes true we do not know the likelihood function and must test for it, sometimes we do. This is a case where we do.

3) Is this really growth rate if I only have initial and final length?

Yes that is a growth rate, it is just a growth rate with limited observations per creature.

Not asked

Is there anything similar to anova or ancova? The answer is "it is unclear." If you will notice there is only one scale parameter in the likelihood and this does not depend on the number of variables. The scale parameter is a composite of the separate scale parameters, but it is not clear that there is any way to take advantage of this.

  • $\begingroup$ I've sent you a message, but for the benefit of others: 1)If I say it's Cauchy, and therefore rule out ANOVA and ANCOVA, what are my options? 2) I've read that Cauchy is problematic if residuals are almost normal. 3) Is this really growth rate if I only have initial and final length? $\endgroup$ Feb 28, 2017 at 16:26
  • 1
    $\begingroup$ You can't claim that all growth rates are from Cauchy. OP didn't even explain what growth rates are in question here. For instance, stable distributions have been tried on financial and economic data, there were a lot of academic papers on the subj. Yet, almost nobody applies stable distributions to these series anymore in the industry. They don't seem to work really $\endgroup$
    – Aksakal
    Feb 6, 2018 at 14:34
  • $\begingroup$ @Aksakal actually I can claim that. I understand your concern. In raw form, all growth rates or rates of returns must be a deformation of the Cauchy distribution, except for single period discount bonds, some specific cases for other bonds and cash-for-stock mergers. There is a missing ingredient in the older work. I wrote a book with a mathematician covering the expanse of methods of estimation and inference for distributions with heavy tails. It ranged from inference with the sign test to estimation using characteristic functions. $\endgroup$ Feb 7, 2018 at 1:16
  • $\begingroup$ @Aksakal what neither of us realized was how fragile these tools were to their assumptions, or where the inferences may not mean what is intended. We pulled the book from electronic availability subsequently. What triggered this was that I build a model to replace Black-Scholes that neither depends upon the underlying distribution nor the existence of a first moment as well as a parametric version and then tested it empirically. My goal was to test it using Pearson and Neyman's Frequentist methods, Fisher's Likelihoodist methods and Bayesian methods. $\endgroup$ Feb 7, 2018 at 1:20
  • $\begingroup$ @Aksakal The MLE doesn't have a closed form solution and I honestly have no clue how to perform a computable numerical estimate, but I decided I would estimate it by treating it as a MAP estimator with a flat prior. What surprised me was that the Frequentist solution did not match the Bayesian with roughly fifty million observations. $\endgroup$ Feb 7, 2018 at 1:22

Your Q-Q plot doesn't look like it has a fat tail. I'll show you how a fat tail looks like: enter image description here

Your tail is like a Victoria Secret's model compared to the above. I wish some of my model residuals had tails like yours has.


It seems that bootstrapping, as suggested by @Tim in comments, might be a good way to proceed. Even if your statistical software doesn't directly support bootstrapping it's not to hard to roll your own. For example, say you have all data for each individual in a single row (species; 3 treatment types; starting and ending length, width, and height; block ID) and you have 304 rows. You set up an index vector of length 304, and for each bootstrap fill it with a random sample with replacement from the integers from 1 to 304. You then take those indexed rows from your full data set for 304 rows total (with some original rows omitted, some taken once, and others 2 or more times). Then do your analysis and store the regression coefficients. Do this 999 times. For each regression coefficient, average the 999 results; then put its values for the 999 repetitions in rank order; the 25th and 975th in order set the 95% confidence limits. Unless your analysis depends heavily on having a balanced design, this should suffice.

This will not work if there is an underlying Cauchy problem, but I'm not convinced that you have this problem with your data set. The Cauchy problem comes from trying to take a ratio of 2 random variables in which you run a risk of dividing by zero or by numbers close to zero. Although this can be a serious issue in the types of economic time series addressed by @DaveHarris in the documents linked from his answer, in your case the lengths, widths and heights are all positive and far from zero so you don't seem to be in that situation. Raw differences between start and finish, or logs of start/finish ratios, should be sufficiently well behaved that you can analyze your data with the bootstrap, which is a well respected way to deal with your type of situation when you can't count on normal distributions of the data of interest.

  • $\begingroup$ This is a little over my head, but it seems promising...I'll see what I can do. $\endgroup$ Mar 1, 2017 at 1:25
  • $\begingroup$ If you are at some type of academic institution, look into what local statistical support might be available. Also ask around for people who use R; this is very straightforward in R. If you intend to do a lot of projects requiring statistical analysis, consider devoting the time to getting up the initially steep learning curve in R; I did so a dozen years ago and am very glad I did. $\endgroup$
    – EdM
    Mar 1, 2017 at 1:48
  • $\begingroup$ I'm meeting with my R-literate lab mates today - they're suggesting glmm. If that doesn't work, I have the option to pay a small fee for someone from the stats department to do it. I still want to try this bootstrapping method - I found an add-in for jmp that does sampling with replacement, so I'll try that. I know that I have to learn R eventually, but I'm a PhD student and I'm trying to finish before I run out of money. $\endgroup$ Mar 1, 2017 at 13:52
  • $\begingroup$ If you have R-literate labmates they should also be able to help you with the boot package in R, which essentially does what I outlined in the first paragraph of my answer. You just write a wrapper function that takes your complete data set and an index vector, and then returns the values of interest (regression coefficients here) for each bootstap sample. Then boot() does the resampling and boot.ci() returns several types of (empirical) confidence intervals. $\endgroup$
    – EdM
    Mar 1, 2017 at 15:06
  • $\begingroup$ Also, if one of your 2 species seems to have qualitatively different behavior from the other in terms of these residuals, it might make sense to analyze them separately. $\endgroup$
    – EdM
    Mar 1, 2017 at 15:08

Your Answer

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

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