As usual, differences between a good statistical approach and a scientific mistake (or fraude) are in the way procedures are conducted and the researcher's motivation.

You can use statistical modelling as an exploratory tool, like in the cited Framingham study, when we did not know anything about the factors to influence the risk of cardiovascular diseases. Once you have found a statiscally significant relationship, you can try to understand the real reason behind that relationship and your research can go further, starting more sensible research methods, like clinical trials, for instance.

To build a predicton model, on the other hand, you can start from any model that fits your data. But simply fitting your data is not enough! You will need a validation of your model in other datasets.

The problem is that, as researchers usually understand about their own fields and do not understand well about statistics, and do not want (can?) to pay a good statistician, they perform statistical analyses as they think to know... Something like me, a statistician, trying to run a RNA extraction and a PCR experiment without proper training. I think "oh! I learned about RNA at school! So, I can do it myself!"...

After that, my RNA extraction results will be evaluated by other statisticians, as trained about RNA as myself (or less), and they will decide if I did it well...

Is not that what happens in peer-review process?

christopolous raised a very valid point - though I feel it is not put across properly. Most organised studies do take the assistance of good stasticians and are able make good predictions. we should be judicious and that is all we can appeal to all the researchers. Dr ramadan's perspective is well taken as well.

All the standard adjustments you mentioned are simple the OLS ( Ordinary Least Squares) model requirements. There are also other regression models, a great collection of such models.

The problem is not the applicability or not of the models, the problem is that models are the easy way: We guess a functional form of the solution and try to fit it with our data.

The investigation for the laws in Physics starts from more fundamental aspects and not on the inspiration involved in model building.

Perhaps we can find a success story of regression at almost every discipline, but is there any study of the failures? How many failures have occurred because of a simple model that didn' t take into account a serious factor? In Medicine for example, is there any article counting the regression failures?

As for the social sciences, just remember that almost every paper in macro-economics that concludes with the need for wage-cutting has several model assumptions and regressions inside.

So, Econometrics has become a powerful tool for the 'scientific' justification of all the austerity measures of the recent years.

Look outside your discipline and you will see the regression trap in science...

@Sandro your view is quite balanced, honour and I am with you.

My point for dialogue is my opposite to the 'religion of regression', not the applications of it. If regression analysis had been used just for estimations of the unknown laws we should probably agree with that usage.

The problem is that all models created to fit data are presented as pure science in a way that you can not disagree with them, according to the authors.

So, a general scientific conversation has to start about the limits in regression usages...

These misconceptions always happen as a result of misunderstanding about the methods...

I work with microarray data analysis. It is a powerfull tool to genome research. But microarrays are, in nature, a hypothesis generetor tool. You take a look into the entire genome and then starts to look closely to some aspects. It is not a hypothesis testing tool at all! But a fast search in the literature will show you a lot of works trying to obtain final understanding about genome interactions from microarray data...

I agree with you, Demetris, that there is a misuse of regression analysis and I truly believe that this is due to the misunderstanding about statistics.

One has to keep in mind that statistics is science of chance. If all "standard adjustments" and requirements of any of regression model or taken care of and fit of the model is highly significant, would the parametres obtained would still be reliable. Perhaps no. It is because the regression line does not pass through all data points. It is simply an approximation and very useful approximation.So far taking into account a "serious factor" (s) is concerned, one should know that know science takes all the factors into account, because it is not possible to know all the factors even in Physics, therefore, there is always under-determination and uncertainty. As long as models work, they are valid, so as long as regression model works and give useful results and estimates of parametres to be used in mathematical model, it is valid and a science.

I agree that we should not be blind worshipers of the "Religion of Regression". Every statistical method makes assumptions, and if a researchers uses the method, they have a responsibility to understand these assumptions and test them if possible. If not possible, they should at least make clear how their interpretation would change if the assumptions of their model were violated-- more false positives? More false negatives? For all terms in the model, or just certain ones? Etc. In my experience that's generally what the statistician recommends, but the main authors want to keep it "brief" and "simple" and ignore serious limitations of their experiment as if they are obscure technical details.

