@Joachim: I don't get what you mean with "polynomial". As Timothy said: polynomial regression is different.

In a polynomial regression you use a polynomial of a metric predictor. The mean is the (very) special case of a polynomial regression of order 0: y = b0x0 = b0.The linear regression is the special case of a polynomial regression of order 1: y = b0 + b1x1 = b0 + b1x.The typical polynomial regression models contain higher powers of x, like y = b0 + b1x + b2x2 + b3x3 + ...

If the aim of the regression is to model, for instance, the quadratic relationship between speed and stopping distance of a car, or between the distance of a light source and the light intensity and such, then the linear term (b1x) is not required - it does not even make much sense to include this term in the model, so that an adequate model would look like y = b0 + b2x2. Under certain conditions it may be neccesary or known that for x=0 y must be 0, too, so that b0 should also be removed from the model.

As another example consider Poiseuille's law of a laminar liquid flow through a pipe of radius r. The flow volume is proportional to r4.  If you measure the flow y depending on the radius x (possibly to ftest the hypothesis that the flow is actually laminar), a model of the form y = b4x4 would be appropriate.

In general linear models the word "interaction" is usually understood as a "difference in diffrence" that allows to model synergistic or antagonistic effects (where the change in the response is more [or less, resp.] than the sum of the effects of the predictors considered separately.

Consider two dichotomous factors A and B, like A: person smokes at least 1 pk/d [yes/no]; B: person drinks at least 20 g/d alcohol [yes/no]; the response might be the blood cholesterole level. The model assuming additive effects of smoking and drinking on the cholesterole levels would be y = m + a*xa + b*xb, where xa and xb are appropriately coded indicator variables for the two factors (x=0 for "no" and x=1 for "yes"). The coefficient m is the cholesterole level in a non-smoker non-drinker. The level of a smoker is changed (increased in this example) by a units, the level of a drinker is changed (also increased) by b units, and the level of a person that smokes and drinks is increased by a+b units. But if smoking and drinking have synergistic effects, the cholesterol levels of a person who smokes and drinks will be systematically larger than only a+b. The extra change can be accounted for by an interaction term in the model:, so that y = m + a*xa + b*xb + c*xa*xb. The coefficient c gives the additional change in cholesterol caused by the combination of smoking and drinking. One may interpret this as follows: a is the smoking effect (a difference "smoker vs.. non-smoker"), but this refers to non-drinkers. For drinkers, the effect is different - it is now a+c. The coefficient c is thus a difference between effects, which are differences in the response -> "difference in difference". You can also start the interpretation from b as the effect of drinking that referes to non-smokers, and this effect will be b+c in smokers.

If you drop one of the main effects, the meanings of the coefficients change. For instance, the model y = m + b*xb + c*xa*xb does not include the main effect for smoking. For non-drinkers, the model predicts the same response for smokers and non-smokers (y = m). It is a question of the subject matter whether or not such a model makes sense. I would not think that makes sense in this particular example, because I would excpect a smoke-effect also in non-drinkers.

If you drop both main effects, the model becomes y = m + c*xa*xb, and this would predict only two possible values: m+c for smokers+drinkers, and m otherwise. Such a model would possibly make sense if a change in the response can happen only when two things must come together, like two enzymes in a metabolic pathway.

@ Mario

Your question requires more explanation, about what interaction are you concerned. The model of your interest is much more interestiong, especially before you treat it with some packages:)

Regards, Joachim

Statistics don't mind what model you fit. Statistics assume that you defined your model in a way sensible for you, adequate for the problem, asnd that you will know the meanings of the coefficients.

Write the model down, look what the coefficients mean in your model. Do this for a model with and without lower-order terms. Start with very simple examples (only main effects and two-way interaction), simulate data and see how coefficients are estimated on that data.

You could envision a chemical reaction that requires a catalyst as a system where there is no main effect of the catalyst, but a strong interaction effect.

In biological systems it is more difficult to make such clear distinctions. So I have three factors A, B, and C. The main effects are all non-significant with p-values of 0.07. The pairwise interactions are all non-significant with p-values of 0.06. The three-way interaction is significant with a p-value of 0.0000000003.

