I performed a hierarchical multiple regression to study the effect of different variables, above and beyond the apriori selected more important variables. The R script used to perform this regression is listed below.  Data &lt;- read.csv("~/Desktop/data.csv", header=T,quote="") Data$catVar.f &lt;- factor(HMRdata$catVar) HRModel1 &lt;- lm(responseV ~ catVar.f, Data) HRModel2 &lt;- lm(responseV ~ catVar.f + v1, Data) HRModel3 &lt;- lm(responseV ~ catVar.f + v1+ v2, Data) HRModel4 &lt;- lm(responseV ~ catVar.f + v1 + v2 + v3, Data) summary (HRModel1) summary (HRModel2) summary (HRModel3) summary (HRModel4) anova (HRModel1, HRModel2, HRModel3, HRModel4)   My supervisor thought that the collinearity within the variables may pose a problem, and so asked to me perform a similar method called sequential regression, in which, instead of adding the new variables in the later regression analyses, I added residuals of these variables regressed against the previous, more important variables. The R script for this regression method is listed below –  RModel1 &lt;- lm(v1 ~ catVar.f, Data, na.action = na.exclude) v1residuals &lt;- resid(RModel1) RModel2 &lt;- lm(v2 ~ catVar.f + v1residuals, Data, na.action = na.exclude) v2residuals &lt;- resid(RModel2) RModel3 &lt;- lm(v3 ~ catVar.f + v1residuals + v2residuals, Data, na.action = na.exclude) v3residuals &lt;- resid(RModel3) SRModel1 &lt;- lm(responseV ~ catVar.f, Data) SRModel2 &lt;- lm(responseV ~ catVar.f + v1residuals, Data) SRModel3 &lt;- lm(responseV ~ catVar.f + v1residuals + v2residuals, Data) SRModel4 &lt;- lm(responseV ~ catVar.f + v1residuals + v2residuals + v3residuals, Data) summary (SRModel1) summary (SRModel2) summary (SRModel3) summary (SRModel4) anova (SRModel1, SRModel2, SRModel3, SRModel4)   To my surprise, the outputs (for summary and anova functions) of both the regression methods are exactly the same. Since I do not know the statistical principals behind these methods, I cannot figure out the reason for the same outputs. Can somebody please explain why both methods gave the same outputs? Also, is there any mistake in my approach?  EDIT - outputs from anova  Hierarchical Regression output  Model 1: responseV ~ catVar.f Model 2: responseV ~ catVar.f + v1   Model 3: responseV ~ catVar.f + v1 + v2 Model 4: responseV ~ catVar.f + v1 + v2 + v3 + v4 Model 5: responseV ~ catVar.f + v1 + v2 + v3 + v4 + v5 + v6 +v7 + v8          Res.Df    RSS Df Sum of Sq        F    Pr(&gt;F)     1     53 7324.2 2     52 2663.8  1    4660.3 111.3874 9.347e-14 *** 3     51 2057.7  1     606.1  14.4874 0.0004238 *** 4     49 1968.2  2      89.5   1.0699 0.3516141 5     45 1882.8  4      85.4   0.5103 0.7283774   Sequential regression output  Model 1: responseV ~ catVar.f Model 2: responseV ~ catVar.f + v1  residuals Model 3: responseV ~ catVar.f + v1residuals + v2residuals Model 4: responseV ~ catVar.f + v1residuals + v2residuals + v3residuals + v4residuals Model 5: responseV ~ catVar.f + v1residuals + v2residuals + v3residuals + v4residuals + v5residuals + v6residuals +v7residuals + v8residuals    Res.Df    RSS Df Sum of Sq        F    Pr(&gt;F)     1     53 7324.2  2     52 2663.8  1    4660.3 111.3874 9.347e-14 *** 3     51 2057.7  1     606.1  14.4874 0.0004238 *** 4     49 1968.2  2      89.5   1.0699 0.3516141 5     45 1882.8  4      85.4   0.5103 0.7283774
27828
Not sure if normalize is the correct word to use here, but I will try my best to illustrate what I am trying to ask. The estimator used here is least squares.  Suppose you have $y=\beta_0+\beta_1x_1$, you can center it around the mean by $y=\beta_0'+\beta_1x_1'$ where $\beta_0'=\beta_0+\beta_1\bar x_1$ and $x_1'=x-\bar x$, so that $\beta_0'$ no longer has any influence on estimating $\beta_1$.  By this I mean $\hat\beta_1$ in $y=\beta_1x_1'$ is equivalent to $\hat\beta_1$ in $y=\beta_0+\beta_1x_1$. We have reduced equation for easier least square calculation.  How do you apply this method in general? Now I have the model $y=\beta_1e^{x_1t}+\beta_2e^{x_2t}$, I am trying to reduce it to $y=\beta_1x'$.
56210
I know that one of the advantages of mixed models is that they allow to specify variance-covariance matrix for the data (compound symmetry, autoregressive, unstructured, etc.) However, lmer function in R does not allow for easy specification of this matrix. Does anyone know what structure lmer uses by default and why there is no way to easily specify it?
88114
I'm not used to using variables in the date format in R. I'm just wondering if it is possible to add a date variable as an explanatory variable in a linear regression model. If it's possible, how can we interpret the coefficient? Is it the effect of one day on the outcome variable?   See my  with an example what I'm trying to do.
17858
A few days ago, a psychologist-researcher of mine told me about his method to select variables to linear regression model. I guess it's not good, but I need to ask someone else to make sure. The method is:      Look at correlation matrix between all variables (including Dependent   Variable Y) and choose those predictors Xs, that correlate most with Y.     He didn't mention any criterion.  Q: Was he right?   [I think that this selection method is wrong, because of many things, like it's the theory that says which predictors should be selected, or even  omitted variable bias (OVB).]
48469
I am using a MANOVA test to compare nine different dependent variables (from neuropsychological and neuropsychiatric assessment) between three groups. The output shows a significant influence from GROUP on my variables ($p &lt; .001$).   Ofcourse, I am interested in how the three groups influence every dependent variable. I have studied Field's "Discovering Statistics Using IBM SPSS Statistics" chapter 16, and he states that the preferred post-hoc analysis is a discriminant analysis, because of the linear combination in which the dependent variables are related to group membership in a MANOVA. Discriminant analysis could account for this linear combination, so Field states.  Otherwise, I read some literature, on basic statistical sites, where is stated that I can use multiple univariate ANOVA's with Bonferroni correction, and use a post-hoc on these univariate ANOVA's when they are significant.   Which of these method's is better? That is, which will make my chance of a Type I or II error the least?
98436
First off, I'll say I am a biologist and new to the statistics side of things so excuse my ignorance  I have a data set that consists of a binary outcome and then a bunch of trinary explanatory variables that looks something like this:  head()  Category block21_hap1 block21_hap2 block21_hap3 block21_check 1        1            1            1            0             2 2        1            2            0            0             2 3        1            1            0            1             2 4        1            1            0            1             2 5        1            1            1            0             2 6        1            1            1            0             2   A quick summary of the data  summary() Category block21_hap1 block21_hap2 block21_hap3 block21_check  1:718    0:293        0:777        0:1026       2:1467         0:749    1:709        1:577        1: 390                               2:465        2:113        2:  51     and another summary grouped by outcome levels  by(hap.ped.final, hap.ped.final$Category, summary) hap.ped.final$Category: 1  block21_hap1 block21_hap2 block21_hap3 block21_check  0:146        0:374        0:518        2:718          1:336        1:286        1:174                       2:236        2: 58        2: 26                      ----------------------------------------------------------------------------  hap.ped.final$Category: 0  block21_hap1 block21_hap2 block21_hap3 block21_check  0:147        0:403        0:508        2:749          1:373        1:291        1:216                       2:229        2: 55        2: 25             So I am trying to run logistic regression on this data. When I do this:  fit = glm(Category~ block21_hap1 + block21_hap2 + block21_hap3, data = hap.ped.final ,family = "binomial") summary(fit)  Deviance Residuals:     Min      1Q  Median      3Q     Max   -1.301  -1.177   1.059   1.177   1.200    Coefficients: (1 not defined because of singularities)                              Estimate Std. Error z value Pr(&gt;|z|) (Intercept)                 -0.039221   0.280110  -0.140    0.889 hap.ped.final$block21_hap11  0.123555   0.183087   0.675    0.500 hap.ped.final$block21_hap12  0.009111   0.295069   0.031    0.975 hap.ped.final$block21_hap21 -0.084334   0.183087  -0.461    0.645 hap.ped.final$block21_hap22 -0.013889   0.337468  -0.041    0.967 hap.ped.final$block21_hap31  0.201113   0.183087   1.098    0.272 hap.ped.final$block21_hap32        NA         NA      NA       NA  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 2033  on 1466  degrees of freedom Residual deviance: 2028  on 1461  degrees of freedom AIC: 2040  Number of Fisher Scoring iterations: 3   So I don't really know what a singularity is or what's going wrong here that is throwing up NA's as a result of my analysis. Is it my data, or what I'm doing to it. I tried googling the warning (or whatever you might call it) and I got some pages talking about collinearity and multilinearity, which I do not understand at all.  Again, sorry for lack of knowledge here. I wish I had done more maths in undergrad.
58570
There is a significant 3-way interaction in a data-set I'm working with.  The interaction involves both categorical and quantitative variables.  I have been pointed towards simple slopes and this  but I find the explanations lacking. I only have a basic background in statistics, and Googling for other examples has not been very helpful.   Any insight as to how to begin to understand this interaction will be very welcome, especially using R rather than SPSS. Thank you very much in advance, from a lower-year college student who has suddenly found himself over his head!
92654
Suppose I have a linear regression model Y ~ b1*x1 + b2*x2 with the sample data:          x1      x2        y 0  -0.7251 -0.6033  10.2508 1  -1.3148  3.7067  10.3054 2  -1.2818 -5.8136   7.2010 3   0.6154  3.6336  11.2807 4   1.4887 -2.8337  10.6299 5  -2.0224 -0.3931   7.4502 6   0.5419 -1.7807  10.7339 7  -1.0745 -5.0296   5.8734 8   3.2116 -0.4681  11.7184 9   1.6767 -4.6040   9.0145 10 -1.3355 -2.6874   6.5014 11  0.8296  2.6793  11.8295 12 -2.8647  4.2916   7.8414 13  1.8445 -3.3924   9.6520 14 -0.1036  2.4273  11.9766 15  1.0886 -3.5636  10.6023 16 -3.3026 -6.9837   4.9962 17  1.4050 -3.8453  10.4630 18 -0.3111 -1.4278   9.2283 19 -3.5454  2.4960   7.3720   And using the following R code to get the Analysis of Variance Table:  &gt; df &lt;- read.csv('f:/data.csv', header=T, sep=",") &gt; model_All &lt;- lm(y ~ x1 + x2, data=df) &gt; anova(model_All) Analysis of Variance Table  Response: y           Df Sum Sq Mean Sq F value    Pr(&gt;F)     x1         1 44.440  44.440  58.277 6.865e-07 *** x2         1 26.640  26.640  34.934 1.711e-05 *** Residuals 17 12.964   0.763                       --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   And I can calculate the F value of hypothesis test H0： b2=0 via:  &gt; model_Drop2 &lt;- lm(y ~ x1, data=df) &gt; SSR_All = sum(resid(model_All)^2) &gt; SSR_Dp2 = sum(resid(model_Drop2)^2) &gt; F_x2 = ((SSR_Dp2 - SSR_All) / 1) / (SSR_All / (20 - 3)) &gt; F_x2 [1] 34.93434   We can see that F_x2 is equal to F value of x2 in Analysis of Variance Table. However, the value F_x1 calculated as follow is not equal to F value of x1, I can't understand why:  &gt; model_Drop1 &lt;- lm(y ~ x2, data=df) &gt; SSR_Dp1 = sum(resid(model_Drop1)^2) &gt; F_x1 = ((SSR_Dp1 - SSR_All) / 1) / (SSR_All / (20 - 3)) &gt; F_x1 [1] 70.30392   So, what is the mathematics definition of the F value in table?  EDIT The result of the features reversed as commented by @Michael M:  &gt; model_All2 &lt;- lm(y ~ x2 + x1, data=df) &gt; anova(model_All2) Analysis of Variance Table  Response: y           Df Sum Sq Mean Sq F value    Pr(&gt;F)     x2         1 17.468  17.468  22.907 0.0001718 *** x1         1 53.612  53.612  70.304 1.914e-07 *** Residuals 17 12.964   0.763                       --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
40756
I found a reference in an article that goes like:     According to Tabachnick &amp; Fidell (1996) the independent variables with   a bivariate correlation more than .70 should not be included in   multiple regression analysis.   Problem: I used in a multiple regression design 3 variables correlated >.80, VIF's at about .2 - .3, Tolerance ~ 4- 5. I cannot exclude any of them (important predictors and outcome). When I regressed the outcome on the 2 predictors which correlated at .80, they remained both significant, each predicted important variances, and these same two variables have   the largest part and semipartial correlation coefficients amongst all 10 variables included (5 controls).   Question: Is my model valid despite high correlations? Any references greatly welcomed!    Thank you for the answers!  I did not use Tabachnick and Fidell as a guideline, I found this reference in an article dealing with high collinearity amongst predictors.   So, basically, I have too few cases for the number of predictors in the model (many categorical, dummy coded control variables- age, tenure, gender, etc) - 13 variables for 72 cases. The Condition Index is ~ 29 with all the controls in and ~23 without them (5 variables).   I cannot drop any variable or use factorial analysis to combine them because theoretically they have sense on their own. It is too late to get more data. Since I am conducting the analysis in SPSS perhaps it would be best to find a syntax for ridge regression (although I haven't done this before and interpreting the results would be new to me).   If it matters, when I conducted stepwise regression, the same 2 highly correlated variables remained the single significant predictors of the outcome.  And I still do not understand if the partial correlations which are high for each of these variables matter as an explanation for why I have kept them in the model (in case ridge regression can't be performed).   Would you say the "Regression diagnostic: identifying influential data and sources of collinearity / David A. Belsley, Edwin Kuh and Roy E. Welsch, 1980" would be helpful in understanding multicollinearity? Or might other references be useful?
23107
I am currently running a multiple linear regression, and I am bit confused in regards to how to properly add interaction terms to the model by hand. All of the variables I am using are continuous and have different scales and units.  So far, the way I have done it has been to   Standardize the observations for each variables Multiply corresponding standardized values from specific variables to create the interaction terms and then add these new variables to the set of regression data Run the regression   Is this the correct way to go about doing this? Should I standardize the interaction term variables also after calculating the 'raw' terms?
113081
Background:  Note: My data set and R code are included below text I wish to use AIC to compare two mixed effects models generated using the lme4 package in R.  Each model has one fixed effect and one random effect.  The fixed effect differs between models, but the random effect remains the same between models.  I've found that if I use REML=TRUE, model2 has the lower AIC score, but if I use REML=FALSE, model1 has the lower AIC score. Support for using ML: Zuur et al. (2009; p. 122) suggest that &quot;To compare models with nested fixed effects (but with the same random structure), ML estimation must be used and not REML.&quot;  This indicates to me that I ought to use ML since my random effects are the same in both models, but my fixed effects differ.  [Zuur et al. 2009.  Mixed Effect Models and Extensions in Ecology with R.  Springer.] Support for using REML: However, I notice that when I use ML, the residual variance associated with the random effects differs between the two models (model1 = 136.3; model2 = 112.9), but when I use REML, it is the same between models (model1=model2=151.5).  This implies to me that I ought instead to use REML so that the random residual variance remains the same between models with the same random variable. Question: Doesn't it make more sense to use REML than ML for comparisons of models where the fixed effects change and the random effects remain the same?  If not, can you explain why or point me to other literature that explains more? # Model2 &quot;wins&quot; if REML=TRUE: REMLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = TRUE) REMLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = TRUE) AIC(REMLmodel1,REMLmodel2) summary(REMLmodel1) summary(REMLmodel2)  # Model1 &quot;wins&quot; if REML=FALSE: MLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = FALSE) MLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = FALSE) AIC(MLmodel1,MLmodel2) summary(MLmodel1) summary(MLmodel2)  Dataset: Response    Fixed1  Fixed2  Random1 5.20    A   A   1 32.50   A   A   1 6.57    A   A   2 24.77   A   B   3 41.69   A   B   3 34.29   A   B   4 1.80    A   B   4 10.00   A   B   5 15.56   A   B   5 4.44    A   C   6 21.65   A   C   6 9.20    A   C   7 4.11    A   C   7 12.52   B   D   8 0.25    B   D   8 27.34   B   D   9 11.54   B   E   10 0.86    B   E   10 0.68    B   E   11 4.00    B   E   11
107346
In the context of regressions, it seems a convention that the HAC estimator should be applied when the residual is serially correlated. But isn't the presence of residual autocorrelations an indication that the model misses lagged dependent variables (or MA terms for that matter)?  For example, suppose one is investigating the impact of variable $x_t$ on variable $y_t$ in a time series context. A simple model is $y_t = a + b x_t + \varepsilon_t$. It is often the case that the resulting error term $\varepsilon_t$ is serially correlated. One response to this problem could be just applying the HAC estimator to correct the standard errors of the parameter estimates, while arguing that the parameter estimates are still consistent in the presence of residual serial correlation.   However, in this case, one could often add lags of $y_t$ to the model to remove the autocorrelations in the residuals, e.g. $y_t = a + b x_t + c y_{t-1} + \varepsilon_t$. Since the lag $y_{t-1}$ could be correlated with the regressor $x_t$, either by economic reasoning ($x_t$ is often serially correlated) or by sample correlation between the lag $y_{t-1}$ and the regressor $x_t$, it would be simply wrong to ignore the lag $y_{t-1}$ while just correcting the standard errors by HAC.  This makes me wonder why one would ever use HAC estimator in regressions? At most, use the White estimator to correct for heteroskedasticidy if it is present. Otherwise, whenever there is residual autocorrelation, investigate adding AR or MA terms but never simply use HAC estimator.
37344
I watched several videos on linear regression, mainly from . As I have no background in statistics, I thought this was a good way to get an idea of the topic. However I'm currently writing a bachelor thesis on Structural Equation Modeling and I want to get a deeper understanding of Regression Analysis.  The Book should roughly include these topics:   linear least squares regression variance, covariance regression coefficient coefficient of determination residual analysis (esp. leverage effect) multiple regression   I especially enjoyed the videos on khan academy, of the proof of minimizing the squared error to the regression line or illustrations like this one from Wikipedia:    I don't like it when books only present a formula and some common rules etc. without any further explanation of these.  F.e. Residuals should be normally distributed.  However I'm also not comfortable with strictly mathematical books, so I'm searching for something in between. A plus would be example code in R, but definitely not a must have.
95307
My question concerns the setup of a mixed model, which is designed to investigate the interaction between two different interventions over time. I will provide a sample data-set in R, but I am generally happy about any theoretical advice.   This is the scenario. It's a repeated measures dataset with two different interventions (food and training), each of them have two levels. All participants underwent each combination of conditions. Hormone samples were collected over the course of the day. Age and BMI are included as covariates.  This is a reproducible data set:  ID &lt;- rep(rep(letters[1:16],each=8),4) TIME &lt;- rep(rep(paste(sprintf("%02i",9:16),"00",sep=":"),16),4) training &lt;- rep(rep(rep(c("T1","T2"),each=8),each=16),2) food &lt;- c(rep(rep(c("F1","F2"),each=8),each=16),rep(rep(c("F2","F1"),each=8),each=16)) hormone &lt;- rnorm(n = 512,mean=5,sd=2) BMI &lt;- as.numeric(rep(rep(sample(18:35,size = 16,replace = T),each=8),4)) age &lt;- as.numeric(rep(rep(sample(40:60,size = 16,replace=T),each=8),4)) DF &lt;- data.frame(ID,TIME,training,food,hormone,BMI,age)  indNA &lt;- which(DF$hormone%in%sample(DF$hormone,23)) indNABMI &lt;- which(DF$ID%in%sample(DF$ID,2)) DF$BMI[indNABMI] &lt;- NA DF$hormone[indNA] &lt;- NA  DF$BMI &lt;- scale(DF$BMI) DF$age &lt;- scale(DF$age)   Now in order to model the effect of Time, I have found various solutions, as for example here:  where it is suggested to include daytime as a circular variable. However, my sampling does not cover the entire day, so am I right in assuming this is not an option for me?  I would then argue for including it as an ordered factor variable as I did here:  DF$TIME &lt;- ordered(DF$TIME) str(DF$TIME) Ord.factor w/ 8 levels "09:00"&lt;"10:00"&lt;..: 1 2 3 4 5 6 7 8 1 2 ...   Next, I am unsure about how to include this in my model. As I said, my research question is about the interaction of the 2 intervention types, so I would come up with the following model:  m1 &lt;- lmer(hormone~training*food*TIME+age+BMI+(1+training*food*TIME|ID),DF)   which I would test against the according NULL-Model: (if you see any flaws in the model setup to begin with, please let me know)  m01 &lt;- lmer(hormone~training*TIME+food*TIME+age+BMI+(1+training*TIME+food*TIME|ID),DF)   However I do not have enough observations to test this model and I therefore get this error-message:  Error: number of observations (=427) &lt;= number of random effects (=448) for term (1 + training * food * TIME | ID); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable   and the NULL-Model does not converge either.  How can I alternatively include the effect of time in my random effects? I would especially like to include it, since hormones display diurnal fluctuations, which might again differ between participants but potentially also between conditions.  In addition, this dataset includes data from a pharmacological test, where I expect the hormone levels to rise and drop again over time. The shape would thus resemble a quadratic shape rather than a linear one?  Would I in that case, include the interaction of conditions and quadratic time (as a numerical predictor)?
119620
Suppose, I have performed multiple regression analysis on the following data set where X1 and X2 are independent variables and Y is the dependent variable.    And achieved the following multiple regression analysis table from where it is clear that X2 is not correlated with Y as the p-value is greater than 0.05 which is 0.108 for X2 variable     So, if p value is greater than 0.05 for a variable, then why it is considered in regression equation  Y = 2.709 + 0.763(X1) + 0.463(X2) which has been shown in any online multiple regression equation calculator.  My question is Why it is not like the following Y = 2.709 + 0.763(X1) omitting the [0.463(X2)] part from the equation, as X2 is not correlated with Y  Please help to clarify the conceptual idea....
73396
I want to estimate a multilevel multinomial logit model but I am struggling with the terminology and notation used by the R-package MCMCglmm. There is documentation available in form of a  and . According to the article, "the [MCMCglmm] syntax used to specify the model closely follows that used by asreml" (Hadfield 2010: 6). Consequently, some abbreviations like us="unstructured" and idh="independent heterogeneous" are never stated explicitly in the MCMCglmm documentation, although they are well presented in formulas. Unfortunately, I am only acquainted with the notation for hierarchical models in the social sciences and I have no experience with ASReml.  The examples of multinomial multilevel models in the MCMCglmm course notes are subsumed under multiresponse models and the course notes (p. 91) show that trait is used as identifier when the dataset is reshaped for estimating multiresponse models.   Is it necessary to use trait in the fixed (or random) formula also for models with a nominal outcome variable like employment status or gender (in contrast to a real multiresponse variable)? Does the inclusion of trait change the estimated formula or is it only a signifier for MCMCglmm to expect reshaped data? Can I define the error structure without using trait, e.g. by specifying the error structure for coefficients from the fixed part?   Here is an (almost working) example code using a typical social science survey dataset:  library(nnet) library(MCMCglmm)  # Load dataset from TraMineR-package data(mvad, package="TraMineR")  ### Create variable "region" for 2nd level attach(mvad) mvad$region[Belfast=="yes"] &lt;- "Belfast"     mvad$region[N.Eastern=="yes"] &lt;- "N.Eastern" mvad$region[Southern=="yes"] &lt;- "Southern"     mvad$region[S.Eastern=="yes"] &lt;- "S.Eastern" mvad$region[Western=="yes"] &lt;- "Western" detach(mvad)  # Add identifier for complete cases  mvad$complete &lt;- complete.cases(mvad[, c("Jul.93", "male", "Grammar",                                           "gcse5eq", "funemp", "region")])  # Multinomial logit model without random effects using nnet-package multinom(Jul.93 ~ male + Grammar + gcse5eq + funemp + region,           data = mvad[mvad$complete==TRUE, ])  # Multinomial logit models with random effects using MCMCglmm mod1 &lt;- MCMCglmm(Jul.93 ~ trait + male + Grammar + gcse5eq + funemp,                 random = ~ region, rcov = ~ us(trait):units,                 family = "categorical", verbose=FALSE,                 data=mvad[mvad$complete==TRUE, ]) # Including trait and defining rcov gives an estimable command  mod2 &lt;- MCMCglmm(Jul.93 ~ male + Grammar + gcse5eq + funemp,                  random = ~ region, rcov = ~ us(trait):units,                  family = "categorical", verbose=FALSE,                  data=mvad[mvad$complete==TRUE, ]) # Command works as well, but is trait now only used for the error terms?  mod3 &lt;- MCMCglmm(Jul.93 ~ male + Grammar + gcse5eq + funemp,                  random = ~ region,                  family = "categorical", verbose=FALSE,                  data=mvad[mvad$complete==TRUE, ]) # Command does not work as definition of error structure using e.g. "us(trait):units" is needed.    Reference: Hadfield, J.D. (2010) 'MCMC Methods for Multi-Response Generalized Linear Mixed Models: The MCMCglmm R Package', Journal of Statistical Software 33(2): 1–22.
41432
As I was working my way through my book about statistics, I came across the topic of linear regression. During the chapter the author begins with explaining that you want to minimize the residuals in order to make your y = a + bx as good a fit as possible: I do understand this, but halfway the chapter all of a sudden the residuals change into sum squares of residuals. Why is this done? I have been googling, but couldn't find the right answer. Who would like to help me understand why the sums of squares of the residuals is used instead of just the sums of residuals?  Kind regards, Bas