But are you going beyond that to say that regression should not be trusted for hypothesis testing at all? Or that formal hypothesis testing should not be done? If so, then how do we distinguish observations that reflect real differences from ones that are likely to have occurred by random chance? Isn't an imperfect tool, whose limitations are at least somewhat understood better than no tool at all? Or do you have a less flawed tool to recommend? If I knew of such an alternative I would start using it tomorrow (or at least writing an R implementation of it tomorrow, so that I could start using it as soon as possible).

@Alex I shall give you a small experiment-example:

Take the attachment xydata.txt file which has 276 (xi,yi) pairs from an unknown function with error and try, by assuming any regression model you wish, to find the inflection point.

Please keep a log file with:

1)the number of models you' ll choose

2)the relevant inflection points you' ll find

3)the total time you' ll expend on this experiment

Then we shall talk again.

Try it!

First I entered in the data and plotted it. Then, I tried taking the estimated first derivative of the data. Then I smacked myself on the forehead and said "He wants the inflection point, not the maximum" so I took a second derivative of the data and regressed that. No, that didn't look like a good model fit either, inflection point estimated at -0.18? I looked at the raw data again, and it reminded me of something I read in Pinheiro and Bates, which turned out to be the four-parameter logistic model. Here is the final result:

##------ Mon Mar 25 05:00:30 2013 ------##

##------ Final final version ------##

> library(nlme)

> foo foonls foo$pred plot(V2~V1,foo,type='l',ylim=c(-1.5,10))

> lines(pred~V1,foo,col='blue')

> abline(h=c(-1.4,9.9))

> abline(v=5.2,lty=2)

> abline(v=5.2+.8,lty=3)

Total models attempted: 3. Relevant inflection points: 1 (the first model could not have found an inflection point, I was just being silly, the second model found a clearly non-relevant inflection point).

Hadn't noticed the time requirement, but I'd estimate an hour as the upper bound on how long this took.

I'm attaching a log of my R session, edited only to delete the screen-spam resulting from my misuse of the scan() function. dydt.R is a script I wrote for taking empirical first derivatives of the data. Not included because it didn't turn out to be useful for the final result. The dataset that matters is the original one I made, named 'foo'. The 'V1' variable is your x and the 'V2' variable is your y. The 'grp' variable is irrelevant, I put it in because my 'dydt' script expects a grouping variable. Neither of those are necessary. All you need to match my results is the data copied from a tab-delimited spreadsheet into your clipboard and the above commands. If read.delim doesn't work with 'clipboard' as an argument on your system, just somehow get the data into a data.frame named foo with columns V1 and V2, and the rest of the above code should work.

So, what would be a non-regression method of attaining the same goal?

Oh, forgot to include the results of the model for people who don't have time to wade through R code:

> foonls

Nonlinear regression model

model: V2 ~ SSfpl(V1, A, B, xmid, scal)

data: foo

A B xmid scal

-1.4075 9.8734 5.1718 0.8117

residual sum-of-squares: 2.236

Number of iterations to convergence: 0

Achieved convergence tolerance: 1.948e-07

> summary(foonls)

Formula: V2 ~ SSfpl(V1, A, B, xmid, scal)

Parameters:

Estimate Std. Error t value Pr(>|t|)

A -1.407487 0.054599 -25.78

By the way, how did you generate the noise for this thing? Was there a sine or cosine function involved, or just random uniform?

Sorry for the multiple messages, the estimated inflection point is 5.1718487 with a 95% confidence interval of (5.1532729, 5.189852)

Darn, another omission. The model uses variables A, B, xmid, and scal which I had assigned. For starting values, I used a visual estimate of: the minimum of V2 for A, the maximum of V2 for B, the midpoint on V1 midway between the minimum and maximum of V2 for xmid. I wildly guessed scal as "1" to have the correct sign for that parameter and let the model figure out its magnitude.