1) a "non-significant" p-value does not mean that there was no effect. It probably means that the effect was too small to be detectable given your data and sample size.

2) How do you keep the ABC interaction yet claim that A had no effect? One possibility is that an unmeasured factor that needs to be in your model? Maybe A has no effect, except that it dilutes B and C? To account for this, you might need a different experimental design.

The problem might be simpler: you need more replicates/treatment.

OK, my question is what is the definition of interaction we are talking about. Is this a polynomial  where  y = a + bx1 + cx2 + dx3 + e x1*x2+ f x1*x3 + g x2*x3 + h x1*x2*x3? Or spmething similar to this dependence? Either  something completely different?

Thanks in advance for any clear answer. Joachim

The way that I have encountered the question (and the way I was assuming it was intended) was in ANOVA where you have main effects (A, B, C) and interactions. Low order interactions are all pairwise, and higher order would be three way, then four .... So the standard (from my background) is to say that you need to keep the model hierarchical. That is, if you have a B*C interaction, then you need to keep the main effects of B and C in the model. If you have an A*B*C interaction then you need all A, B, and C main effects and all pairwise interaction terms as well.

I was trying to figure out if there was any way that this would not make sense. It would require some theoretical model that explains why B should not have a main effect, but will influence the effects of A and C. I could not think of any biological examples, and I can think of several ways that the catalyst example could fail. One way is that the catalyst dilutes the other reactants.

Polynomial regression is different.

@Joachim: I don't get what you mean with "polynomial". As Timothy said: polynomial regression is different.

In a polynomial regression you use a polynomial of a metric predictor. The mean is the (very) special case of a polynomial regression of order 0: y = b0x0 = b0.The linear regression is the special case of a polynomial regression of order 1: y = b0 + b1x1 = b0 + b1x.The typical polynomial regression models contain higher powers of x, like y = b0 + b1x + b2x2 + b3x3 + ...

If the aim of the regression is to model, for instance, the quadratic relationship between speed and stopping distance of a car, or between the distance of a light source and the light intensity and such, then the linear term (b1x) is not required - it does not even make much sense to include this term in the model, so that an adequate model would look like y = b0 + b2x2. Under certain conditions it may be neccesary or known that for x=0 y must be 0, too, so that b0 should also be removed from the model.

As another example consider Poiseuille's law of a laminar liquid flow through a pipe of radius r. The flow volume is proportional to r4.  If you measure the flow y depending on the radius x (possibly to ftest the hypothesis that the flow is actually laminar), a model of the form y = b4x4 would be appropriate.

