I'm fitting a linear model to a dataset. I've reduced the model from 7 variables to 3 and as part of that I've transformed one of my variables (wt) in order to get a better fit.  For my final model, the summary is:  summary(lm22) Call: lm(formula = mpg ~ poly(wt, 2) + year + origin, data = redautoi)  Residuals:     Min      1Q  Median      3Q     Max  -6.4126 -1.5029 -0.0003  1.5994  7.6032   Coefficients:                Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)   -33.94088    2.84630 -11.925  &lt; 2e-16 *** poly(wt, 2)1 -100.09492    3.14080 -31.869  &lt; 2e-16 *** poly(wt, 2)2   27.18431    2.61085  10.412  &lt; 2e-16 *** year            0.74125    0.03736  19.839  &lt; 2e-16 *** origin          0.55447    0.20314   2.729  0.00665 **  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 2.477 on 367 degrees of freedom Multiple R-squared:  0.8842,    Adjusted R-squared:  0.883  F-statistic: 700.7 on 4 and 367 DF,  p-value: &lt; 2.2e-16   My question is how do I interpret these estimates? Especially  poly(wt,2)1 and poly(wt,2)2  Substituting values from the original dataset using the above estimates gives a wild fitted value. (but using fitted(lm22) appears to give reasonable estimates)
84506
I'm trying to create a second order polynomial fit to some data I have. Let's say I plot this fit with ggplot():  ggplot(data, aes(foo, bar)) + geom_point() +         geom_smooth(method="lm", formula=y~poly(x, 2))   I get:    So, a second order fit works quite well. I calculate it with R:  summary(lm(data$bar ~ poly(data$foo, 2)))   And I get:  lm(formula = data$bar ~ poly(data$foo, 2)) # ... # Coefficients: #                     Estimate Std. Error t value Pr(&gt;|t|)     # (Intercept)         3.268162   0.008282 394.623   &lt;2e-16 *** # poly(data$foo, 2)1 -0.122391   0.096225  -1.272    0.206 # poly(data$foo, 2)2  1.575391   0.096225  16.372   &lt;2e-16 *** # ....   Now, I would assume the formula for my fit is:  $$ \text{bar} = 3.268 - 0.122 \cdot \text{foo} + 1.575 \cdot \text{foo}^2 $$  But that just gives me the wrong values. For example, with $\text{foo}$ being 3 I would expect $\text{bar}$ to become something around 3.15. However, inserting into above formula I get:    $$ \text{bar} = 3.268 - 0.122 \cdot 3 + 1.575 \cdot 3^2 = 17.077 $$  What gives? Am I incorrectly interpreting the coefficients of the model?
106714
After searching here:    I am still curious about this.  Is it possible to derive the standard error of a regression coefficient from summary data alone?   E.g., assume we are given the following variance-covariance matrix.  $\begin{bmatrix}Var(X) &amp; Cov(X,Y)\\Cov(X,Y) &amp; Var(Y)\end{bmatrix}$  We can derive the regression coefficient $\beta_{XY} = Cov(X,Y) / Var(X)$.  Given a specific n, is it possible to derive the standard error of $\beta$ as well? If so, which formula is being used? It appears that all formulas for regression standard errors that I could find assume that you know the variance of residuals of the regression, which we don't know from summary data alone.  ---UPDATE--- Based on feedback by Greg and whuber, I am now at a point, where I know that  $\sigma^2(b) = \sigma^2(X'X)^{-1}=\sigma^2(Var(X)(n-1))^{-1}$,  but typically matrix formulas for the standard error of $b$ then note that $\sigma^2$, the unknown variance of the errors, is estimated with MSE.   Based on comments by whuber, sum of squares residuals is:  $Y'Y - bX'X$, which according to Greg can be spelled out as  $Var(Y)(n-1) - b(Var(x)(n-1))$  And to go from here to MSE, I think (?) all that is left is divide by n-2,  so $(Var(Y)(n-1) - b(Var(x)(n-1))) / (n-2)$  And plugging this all in yields:  $\sigma^2(b) =((Var(Y)(n-1) - b(Var(X)(n-1))) / (n-2))\times(Var(X)(n-1))^{-1}$  I will try to see if this works out in R, and will report back. If anyone sees something blatantly wrong, I appreciate the heads up.
111093
I built a logistic model with 7 continuous independent variable and 8 categorical (dummy) variables in R. I tested for the linearity assumption using the Box Tidwell test and all my variables show high significance.  I am wondering now, if the assumption is violated what it would mean for the coefficient estimates and if I still can run the model with caution.   I read about the approach of using regression splines but I am confused how I would determine those and interpret the model then.
41967
In R, I'm running a linear model of the form:  fit &lt;- lm(cbind(Y1,Y2,Y2) ~ X1 + X2 + X3, data=regression_table);   Then I'm using predict to generate predictions;  predicted_vals &lt;- predict(fit, test_data);   Does any one have suggestions on how I should go about evaluating the training error of this multivariate model? I understand the notion of cross validation for univariate multiple models, though I'm not sure how to go about go about doing something similar for a multivariate model.
16998
I have a set of experimental data points. I performed the measurements in triplicate, for each of the point of the data set. Therefore, I can draw each data point with the standard deviation of each triplicate. See the picture attached.  In experimental sciences, it is common to report a value with its standard deviation. Ex: a mean, +/- the std.  I can calculate a linear regression for the data set. If I have the equation of the linear regression, I can calculate x for any y. Let's take y=50 -&gt; x=11.69  Now, is there a way to evaluate the "dispersion" of this extrapolated point ? Something like 11.69 +/- something.  I know it should be the other way around, like 50 +/- something for x=11.69, but then I could use the equation to transform it to x.  Basically what I'm asking: is there a global "standard deviation" for a complete linear regression ?  EDIT: When I say "any y", I mean that y will not be an experimental value. I choose it to be 50.
17738
I'm using R to run some logistic regression. My variables were continuous, but I used cut to bucket the data. Some particular buckets for these variables always result in dependent variable being equal to 1. As expcted, the coefficient estimate for this bucket is very high, but the p-value is also high. There are about ~90 observations in either these buckets, and around 800 total observations, so I don't think it's a problem of sample size. Also, this variable should not be related to other variables, which would naturally reduce their p-values.  Are there any other plausible explanations for the high p-value?  Example:  myData &lt;- read.csv("application.csv", header = TRUE) myData$FICO &lt;- cut(myData$FICO, c(0, 660, 680, 700, 720, 740, 780, Inf), right = FALSE) myData$CLTV &lt;- cut(myData$CLTV, c(0, 70, 80, 90, 95, 100, 125, Inf), right = FALSE) fit &lt;- glm(Denied ~ CLTV + FICO, data = myData, family=binomial())   Results are something like this:  Deviance Residuals:       Min        1Q    Median        3Q       Max   -1.53831  -0.77944  -0.62487   0.00027   2.09771    Coefficients:                 Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)     -1.33630    0.23250  -5.747 9.06e-09 *** CLTV(70,80]     -0.54961    0.34864  -1.576 0.114930     CLTV(80,90]     -0.51413    0.31230  -1.646 0.099715 .   CLTV(90,95]     -0.74648    0.37221  -2.006 0.044904 *   CLTV(95,100]     0.38370    0.37709   1.018 0.308906     CLTV(100,125]   -0.01554    0.25187  -0.062 0.950792     CLTV(125,Inf]   18.49557  443.55550   0.042 0.966739     FICO[0,660)     19.64884 3956.18034   0.005 0.996037     FICO[660,680)    1.77008    0.47653   3.715 0.000204 *** FICO[680,700)    0.98575    0.30859   3.194 0.001402 **  FICO[700,720)    1.31767    0.27166   4.850 1.23e-06 *** FICO[720,740)    0.62720    0.29819   2.103 0.035434 *   FICO[740,780)    0.31605    0.23369   1.352 0.176236     --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 1037.43  on 810  degrees of freedom Residual deviance:  803.88  on 798  degrees of freedom AIC: 829.88  Number of Fisher Scoring iterations: 16   FICO in the range [0, 660) and CLTV in the range (125, Inf] indeed always results in Denial = 1, so their coefficients are very large, but why are they also "insignificant"?
13759
Suppose we have regression problem:  $$Y = \beta_{0}+\beta_{1}X_{1} + \beta_{2}X_{2} + \epsilon \text{, } \epsilon \sim \mathcal{N}(0,1),$$  where $X_1 \sim U(0,1)$,$X_{2} \sim U(1,2)$, and we suppose our model with regression function  $$m:=E[Y|X_1,X_2] = \beta_{0}+\beta_{1}X_{1} + \beta_{2}X_{2}. $$  We apply OLS estimation of coefficient and receive estimate of regression coeficeint. Let's do this in $\verb|R|$:  X1 &lt;- runif(1000,0,1) X2 &lt;- runif(1000,1,2) Y &lt;- X1 + X2 + rnorm(1000,0,1) model &lt;- lm(Y ~ X1 + X2) summary(model)  Coefficients:              Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)     0.1401     0.1828   0.766    0.444     X1              1.0250     0.1150   8.909  &lt; 2e-16 *** X2              0.9247     0.1139   8.119 1.37e-15 *** ---  Signif. codes:   0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 1.026 on 997 degrees of freedom Multiple R-squared:  0.1276,    Adjusted R-squared:  0.1259  F-statistic: 72.94 on 2 and 997 DF,  p-value: &lt; 2.2e-16   Next, we want to test, if the regression coeficient $\beta_2$ is significant or not, so here is T-test for this:  $$H_{0}: \beta_{2} = 0$$ $$H_{1}: \beta_{2} \neq 0$$  Aproppiate test statistic is:  $$T_n = \frac{\hat{\beta_{2}}}{S.E(\hat{\beta_{2}})} \sim t_{n-r},$$  and p-value:$2min\{CDF_{t}(t_{0}), 1 - CDF_{t}(t_0)\}$, where $CDF_t$ is cummulative distribution function o $t$-distribution with $n-r$, in ou case $1000-3$ degrees of freedom. Let's suppose, that we also made test about submodel:  $$H_{0} Y \sim X_{1} \text{ holds} (M^{0})$$ $$H_{1} Y \sim X_{1} + X_{2} \text{ holds} (M).$$  Test statistics is $F = \frac{\frac{SSe^0 - SSe}{r - r_{0}}}{MSe} \sim F_{r-r_0,n-r}$ and p-value:$1-CDF_f(f_0)$, where $CDF_f$ is cummulative distribution function of F-distribution with $r-r_0$ and $n-r$ degrees of freedom. In $\verb|R|$ let's make an test about submodel:  m0 &lt;- lm(Y ~ X1) anova(m0,model)  Analysis of Variance Table  Model 1: Y ~ X1  Model 2: Y ~ X1 + X2     Res.Df    RSS Df Sum of Sq      F    Pr(&gt;F)      1    998 1118.2                                    2    997 1048.8  1    69.352 65.926 1.373e-15 ***  --- Signif. codes:   0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Now there is a question, i know that if $X \sim t_{n-r}$, then $X^2 \sim F_{1,n-r}$. So if we squared t-value in our summary table, in line where regression coeficient $\hat{\beta}_2$ is, we geta f-value in output in test of submodel, so that's okay, but we have also same p-values, but p-value in summary table is computed from t-distribution, not squared distribution, and this numbers implies that $2min\{CDF_{t}(t_{0}), 1 - CDF_{t}(t_0)\} = 1-CDF_f(f_0)$, but i think this is not really true. What am i missing, how it is possible, that p-value which are computed from two different distributions are same ? (and this is not only for this case, but for every liner regression problems). Why p-values are same ? Please help.  In all cases i asuume normal linear model, that's why i assume those distribution of test statistics under $H_0$ hypothesis.
58523
I want to be able to calculate the confidence interval from the estimated coefficient and respective standard errors.  I have a linear regression model which can be summarized (in R):  summary(fit1)  Call: lm(formula = bwt ~ height + weight + parity, data = data)  Residuals:     Min      1Q  Median      3Q     Max  -66.913 -10.624   0.991  10.979  55.621   Coefficients:             Estimate Std. Error t value Pr(&gt;|t|)     (Intercept) 31.19217   13.56879   2.299   0.0217 *   height       1.24964    0.23083   5.414 7.48e-08 *** weight       0.06781    0.02823   2.402   0.0164 *   parity1     -1.83309    1.19838  -1.530   0.1264     --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 17.9 on 1170 degrees of freedom Multiple R-squared:  0.04898,   Adjusted R-squared:  0.04654  F-statistic: 20.08 on 3 and 1170 DF,  p-value: 1.071e-12   With this model I can calculate the respective confidence intervals:  &gt; confint(fit1)                   2.5 %     97.5 % (Intercept)  4.57029503 57.8140351 height       0.79676227  1.7025207 weight       0.01243198  0.1231932 parity1     -4.18429933  0.5181151   I would expect the intervals of the predictor height to be given by  $$ 1.24964 \pm (1.96*0.23083) = [0.7972132,1.702067] $$  where 1.24964 is the estimated value for the coefficient and 0.23083 is the standard error for this coefficient. The numbers are close but not quite the same.  What am I doing wrong?
79370
I have carried out a linear regression. The plot below shows the distribution of the response variable:    I believe the response variable is beta distributed, therefore virtually the exact opposite of normally distributed. However, when including all my predictor variables in the linear regression, the residuals turn out to be quite normally distributed, as shown in this plot:    Has my model satisfied the assumptions of linear regression? Might there be a better model to use?
101841
I have a linear regression and two estimates, say A and B and their standard errors. I need to find the standard error of a ratio A/B [or A/(1-B)].  I guess the main problem is that I don't know the correlation between A and B. I also guess that makes this unsolvable, since the correlation plays a major role. Is this correct?  If I knew it, what would be the way to calculate the standard error?
64709
I have conducted an experiment with multiple (categorical) conditions per subject, and multiple subject measurements.  My data-frame in short: A subject has one property, is_frisian which is either 0 or 1 depending on the subject. And it is tested for two conditions, person and condition. The measurement variable is error, which is either 0 or 1.  My mixed linear model in R is:  &gt; model &lt;- lmer(error~is_frisian*condition*person+(1|subject_id), data=output)   However, the residuals plot of this model gives an unexpected (?) result.    I was taught that this plot should show randomly scattered points, and they should be normal distributed. When plotting the density of the fitted and the residuals, it shows a reasonable normal distribution. The lines you can see in the graph, however, how is this to be explained? And is this okay?  The only thing I could come up with is that the graph has two lines due to the categorical variables. The output variable error is either 0 or 1. But I do not have that much knowledge of the underlying system to confirm this. And then again, the lines also seem to have a low negative slope, is this then perhaps a problem?  UPDATE:  &gt; model &lt;- glmer(error~is_frisian*condition*person + (1|subject_id), data=output, family='binomial') &gt; binnedplot(fitted(model),resid(model))   Gives the following result:    FINAL EDIT:  The density-plots have been omitted, they have nothing to do with satisfaction of assumptions in this case. For a list of assumptions on logistic regression (when using family=binomial),
78874
I am trying to run the following model:         I(week^2):mutation_status +         (week + I(week^2) | subject_id) ,          data = sim_dat)  This is the output I get from this model. The correlation between week and I(week)^2 is rather high (-0.95) and I was curious to know how I can change my data simulation code to lower that value. I am also having issues with convergence and would love to get thoughts on how to avoid issues regarding convergence. Thank you! Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerModLmerTest]  Formula: fetal_weight ~ week + mutation_status + week:mutation_status + I(week^2) +  I(week^2):mutation_status + (week + I(week^2) | subject_id) Data: sim_dat  REML criterion at convergence: 114  Scaled residuals:      Min      1Q  Median      3Q     Max  -1.9343 -0.3980 -0.0075  0.4168  1.9928   Random effects:  Groups     Name        Variance  Std.Dev. Corr         subject_id (Intercept)   1.51524  1.2309                          week         61.01832  7.8114   0.16                   I(week^2)   496.20236 22.2756  -0.19 -0.95  Residual                 0.05892  0.2427              Number of obs: 100, groups:  subject_id, 20  Fixed effects:                            Estimate Std. Error       df t value       Pr(&gt;|t|)     (Intercept)                  3.1331     0.4226  17.9961   7.413 0.000000713699 *** week                        -1.2465     3.5210  18.0053  -0.354          0.727     mutation_statusY             0.5061     0.5977  17.9961   0.847          0.408     I(week^2)                   34.5706    10.8146  18.0084   3.197          0.005 **  week:mutation_statusY       -1.0697     4.9795  18.0053  -0.215          0.832     mutation_statusY:I(week^2) 202.5521    15.2942  18.0084  13.244 0.000000000101 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Correlation of Fixed Effects:             (Intr) week   mttn_Y I(w^2) wk:m_Y week        -0.160                             muttn_sttsY -0.707  0.113                      I(week^2)    0.145 -0.966 -0.103               wk:mttn_stY  0.113 -0.707 -0.160  0.683        mtt_Y:I(^2) -0.103  0.683  0.145 -0.707 -0.966  Code for data simulation: set.seed(123)  # set the seed J &lt;- 20  # number of individuals (clusters) cs &lt;- 5  # number of time points (cluster size) gam &lt;- c(0, 0.75, 1.34)  # fixed effects G &lt;- matrix(c(1.75, 0, 0,               0, 0.0125, 0,               0,0,0.0625),              nrow = 3)  # random effect variances (G-matrix) sigma2 &lt;- 0.08  mutation_status &lt;- rep(0:1, each = 5) mutation_status &lt;- rep(mutation_status, times = 10)  X &lt;- cbind(1, seq_len(cs),             (seq_len(cs))^2)  # for each individual X &lt;- X[rep(seq_len(cs), J), ]  X &lt;- X[rep(seq_len(cs), J), ]  X[,2] &lt;- X[,2]*0.05 X[,3] &lt;- X[,3]*0.05  # repeat each row cs times pid &lt;- seq_len(J)  # individual id pid &lt;- rep(pid, each = cs)   # Generate person-level (lv-2) random effects uj &lt;- lmf::rmnorm(J, mean = rep(0, 2), varcov = G)  # Generate repeated-measure-level (lv-1) error term eij &lt;- rnorm(J * cs, sd = sqrt(sigma2))  # Compute beta_j's betaj &lt;- matrix(gam, nrow = J, ncol = 3, byrow = TRUE) + uj  # Compute outcome: y &lt;- rowSums(X * betaj[pid, ]) + eij + mutation_status*X[ , 3]*10 + 3  # Form a data frame sim_dat &lt;- tibble(y, time = X[ , 2], pid, mutation_status) sim_dat &lt;- sim_dat[,c(3,2,4,1)] colnames(sim_dat) &lt;- c(&quot;subject_id&quot;, &quot;week&quot;, &quot;mutation_status&quot;, &quot;fetal_weight&quot;) sim_dat$week_2 &lt;-sim_dat$week*5 + 15  sim_dat$mutation_status &lt;- as.factor(ifelse(sim_dat$mutation_status == 0, &quot;N&quot;, &quot;Y&quot;)) sim_dat$subject_id &lt;- as.factor(sim_dat$subject_id ) #Model Building
100272
I am a little bit confused on what the assumptions of linear regression are. So far I checked whether:  all of the explanatory variables correlated linearly with the response variable. (This was the case) there was any collinearity among the explanatory variables. (there was little collinearity). the Cook's distances of the datapoints of my model are below 1 (this is the case, all distances are below 0.4, so no influence points). the residuals are normally distributed. (this may not be the case)  But I then read the following:  violations of normality often arise either because (a) the distributions of the dependent and/or independent variables are themselves significantly non-normal, and/or (b) the linearity assumption is violated.  Question 1 This makes it sound as if the independent and depend variables need to be normally distributed, but as far as I know this is not the case. My dependent variable as well as one of my independent variables are not normally distributed. Should they be? Question 2 My QQnormal plot of the residuals look like this:   That slightly differs from a normal distribution and the shapiro.test also rejects the null hypothesis that the residuals are from a normal distribution: &gt; shapiro.test(residuals(lmresult)) W = 0.9171, p-value = 3.618e-06  The residuals vs fitted values look like:   What can I do if my residuals are not normally distributed? Does it mean the linear model is entirely useless?
1017
For my own understanding, I am interested in manually replicating the calculation of the standard errors of estimated coefficients as, for example, come with the output of the lm() function in R, but haven't been able to pin it down. What is the formula / implementation used?
20940
I used the inbuilt dataset stackloss here. I used following R code for creating the multiple regression model.  &gt; head(stackloss)   Air.Flow Water.Temp Acid.Conc. stack.loss   1        80         27         89         42   2        80         27         88         37   3        75         25         90         37   4        62         24         87         28   5        62         22         87         18   6        62         23         87         18  &gt; fit &lt;- lm(stack.loss ~ .,data=stackloss) &gt; plot(fit)   And when I plot the model, the four following graphs are displayed.      How to interpret these graphs? What are the significance of these graphs? How to use these graphs?
92298
I have just read a paper in which the authors carried out a multiple regression with two predictors. The overall r-squared value was 0.65. They provided a table which split the r-squared between the two predictors. The table looked like this:              rsquared beta    df pvalue whole model     0.65   NA  2, 9  0.008 predictor 1     0.38 1.01 1, 10  0.002 predictor 2     0.27 0.65 1, 10  0.030   In this model, ran in R using the mtcars dataset, the overall r-squared value is 0.76.  summary(lm(mpg ~ drat + wt, mtcars))  Call: lm(formula = mpg ~ drat + wt, data = mtcars)  Residuals:     Min      1Q  Median      3Q     Max  -5.4159 -2.0452  0.0136  1.7704  6.7466   Coefficients:             Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)   30.290      7.318   4.139 0.000274 *** drat           1.442      1.459   0.989 0.330854     wt            -4.783      0.797  -6.001 1.59e-06 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 3.047 on 29 degrees of freedom Multiple R-squared:  0.7609,    Adjusted R-squared:  0.7444  F-statistic: 46.14 on 2 and 29 DF,  p-value: 9.761e-10   How can I split the r-squared value between the two predictor variables?
41097
I am interested in getting an unbiased estimate of $R^2$ in a multiple linear regression.  On reflection, I can think of two different values that an unbiased estimate of $R^2$ might be trying to match.    Out of sample $R^2$: the r-square that would be obtained if the regression equation obtained from the sample (i.e., $\hat{\beta}$) were applied to an infinite amount of data external to the sample but from the same data generating process. Population $R^2$: The r-square that would be obtained if an infinite sample were obtained and the model fitted to that infinite sample (i.e., $\beta$) or alternatively just the R-square implied by the known data generating process.   I understand that  is designed to compensate for the overfitting observed in sample $R^2$. Nonetheless, it's not clear whether adjusted $R^2$ is actually an unbiased estimate of $R^2$, and if it is an unbiased estimate, which of the above two definitions of $R^2$ it is aiming to estimate.  Thus, my questions:   What is an unbiased estimate of what I call above out of sample $R^2$? What is an unbiased estimate of what I call above population $R^2$? Are there any references that provide simulation or other proof of the unbiasedness?