So, that's where A, B, xmid, and scal come from in the nls() expression above. I actually tried a few departures from these "reasonable" starting values and the model still returns the same final estimates. It does, however, seem to require that these variables already be set before SSfpl() is invoked. Attempting to give the raw starting numbers as arguments gives an error.

@Alex, can you give the estimated function in form of a simple type:

f(x) = ...

just to check how 'close' is your answer to the true function used for data generation?

OK, I found it:

f(x)=A+(B-A)/(1+exp((xmid-x)/scal))

Sorry, but my initial function was not this one!

Any other suggestion?

That's irrelevant. What you asked me for was the inflection point. What was the real inflection point, and how would you find it using the information given without relying on regression?

@Alex, the answer:

The true inflection point was p=5.

The function that generated data was a typical non-symmetrical Gompertz sigmoid:

f(x) = 10*exp(-exp(5)*exp(-x))

The error added to data came from Uniform(-0.10,0.10)

Your answer after simplification is:

f(x) = -1.407487+11.280881/(1+exp(6.371846018-1.232024759*x))

was a symmetrical sigmoid with negative values for the interval x

Wow, cool. Thank you for teaching me something. But, if you had two of these datasets, how would you determine your confidence level that they were obtained from populations with different inflection points versus the same population and the difference in inflection points is due to random variation? Or, if you needed to predict the inflection point of a not yet observed dataset, how would you calculate the confidence interval? Finally (and maybe this will be clear to me after I download the package and read the docs) does this estimate only the inflection point or can it recover the functional form and parameters of the process that generated the data? Again, thank you. Rare to have a substantive debate like this.

In the meantime, I ran some simulations to see how the two methods compare as sampling variability increases.

Let's say you already have the dataset from before, named 'dd'...

#We create 100 additional versions of y, each with increasing random uniform variance, from (-0.2,0.2) to (-10,10):

ddsim

Of course, this was a small number of simulations, they were done on a tiny sliver of the parameter space, and the effect of sample size was not explored. Moreover, a uniform error is implausible for any real-life sample I can think of. It would be interesting to see what the overall behavior is over various error sizes and distributions, and various underlying functions. You could get empirical distributions of the inflection point by resampling findiplist() so then you would have a method that might be less biased than nls at low errors and can be used for predictions and hypothesis testing. Do you have a package for finding extrema? In my field, that seems to be a more common need than inflection points.

Your docs have answered one of my earlier questions:

findiplist() only finds the inflection point, like it says. But doesn't that imply a very strong assumption about the underlying function? What if it was given a trig function? Or a high-order polynomial? Or a straight line? As far as I can tell, it returns "NaN" or an estimate, without any way to validate the assumptions it makes.

Not that there is anything wrong with a method that specializes in estimating one thing only. But then, that's the same "all models are wrong but some are useful" logic motivating regression: the model fits the data within the observed range, you extrapolate beyond the observed range at your own risk, and there is always the possibility that a better model will be found (and there are formal tests for determining if that model is better).

This is a recapitulation of the scientific method itself-- we don't claim to have the truth, just a process that overall trends toward closer and closer approximations of it.

@Alex, I'm impressed.

I' d like to thank you for the valuable time you spend on this debate and I agree that it is a high quality one.

I' d also like to apologize for the demo-like output of inflection package, at the next version I' ll add an option for a single output.

1)ESE method always gives an answer and EDE method many times is not applicable.

I have tested them with data of size 101 for 1000 times and I have produced the five point summary for them.

2)I have other methods, unpublished til now, that can hundle the problem of local extremums, minimum or maximum.

3) There is not any kind of 'very strong assumption about the underlying function' in findiplist(). It has been tested for a big variety of functions: polynomials, trig functions, rational functions etc. For the case of 3rd order polynomial, a renormalization formula also exists as a correction to ESE-bias.

One crucial concept is that both methods give a consistent estimator of inflection point, so statistics is not absent.