In general linear models the word "interaction" is usually understood as a "difference in diffrence" that allows to model synergistic or antagonistic effects (where the change in the response is more [or less, resp.] than the sum of the effects of the predictors considered separately.

Consider two dichotomous factors A and B, like A: person smokes at least 1 pk/d [yes/no]; B: person drinks at least 20 g/d alcohol [yes/no]; the response might be the blood cholesterole level. The model assuming additive effects of smoking and drinking on the cholesterole levels would be y = m + a*xa + b*xb, where xa and xb are appropriately coded indicator variables for the two factors (x=0 for "no" and x=1 for "yes"). The coefficient m is the cholesterole level in a non-smoker non-drinker. The level of a smoker is changed (increased in this example) by a units, the level of a drinker is changed (also increased) by b units, and the level of a person that smokes and drinks is increased by a+b units. But if smoking and drinking have synergistic effects, the cholesterol levels of a person who smokes and drinks will be systematically larger than only a+b. The extra change can be accounted for by an interaction term in the model:, so that y = m + a*xa + b*xb + c*xa*xb. The coefficient c gives the additional change in cholesterol caused by the combination of smoking and drinking. One may interpret this as follows: a is the smoking effect (a difference "smoker vs.. non-smoker"), but this refers to non-drinkers. For drinkers, the effect is different - it is now a+c. The coefficient c is thus a difference between effects, which are differences in the response -> "difference in difference". You can also start the interpretation from b as the effect of drinking that referes to non-smokers, and this effect will be b+c in smokers.

If you drop one of the main effects, the meanings of the coefficients change. For instance, the model y = m + b*xb + c*xa*xb does not include the main effect for smoking. For non-drinkers, the model predicts the same response for smokers and non-smokers (y = m). It is a question of the subject matter whether or not such a model makes sense. I would not think that makes sense in this particular example, because I would excpect a smoke-effect also in non-drinkers.

If you drop both main effects, the model becomes y = m + c*xa*xb, and this would predict only two possible values: m+c for smokers+drinkers, and m otherwise. Such a model would possibly make sense if a change in the response can happen only when two things must come together, like two enzymes in a metabolic pathway.

@ Timothy, @ Jochen

OK, I like your explanation since it CONFIRMS my guess:  the dependence with interaction IS polynomial of more that one variable. I named them  x1, x2, x3, Timothy named them A, B and C  and Jochen named them  xa, xb.  His formula

y = m + a*xa + b*xb + c*xa*xb

is exactly of polynomial form of total degree at most 2 and of degree  1  in each variable.Also my formula

y = a + bx1 + cx2 + dx3 + e x1*x2+ f x1*x3 + g x2*x3 + h x1*x2*x3

is polynomal of three variables, of total degree 3 and of degree  1  in each variable.

My question is, what is exactly the meaning of Timothy's  formula in the Jochen presentation if the variables are binary.

What I mean is as follows: for binary variables A,B,D the model described by  

D  =  a A  + b B  +  c A*B 

means that the probabilities (frequencies(?)) satisfy

P(D=1) = a P(A=1) + b P(B=1) + c P(A=B=1). 

Please tell me if I am right. If not, then please continue teaching me. In any case - of continual variables, of binary variables or of 3- and more-valued variables I can suggest some solutions. However, sorry, with regard to ANOVA no advice can be expected from me.

Best regards, Joachim

Good answers Timothy A Ebert and Jochen Wilhelm, especially the last sentence "Such a model would possibly make sense if a change in the response can happen only when two things must come together, like two enzymes in a metabolic pathway."  In some areas of research, the presence or absence of a minor factor results in dramatic outcome changes based on the synthetic (and sometimes synergistic) interactions of the factors.  Excellent conversation here, thanks for placing it in the questions section.

Cheers!

Rich

@Joachim,

 No, not even a little.

In a polynomial I can have a term that is A2, A3, or A4.

In ANOVA (=AOV, = analysis of variance) this is nonsensical. "A" cannot interact with itself, let alone interact with itself multiple times.

In polynomial regression I can have an A4B9C3 interaction term. In ANOVA the highest interaction term is the ABC interaction (assuming that I have three factors).

Thanks to all of you, the clarification is really good.

Sure, the model depends on what is my answer, this I had it clear. But in some places a read that there is a "computational" constraint when defining the model. 

I have found an explanation (see link below) on why mathematically is more recommended ot create a new variable enconding the interaction (e.g. x3=x1:x2) and run the model with this new variable (Y=a+dx3) than to run a model including the interaction and not the main effects (Y=a+cx1x2).

what do you think about that?

thanks again. I find this to be a very interesting discussion

https://stats.stackexchange.com/questions/27724/do-all-interactions-terms-need-their-individual-terms-in-regression-model

@ Timothy

Dear sir, your are not right with respect to the definition of polynomal. Simply, your dependence IS a multivariate polynomial with restricted set of non-zero coefficients (see my last but one answer above).  It is worth to add, that  function  ax +by + cxy  is the same as polynomal  f +  ax + by + cxy + dx^2 + ey^2   with  restriction d=e=f=0. However, OK, in the handbooks for programmers, in particular within software packages, these polynomials are called MPolynomials (afret multivariate polynomials). In pure mathematics, however, the adding of "multi-"  is used mainly for stressing their multidimensional domain;  consult, please:

http://dgtal.org/doc/0.8/modulePolynomial.html

http://www.siam.org/books/ot132/AMB_Slides_7_1_Polynomials.pdf

Moreover,  if  X  is a binary variable (like in case discussed here), with vales  0  or  1, only, then there is really senseless  to use  products of powers of  X  and  Y   because  X^n*Y^m=X*Y  for  all  n, m =1,2,...But, the result still is a polynomial of two variables.

Best regards, Joachim

@Joachim,

The (general linear) model does neither predict (or fit) frequencies nor probabilities. It predicts the values of a continous random variable, given the observed values of the predictors. The values of the "best" coefficients are determined by maximizing the likelihood of the observed data, given the model together with a stochastic part that is needed to judge what "best" means. The stochastic part in a general linear model sais that Y ~ Normal(µ,s²), where µ is given by your function f(X) relating the predictor values X to the response, and s² is a nuisance parameter.

If the response is binary, the the stochastic part is usually modelled as Y ~ Bernoulli(p), where p is given by your function f(X). To ensure that your function will return only permissive values for p (because p must be in [0,1]), one uses a link function g(p) = f(x) --> p = g-1(f(x)) that guarantees that p is in [0,1]. The common link function here is the g=logit. The intercept is then a log odds, and the main effect coefficients are log odds ratios (and frist-order interactions are differences between log odds ratios, etc.). Such models are called logistic (regression) models.

If the reponse would be a probability (or a propordion) itself (that is, Y can take all real values in [0,1]), it is often useful to assume that Y ~ beta(µ,φ), where again µ is given by f(x), usually again with the logit-link to make sure that µ is in [0,1]. Such models are called "beta (regression) models". The determination of the coefficients is still based on maximizing the likelihood of the observed data, but the solution can not be obtained within the usual (generalized) linear model framework.

To your example:

"for binary variables A,B,D the model described by  D  =  a A  + b B  +  c A*B"

If D (the response) is binary, we talk about a logistic model. One would not model P(D=1), but D (which can take the values 0 or 1). The model assumes that

D ~ Bernoulli( logit-1(aA+bB+cAB) )

Using this, the model can search those values for a, b, c, for which the likelihood of the observed D-values is maximal. Each observed D is associated with a set of values of the predictors A and B. When these are binary, A as well as B can also only be either 0 or 1 (usually - if not coded differently). The model does not care about that. For any observation one has a value for D and corresponding values for A and B. Chosing some (arbitrary) values of a, b, c you can calculate f = aA+bB+cAB and thus p = logit-1(f). Given this p you get the probability of getting your observed D from the Bernoulli distribution with parameter p: P(D) = pD(1-p)(1-D). This will be p if D=1 and (1-p) if D=0. The likelihood is the probability of getting all the observed D's together, so it is the product of the probabilities for all the observations. This likelihood refers to the chosen values of a,b,c. The fit is about finding those values for a,b,c, that gives the highest possible likelihood (the nice thing about linear models is that the solution can be found analytically, so there is no need to actually search the optimum).

Example:

the data is:

D = (0, 1, 1, 0)

A = (0, 1, 0, 1)

B = (0, 0, 1, 1)

Now let's chose a=2, b=-1.5 and c=0.2.

The values of f (in the same order as above) are then:

2*0 - 1.5*0 + 0.2*0*0 = 0

2*1 - 1.5*0 + 0.2*1*0 = 2

2*0 - 1.5*1 + 0.2*0*1 = -1.5

2*1 - 1.5*1 + 0.2*1*1 = 0.7

The values for the parameter p of the Bernoulli distribution are obtained by the inverse logit transformation of these values, that is:

logit-1(0) = 0.5

logit-1(2) = 0.881

logit-1(1.5) = 0.182

logit-1(0.7) = 0.668

The probabilities of the observed D's are:

Pr(D=0 | p=0.5) = 1-0.5 = 0.5

Pr(D=1 | p=0.811) = 0.881

Pr(D=1 | p=0.182) = 0.182

Pr(D=0 | p=0.668) = 1-0.668 = 0.332

The likelihood of observing all D's together is thus

0.5 * 0.881 * 0.182 * 0.332 = 0.027.

Searching other values for a,b,c reveals that a=20.57 b=20.57 and c=-61.70 give the highest possible likelihood, which is 0.5 in this example. For this combination of coefficients,

f = (0.00, 20.57, 20.57, -20.56),

p = logit-1(f) = (0.5, 1, 1, 0) (up to 9 decimal digits), and

Pr(D | p) = (0.5, 1, 1, 1)

so the likelihood is 0.5 * 1 * 1 * 1 = 0.5.

No other combination of a, b and c will give you a higher likelihood (for this data), therefore we say that these values are the "best".

Jochen, Thanks a lot for the detailed explanation of the usually applied model. Despite this probably very popular method, I still think that some algorithms joining the probabilities can be accepted as well. However, really, right now I do not have ready any particular mathematically justified example. Hopefully, I'll get it soon.

Again, cordial thanks for your lecture "for the beginners". I appreciate very much your contribtion, useful, probably for many othe RGaters.

All  best wishes to you, Joachim

Dear Joachim,

as I said: if D is real valued in [0,1], like for instance probability values, then a beta regression is an option. It is not relevant if the predictors are binary. The predictors can also be real values in [0,1] (these can also be proabilities, if you wish). In such a case I would use the logit-transformed values as predictors, so that the coefficients refer to a unit change in the log odds rather than to a unit change in probability.

To your example:

"P(D=1) = a P(A=1) + b P(B=1) + c P(A=B=1)"

Lets write P(D=1) = D', P(A=1) = A' and P(B=1) = B', so the model is

D'  =  a A'  + b B'  +  c A'*B'

Here, A', B' and D' are all real valued in [0,1]. The response D' can be modelled as a Beta-distributed random variable: D' ~ Beta(µ, φ) where µ ist given by logit-1(aA'+bB'+cA'B') and φ is a nuisance parameter (that is also estimated from the observed data). Inference on µ (and, hence on a, b and c) is done via the profile-likelihood. As before, the model does not care about what values A' and B' are allowed to take, as only observed (and,thus evidently allowed) values are used.

Not sure if this is what you meant...

I'll try again using numbers. The data are in the excel file. There are two independent variables: Date and trt. There is one dependent variable: Dep. In this example Date and trt are fixed effects in the ANOVA model.

I run ANOVA to see if date or treatment has an effect on the dependent variable.

I get a table like this

Type 3 Tests of Fixed Effects

Effect     Num DF    Den   DF    F Value    Pr>F

Date           3                72            17.07       F

Date           3                72            17.07       F

Date           3                72            17.07      

Jochen,

Thanks again for further details. My formal answer requires time, but some general lines can be formulated already as follows.

1. Yes, my first intention was to model the dependence by

D' = a A' + b B' + c A'*B',

treating this as a simple expression which can be meaningful for a rather wide range of the variables and the coefficients; the imagined way of determining the coefficients was the corresponding system of (linear) equations, where the observables are simply fractions of the sample with the possible 8 events (A=0 or 1, B=0 or 1, D=0 or 1), which should be treated e.g. by the least square method (with some additional improvments like weights etc.). Currently I have no idea if it can be formulated within a consistent mathematical model

2. Your excellent description of the way with the Dirichlet model of the a'priori distribution of the probabilities is very interesting and remainds me some Bayesian procedures. I will look at this topic much closer in few days. Perhaps the selfinvestigation will lead me to the algorithms you have presented within this current thread.

Best regards,

Joachim

@ Timothy

Let me share my deeply justified impression, that we are talking on polynomials from completely different worlds. A polynomial in the regression model is an ALGEBRAIC formula showing dependence on the independent numerical variables. Thus in your example the result of such sttistical analysis should be coefficients in the formula like Dep= a Date + b trt +  c Date*trt  or in one of the competing model quite another coefficients   of   Dep = A Date + B trt^2 +  c Date* trt^2.  The results must be different, at least based on the following reason:   if  the dimension of  trt is e.g meter, then   b  is of dimension  1/meter,  and  the dimension of   B  is  1/meter^2.

Contrarily, your undertanding of polynomial is herited from the symbolic notation of ANOVA, which means, that the result of statistical analysis are  power  of the influence of   the varibles listed  with  +   as a separator,  and  *  as an indicator, that e.g.   the third  independent varialble is the JOINT influence of the fixed variables (say  Date And trt)

Therefore our meaning have to be opposite!

Let me ask two questions: 

Q1. Shouldn't be the impact of joint influence of  Date*trt  greater than just of  each one separately ( in the first table )

Q2. Sholdn't be the impact factors in the tables  1  and  2 all the same?

Regards, Joachim

@Joachim,

In algebraic notation (exclusively) yes, all of these are polynomials. However, this gets very confusing when mixing algebraic notation and statistical notation.

So I have a categorical variable that can take states 1, 2, 3.

Since these are numeric I can take the square of these and get 1, 4, and 9.

However, the distance between 1 and 2 versus the distance between 1 and 4 has no bearing on the ANOVA model for a categorical variable. The only thing that matters is that I have three categories. It does not make any difference if those categories are "1", "2", and "3" or "red", "AK47", and "smurfs". The ANOVA model only cares about the existence of three categories. You can name them whatever makes you happy. The mathematical/statistical part will have the same outcome no matter what you choose.

Please do not assume anything about my understanding. I used ANOVA in one specific way to help show that the notation does not conform to the algebraic equivalent that would appear to be polynomial.

Notation: I use a + to separate variables as is standard practice. I can use either AB or A*B to indicate an interaction term. I have seen it both ways. I sometimes prefer the A*B format because sometimes variable names are more descriptive. So maybe I have a variable called Date and another variable called Time, I then have a variable called DateTime and the spelling checker has fits, while Date*Time is just fine. A*B is not the joint action of both A and B, it is the influence of factor A on the outcome of factor B.

A1. no. The two are independent. Just because factor A is significant has no bearing on whether factor A influences the outcome of factor B or to what degree. The notation is polynomial in form only, but has no correspondence to the algebraic concept.

A2. Curses! The problem with manually trying to type in things when you can't copy and paste. I edited the previous post. You are correct that the two tables should be exactly the same. Thank you for helping.

What might help you a great deal is to find a very small data set. Then find a basic statistics textbook, and try making all the calculations by hand. Alternatively, use a program like minitab, mathematica, or matlab but do everything manually, using only the program's linear algebra functions.

A final note: with the exception of acronyms like ANOVA, putting things in all upper case is yelling (or at least that is the custom here).

Dear Timothy, thanks a lot for your detailed explanation. I understand that "polnomial" wrt to ANOVA i ssenseless, I also undertand that the catecorial data will not change if their values are replaced  in a 1-1 way by arbitrary symbols.  And we now know, that your symbolic meaning of  + and  *  are different from the algebraic style. THAT'S OK. My problem was and is, that - frankly - I don't know the system ANOVA; I have supposed that this is one of packages which for estimating dependence/independence uses some kind of regression, only.  After some readings, in particular about AOV I have got closer to the algorithms of ANOVA. 

Moreover, and this is the basic source of my unsuitable answers at the beginning, in the original question there is no word about the algorithms used by Mario, only some effects are mentioned. Thus, I hope it is natural that I could  think that this is about effects measured by regression coefficients.

Finally, all my answers are motivated by a wish to explain first to the readers all I can say to the best of my current knowledge, adding also some original deriviations; next - with respect to myself, to get as much as possible to know within my interest from other participants of this particularly loyal and cooperating community of RG.

Best regards, Joachim

@Joachim,

I think (hope) that all RG contributors answer questions with the desire to help. Even though I think that Mario's question was about keeping the ANOVA model hierarchical, there is always the possibility that it was about a polynomial regression problem. The question did not provide enough detail to tell for certain, but the discussion has covered the two most likely interpretations for the question as asked.

Take care, Tim

Tim, thanks for understanding my point of view, and - hopefully - also my possibly really slightly different disires :)  Best to you, JoD

Hi Mario,

here a reference to one of apparently solid considerations of this issue (retaining lower level terms, model hierarchy), and many other important things.

Mundry R. 2014. Statistical issues and assumptions of phylogenetic generalized least squares. In: Modern phylogenetic comparative methods and their application in evolutionary biology (Garamszegi L, ed.), pp. 131–153. Berlin, Germany: Springer.

I attach a pdf just in case.

Evgeny