After a series of tests and exploratory analysis, I have design a linear model in R, the model is the following:  summary(gn)  Call: lm(formula = NA. ~ I(PC^0.25) + I(((PI)^2)), data = DSET)  Residuals:     Min      1Q  Median      3Q     Max  -425.22  -87.46   -2.30   79.11  396.14   Coefficients:               Estimate Std. Error t value Pr(&gt;|t|)     (Intercept) -2.047e+03  1.094e+02  -18.71  &lt; 2e-16 *** I(PC^0.25)   1.206e+03  4.231e+01   28.52  &lt; 2e-16 *** I(((PI)^2)) -5.242e-02  1.233e-02   -4.25 2.81e-05 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 127.2 on 319 degrees of freedom Multiple R-squared:  0.7475,    Adjusted R-squared:  0.746  F-statistic: 472.3 on 2 and 319 DF,  p-value: &lt; 2.2e-16   Then, for example, I calculate the interval confidence for a certain input:  &gt; a&lt;-data.frame(PC=58,PI=12) &gt; a&lt;-predict(gn,newdata=a,interval="prediction",level=0.95) &gt; a        fit      lwr      upr 1 1274.515 1022.976 1526.054   Finally, given the coefficients of the model and the same input, I try to calculate the same confidence interval in EXCEL. Using the formula for the 95% interval of each coefficient, &lt;estimate coefficient - 1.96*error, estimate coefficient + 1.96*error&gt;, the fitted value is reasonable but the extremes are far different.  So, the results are : fitted value=1274.499125 , lower=826.5811979 , upper=1722.417051. As you see, the extremes are different from the output of R.  What's my error? What should I do to get the same intervale I got in R?
27355
Let's take the following example:  set.seed(342) x1 &lt;- runif(100) x2 &lt;- runif(100) y &lt;- x1+x2 + 2*x1*x2 + rnorm(100) fit &lt;- lm(y~x1*x2)   This creates a model of y based on x1 and x2, using a OLS regression.  If we wish to predict y for a given x_vec we could simply use the formula we get from the summary(fit).  However, what if we want to predict the lower and upper predictions of y? (for a given confidence level).  How then would we build the formula?
107411
I am interested in comparing standardized regression coefficients within the same model. A linear hypothesis test (e.g., car::linearHypothesis in R) can tell me whether the difference between them is significant, and I can calculate the difference to get a point estimate. However, I am not sure how to get a CI around that. Currently, I believe that I can use $\beta_\delta$ + 1.96 * $\sqrt{SE1^2+SE2^2}$ in the same way that I would deal with the difference between two means, but I am not sure if that is applicable in this case?
37989
I have following simple X and Y vectors:   &gt; X [1] 1.000 0.063 0.031 0.012 0.005 0.000 &gt; Y [1] 1.000 1.000 1.000 0.961 0.884 0.000 &gt;  &gt; plot(X,Y)     I want to do regression using log of X. To avoid getting log(0), I try to put +1 or +0.1 or +0.00001 or +0.000000000000001 :   &gt; summary(lm(Y~log(X))) Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :    NA/NaN/Inf in 'x' &gt; summary(lm(Y~log(1+X)))  Call: lm(formula = Y ~ log(1 + X))  Residuals:        1        2        3        4        5        6  -0.03429  0.22189  0.23428  0.20282  0.12864 -0.75334   Coefficients:             Estimate Std. Error t value Pr(&gt;|t|)   (Intercept)   0.7533     0.1976   3.812   0.0189 * log(1 + X)    0.4053     0.6949   0.583   0.5910   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.4273 on 4 degrees of freedom Multiple R-squared:  0.07838,   Adjusted R-squared:  -0.152  F-statistic: 0.3402 on 1 and 4 DF,  p-value: 0.591  &gt; summary(lm(Y~log(0.1+X)))  Call: lm(formula = Y ~ log(0.1 + X))  Residuals:        1        2        3        4        5        6  -0.08099  0.20207  0.23447  0.21870  0.15126 -0.72550   Coefficients:              Estimate Std. Error t value Pr(&gt;|t|)   (Intercept)    1.0669     0.3941   2.707   0.0537 . log(0.1 + X)   0.1482     0.2030   0.730   0.5058   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.4182 on 4 degrees of freedom Multiple R-squared:  0.1176,    Adjusted R-squared:  -0.103  F-statistic: 0.5331 on 1 and 4 DF,  p-value: 0.5058  &gt; summary(lm(Y~log(0.00001+X)))  Call: lm(formula = Y ~ log(1e-05 + X))  Residuals:        1        2        3        4        5        6  -0.24072  0.02087  0.08796  0.13872  0.14445 -0.15128   Coefficients:                Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)     1.24072    0.12046  10.300 0.000501 *** log(1e-05 + X)  0.09463    0.02087   4.534 0.010547 *   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.1797 on 4 degrees of freedom Multiple R-squared:  0.8371,    Adjusted R-squared:  0.7964  F-statistic: 20.56 on 1 and 4 DF,  p-value: 0.01055  &gt;  &gt; summary(lm(Y~log(0.000000000000001+X)))  Call: lm(formula = Y ~ log(1e-15 + X))  Residuals:         1         2         3         4         5         6  -0.065506  0.019244  0.040983  0.031077 -0.019085 -0.006714   Coefficients:                Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)     1.06551    0.02202   48.38 1.09e-06 *** log(1e-15 + X)  0.03066    0.00152   20.17 3.57e-05 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.04392 on 4 degrees of freedom Multiple R-squared:  0.9903,    Adjusted R-squared:  0.9878  F-statistic: 406.9 on 1 and 4 DF,  p-value: 3.565e-05   The output is different in all cases. What is the correct value to put to avoid log(0) in regression? What is the correct method for such situations.  Edit: my main aim is to improve prediction of the regression model by adding log term, i.e.: lm(Y ~ X + log(X))
66581
I am a user more familiar with R, and have been trying to estimate random slopes (selection coefficients) for about 35 individuals over 5 years for four habitat variables. The response variable is whether a location was "used" (1) or "available" (0) habitat ("use" below).  I am using a Windows 64-bit computer.  In R version 3.1.0, I use the data and expression below. PS, TH, RS, and HW are fixed effects (standardized, measured distance to habitat types). lme4 V 1.1-7.   str(dat) 'data.frame':   359756 obs. of  7 variables:  $ use     : num  1 1 1 1 1 1 1 1 1 1 ...  $ Year    : Factor w/ 5 levels "1","2","3","4",..: 4 4 4 4 4 4 4 4 3 4 ...  $ ID      : num  306 306 306 306 306 306 306 306 162 306 ...  $ PS: num  -0.32 -0.317 -0.317 -0.318 -0.317 ...  $ TH: num  -0.211 -0.211 -0.211 -0.213 -0.22 ...  $ RS: num  -0.337 -0.337 -0.337 -0.337 -0.337 ...  $ HW: num  -0.0258 -0.19 -0.19 -0.19 -0.4561 ...  glmer(use ~  PS + TH + RS + HW +      (1 + PS + TH + RS + HW |ID/Year),      family = binomial, data = dat, control=glmerControl(optimizer="bobyqa"))   glmer gives me parameter estimates for the fixed effects that make sense to me, and the random slopes (which I interpret as selection coefficients to each habitat type) also make sense when I investigate the data qualitatively. The log-likelihood for the model is -3050.8.  However, most research in animal ecology do not use R because with animal location data, spatial autocorrelation can make the standard errors prone to type I error. While R uses model-based standard errors, empirical (also Huber-white or sandwich) standard errors are preferred.   While R does not currently offer this option (to my knowledge - PLEASE, correct me if I am wrong), SAS does - although I do not have access to SAS, a colleague agreed to let me borrow his computer to determine if the standard errors change significantly when the empirical method is used.  First, we wished to ensure that when using model-based standard errors, SAS would produce similar estimates to R - to be certain that the model is specified the same way in both programs. I don't care if they are exactly the same - just similar. I tried (SAS V 9.2):  proc glimmix data=dat method=laplace;    class year id;    model use =  PS TH RS HW / dist=bin solution ddfm=betwithin;    random intercept PS TH RS HW / subject = year(id) solution type=UN; run;title;   I also tried various other forms, such as adding lines  random intercept / subject = year(id) solution type=UN; random intercept PS TH RS HW / subject = id solution type=UN;   I tried without specifying the   solution type = UN,   or commenting out  ddfm=betwithin;   No matter how we specify the model (and we have tried many ways), I cannot get the random slopes in SAS to remotely resemble those output from R - even though the fixed effects are similar enough. And when I mean different, I mean that not even the signs are the same. The -2 Log Likelihood in SAS was 71344.94.   I can't upload my full dataset; so I made a toy dataset with only the records from three individuals. SAS gives me output in a few minutes; in R it takes over an hour. Weird. With this toy dataset I'm now getting different estimates for the fixed effects.   My question: Can anyone shed light on why the random slopes estimates might be so different between R and SAS? Is there anything I can do in R, or SAS, to modify my code so that the calls produce similar results? I'd rather change the code in SAS, since I "believe" my R estimates more.   I'm really concerned with these differences and want to get to the bottom of this problem!  My output from a toy dataset that uses only three of the 35 individuals in the full dataset for R and SAS are included as jpegs.         EDIT AND UPDATE:  As @JakeWestfall helped discover, the slopes in SAS do not include the fixed effects. When I add the fixed effects, here is the result - comparing R slopes to SAS slopes for one fixed effect, "PS", between programs: (Selection coefficient = random slope). Note the increased variation in SAS.
23704
I always use lm() in R to perform linear regression of $y$ on $x$. That function returns a coefficient $\beta$ such that $$y = \beta x.$$  Today I learned about total least squares and that princomp() function (principal component analysis, PCA) can be used to perform it. It should be good for me (more accurate). I have done some tests using princomp(), like:  r &lt;- princomp( ~ x + y)   My problem is: how to interpret its results? How can I get the regression coefficient? By "coefficient" I mean the number $\beta$ that I have to use to multiply the $x$ value to give a number close to $y$.
23544
I have a simple linear OLS regression $Y_i = \alpha+ \beta_1 X_{1i} + \beta_2 X_{2i} + e_i$ where $e_i \sim N(0,\sigma)$. I have estimated the regression from the data and obtained estimates for my coefficients as well as the corresponding covariance matrix. Assume my dataset is large (>500)  Now I would like to construct the 95% confidence interval for $Z = a_1 \beta_1 + a_2 \beta_2$. My first thought was that I could assume that Z follows a normal distribution because my regression coefficients are approximately normal and a linear combination of normal variables is also a normal variable. I do find some sources where they treat $Z$ as a normal random variable.  But now I am in doubt because I found that a linear combination of normal random variables is only normal if they are independent. So is it still save to assume $Z$ is random?
16711
The help pages in R assume I know what those numbers mean, but I don't. I'm trying to really intuitively understand every number here. I will just post the output and comment on what I found out. There might (will) be mistakes, as I'll just write what I assume. Mainly I'd like to know what the t-value in the coefficients mean, and why they print the residual standard error.   Call: lm(formula = iris$Sepal.Width ~ iris$Petal.Width)  Residuals:      Min       1Q   Median       3Q      Max  -1.09907 -0.23626 -0.01064  0.23345  1.17532    This is a 5-point-summary of the residuals (their mean is always 0, right?). The numbers can be used (I'm guessing here) to quickly see if there are any big outliers. Also you can already see it here if the residuals are far from normally distributed (they should be normally distributed).  Coefficients:                  Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)       3.30843    0.06210  53.278  &lt; 2e-16 *** iris$Petal.Width -0.20936    0.04374  -4.786 4.07e-06 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1    Estimates $\hat{\beta_i}$, computed by least squares regression. Also, the standard error is $\sigma_{\beta_i}$. I'd like to know how this is calculated.  I have no idea where the t-value and the corresponding p-value come from. I know $\hat{\beta}$ should be normal distributed, but how is the t-value calculated?  Residual standard error: 0.407 on 148 degrees of freedom   $\sqrt{ \frac{1}{n-p} \epsilon^T\epsilon }$, I guess. But why do we calculate that, and what does it tell us?  Multiple R-squared: 0.134,  Adjusted R-squared: 0.1282    $ R^2 = \frac{s_\hat{y}^2}{s_y^2} $, which is $ \frac{\sum_{i=1}^n (\hat{y_i}-\bar{y})^2}{\sum_{i=1}^n (y_i-\bar{y})^2} $. The ratio is close to 1 if the points lie on a straight line, and 0 if they are random.  What is the adjusted R-squared?  F-statistic: 22.91 on 1 and 148 DF,  p-value: 4.073e-06    F and p for the whole model, not only for single $\beta_i$s as previous. The F value is $ \frac{s^2_{\hat{y}}}{\sum\epsilon_i} $. The bigger it grows, the more unlikely it is that the $\beta$'s do not have any effect at all.
5378
I originally posted this in stackoverflow (as given ) but was told to try here since it might be more relevant here. I am very new to statistics and R in general so my question might be a bit dumb, but since I cannot find my solutions online I thought I should try ask it here.  I have a data frame dataset of a whole lot of different variables very similar to as follows:   Item | Size   | Value    | Town ---------------------------------- A     |  10    |   800   | 1 B     |  11    |   100   | 2 A     |  17    |   900   | 2 D     |  13    |   200   | 3 B     |  15    |   500   | 1 C     |  12    |   250   | 3 E     |  14    |    NA   | 2 A     |        |   800   | 1 C     |        |   800   | 2   Basically, I have to try and 'guess' the Size based on the type of Item, it's Value, and the Town it was sold in, so I think a regression method would be a good idea.  I now use a regression as follows:  lm(Size ~ factor(Item) + factor(Town) + Value,...)  The problem however, occurs when I try and predict the Size using this model. I have the following lines of code:  pmodel &lt;- lm(Size ~ factor(Item) + factor(Town) + Value,...) prediction &lt;- predict(pmodel, dataset2)   (where dataset2 is a subset of dataset which has all the empty "Size" values which I want to predict)  But this now comes up with the error:  Error in model.frame.default(Terms, newdata, na.action=na.action, xlev = object$xlevels): factor factor(Town) has new levels   Is there any way around this, or is there any way I can get the model to make prediction based on the other values if it encounters a 'level' which is not in the original dataset?
78243
The regression model y= b0+ b1 x + b2 x^2 + b3 x^3  and the second regression model y = b0 +b1 (x-u) + b2 (x-u)^2 + b3 (x-u)^3 where u is the mean of x  These two models lead to the same curves, or says the same fitted values. I understand centering the variables for polynomial or interaction regressors will lower the correlations of regressors and therefore reduce the multicolinearity, variance inflation, prediction interval... But they are still the same model, does the difference between the variances of two identical models really mean something important?  R code example:  y &lt;- rnorm(100) x &lt;- runif(100) u &lt;- mean(x)  model1&lt;- lm(y~x+I(x^2)+I(x^3)) model2 &lt;- lm(y~x+I((x-u)^2)+I((x-u)^3))  summary(model1) summary(model2)   They have the same sd error, R squared.... Their fitted values are the same as well.
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?
292
Imagine that you repeat an experiment three times. In each experiment, you collect triplicate measurements. The triplicates tend to be fairly close together, compared to the differences among the three experimental means. Computing the grand mean is pretty easy. But how can one compute a confidence interval for the grand mean?  Sample data:  Experiment 1:   34, 41, 39  Experiment 2:   45, 51, 52  Experiment 3:   29, 31, 35  Assume that the replicate values within an experiment follow a Gaussian distribution, as does the mean values of each experiment. The SD of variation within an experiment is smaller than the SD among experimental means. Assume also that there is no ordering of the three values in each experiment. The left-to-right order of the three values in each row is entirely arbitrary.  The simple approach is to first compute the mean of each experiment: 38.0, 49.3, and 31.7, and then compute the mean, and its 95% confidence interval, of those three values. Using this method, the grand mean is 39.7 with the 95% confidence interval ranging from 17.4 to 61.9.   The problem with that approach is that it totally ignores the variation among triplicates. I wonder if there isn't a good way to account for that variation.
6284
I am asked to determine using an appropriate test whether the internal concentration of isoleucine is statistically different between appropriate experiments. So I constructed a data frame in R as follow: dat &lt;- data.frame(experiment = c('no glu', 'glu', 'gluFCCP', 'gluV', 'gluN', 'gluNV'),internal_concentration = c(2.38, 8.57, 3.42, 6.17, 4.58,3.51))  print(dat) experiment     internal_concentration 1     no glu                   2.38 2        glu                   8.57 3    gluFCCP                   3.42 4       gluV                   6.17 5       gluN                   4.58 6      gluNV                   3.51  I chose to perform an ANOVA test: model&lt;-lm(dat$internal_concentration~dat$experiment)  anova(model)  Analysis of Variance Table  Response: dat$internal_concentration                Df Sum Sq Mean Sq F value Pr(&gt;F) dat$experiment  5 25.558  5.1117                Residuals       0  0.000    I need the F and p values to determine if there's a significant difference, but this result is not working. How should I interpret this error and any suggestion on how to change the code/the test used? I just have one value for the concentration for each experiment.
41895
I try to reproduce with optim the results from a simple linear regression fitted with glm or even nls R functions. The parameters estimates are the same but the residual variance estimate and the standard errors of the other parameters are not the same particularly when the sample size is low. I suppose that this is due differences in the way the residual standard error is calculated between Maximum Likelihood and Least square approaches (dividing by n or by n-k+1 see bellow in the example). I understand from my readings on the web that optimization is not a simple task but I was wondering if it would be possible to reproduce in a simple way the standard error estimates from glm while using optim.   Simulate a small dataset  set.seed(1) n = 4 # very small sample size ! b0 &lt;- 5 b1 &lt;- 2 sigma &lt;- 5 x &lt;- runif(n, 1, 100) y =  b0 + b1*x + rnorm(n, 0, sigma)    Estimate with optim  negLL &lt;- function(beta, y, x) {     b0 &lt;- beta[1]     b1 &lt;- beta[2]     sigma &lt;- beta[3]     yhat &lt;- b0 + b1*x     likelihood &lt;- dnorm(y, yhat, sigma)     return(-sum(log(likelihood))) }  res &lt;- optim(starting.values, negLL, y = y, x = x, hessian=TRUE) estimates &lt;- res$par     # Parameters estimates se &lt;- sqrt(diag(solve(res$hessian))) # Standard errors of the estimates cbind(estimates,se)       &gt; cbind(estimates,se)       estimates         se b0     9.016513 5.70999880 b1     1.931119 0.09731153 sigma  4.717216 1.66753138   Comparison with glm and nls  &gt; m &lt;- glm(y ~ x) &gt; summary(m)$coefficients             Estimate Std. Error   t value    Pr(&gt;|t|) (Intercept) 9.016113  8.0759837  1.116411 0.380380963 x           1.931130  0.1376334 14.030973 0.005041162 &gt; sqrt(summary(m)$dispersion) # residuals standard error [1] 6.671833 &gt;  &gt; summary(nls( y ~ b0 + b1*x, start=list(b0 = 5, b1= 2)))  Formula: y ~ b0 + b1 * x  Parameters:    Estimate Std. Error t value Pr(&gt;|t|)    b0   9.0161     8.0760   1.116  0.38038    b1   1.9311     0.1376  14.031  0.00504 ** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 6.672 on 2 degrees of freedom   I can reproduce the different residual standard error estimates like this :   &gt; # optim / Maximum Likelihood estimate &gt; sqrt(sum(resid(m)^2)/n) [1] 4.717698 &gt;  &gt; # Least squares estimate (glm and nls estimates) &gt; k &lt;- 3 # number of parameters &gt; sqrt(sum(resid(m)^2)/(n-k+1)) [1] 6.671833
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?
3446
I have noticed that the confidence interval for predicted values in an linear regression tends to be narrow around the mean of the predictor and fat around the minimum and maximum values of the predictor. This can be seen in plots of these 4 linear regressions:    I initially thought this was because most values of the predictors were concentrated around the mean of the predictor. However, I then noticed that the narrow middle of the confidence interval would occur even if many values of were concentrated around the extremes of the predictor, as in the bottom left linear regression, which lots of values of the predictor are concentrated around the minimum of the predictor.  is anyone able to explain why confidence intervals for the predicted values in an linear regression tend to be narrow in the middle and fat at the extremes?
119400
I'm using the standard formula for the confidence interval for standard deviation:  SQRT(((N-1)*s^2)/(CHISQ.INV((1-Confidence)/2, N-1)))   By mistake I entered a very low confidence level and I noticed that the resulting interval did not include the calculated standard deviation. This happens when the confidence is sufficiently low (e.g. 10%).  Can someone explain me the reason? I've double checked the results with Minitab and they are correct.
95631
When using factors in a linear model,  I would like to retrieve the 'average' effect as baseline (intercept), rather than level 1 of all factors. Is this possible?  tab &lt;- data.frame(gl(2,5),rnorm(10)+as.numeric(gl(2,5))) colnames(tab) &lt;- c('pred','dat') mylm &lt;- lm(dat~pred,data=tab) coefficients(mylm)   returns something like  (Intercept)       pred2    0.5879712   1.8542153    coefficients are on average 1 for both intercept and pred2. I would like to rather have intercept = 1.5, pred1 = -0.5 and pred2 = +0.5  model.matrix(mylm)   returns           (Intercept) pred2     1            1     0     2            1     0     3            1     0     4            1     0     5            1     0     6            1     1     7            1     1     8            1     1     9            1     1     10           1     1     attr(,"assign")     [1] 0 1     attr(,"contrasts")     attr(,"contrasts")$pred     [1] "contr.treatment"   I would like a matrix like this            (Intercept)  pred1   pred2     1            1    1     0     2            1    1     0     3            1    1     0     4            1    1     0     5            1    1     0     6            1    0     1     7            1    0     1     8            1    0     1     9            1    0     1     10           1    0     1   Does that make sense? if not, why? Is there some command to generate such a design matrix? I've been looking for a way to do this without success...