An other important aspect is that a classification of curves and x-range is done, so we are able to expect what our method theoretically will give as an output when errors do not exist.

As for our small experiment try this:

rm(list=ls())

gomp=function(x,L,b,k)

{

L*exp(-b*exp(-k*x))

}

xy

You can try also next data sets:

xydata150no.txt (without error)

xydata150.txt (with error U(-0.05,0.05) )

Did you find easy the inflection point?

(Answer=0.6883008876)

For data without error, here are all available iterations with 95% confidence intervals:

> xy xx yy a t(a$aese);

ESE iterations

[1,] 0.6951677856

[2,] 0.7369463089

[3,] 0.6784563760

[4,] 0.7035234901

[5,] 0.6868120807

> t(a$aede);

EDE iterations

[1,] 0.8372147654

[2,] 0.7369463091

[3,] 0.7035234901

[4,] 0.6868120807

[5,] 0.6868120807

[6,] 0.6951677854

[7,] 0.6868120807

> a$first;#a$first[,3]

i1 i2 chi_S,D

ESE 14 70 0.6951677856

EDE 13 88 0.8372147654

> a$aesmout;a$aesmout[1]

mean sdev 95%(l) 95%(r)

ESE method 0.7001812083 0.0225758709 0.6721495621 0.7282128544

[1] 0.7001812083

> a$besmout;a$besmout[1]

mean sdev 95%(l) 95%(r)

ESE from EDE iters 0.6918255035 0.004576607948 0.687022648 0.696628359

[1] 0.6918255035

> a$aedmout;a$aedmout[1]

mean sdev 95%(l) 95%(r)

EDE method 0.7190412274 0.05509457343 0.6680872455 0.7699952094

[1] 0.7190412274

> a$bedmout;a$bedmout[1]

mean sdev 95%(l) 95%(r)

EDE from ESE iters 0.7014345639 0.01427008245 0.6787276783 0.7241414495

[1] 0.7014345639

> a$esmout;a$esmout[1]

mean sdev 95%(l) 95%(r)

ESE method, all results 0.6960033559 0.01597569562 0.6852707375 0.7067359742

[1] 0.6960033559

> a$edmout;a$edmout[1]

mean sdev 95%(l) 95%(r)

EDE method, all results 0.7126388043 0.04428595156 0.6828870973 0.7423905114

[1] 0.7126388043

> a$ipall;a$ipall[1]

mean stdev 95%(l) 95%(r)

All methods, all results 0.7047171622 0.0341753138 0.689564677 0.7198696475

[1] 0.7047171622

> t(a$aese);

ESE iterations

[1,] 0.6951677856

[2,] 0.7369463089

[3,] 0.6784563760

[4,] 0.7035234901

[5,] 0.6868120807

> a$bese;

[,1] [,2] [,3] [,4] [,5] [,6]

ESE results from EDE iters 0.6951677854 0.6951677854 0.6868120807 0.6951677854 0.6868120807 NaN

> t(a$bese);

ESE results from EDE iters

[1,] 0.6951677854

[2,] 0.6951677854

[3,] 0.6868120807

[4,] 0.6951677854

[5,] 0.6868120807

[6,] NaN

> t(a$esm);

ESE all iterations

[1,] 0.6951677856

[2,] 0.7369463089

[3,] 0.6784563760

[4,] 0.7035234901

[5,] 0.6868120807

[6,] 0.6951677854

[7,] 0.6951677854

[8,] 0.6868120807

[9,] 0.6951677854

[10,] 0.6868120807

[11,] NaN

> t(a$bede);

EDE results from ESE iters

[1,] 0.7202348995

[2,] 0.7035234901

[3,] 0.6951677854

[4,] 0.6868120807

> t(a$edm);

EDE all iteratio...

The message was too large.

Anyway, very soon the iterative methods will be available by updating R package inflection.

All the above results can be found easily with the command:

....

>library(inflection)

....

....

> a