I want to fit a multiple regression to test the relationship between several independent variables and a dependent variable. A few independent variables are standard demographic variables, the rest of the IVs are summed scores of responses to likert items. The dv is also a sum of responses to likert items.   The books I have read before and during my analysis are Discovering statistics using SPSS, by Field, 4th ed., and R in action, data analysis and graphics with r.    The problem is that the homoscedasticity assumption has been violated. To test homoscedasticity, I used the ncvTest() available in Car package of R, and the Breusch and Pagan test available in olsrr package for R. Both  show violation of homoscedasticity.   Both books I have mentioned above, suggest that bootstrapping for robust standard errors and confidence intervals should be used when the homoscedasticity assumption is not met. In his book, Field suggests that when there is  heteroscedasticity, the tests of hypotheses may not be valid.   I have bootstraped using the Boot function in R package, as well as  using SPSS. My problem is that I  am unable to work out how to present my results after a  bootstrap, and which values may or should change other than standard errors as a result of the bootstrap.   I’ll use the output of both R and SPSS as examples to elaborate my question. In R,  I know that rlm(), Boot(), summary(ofboot) will give me the bootstrap standard errors. But, when I am going to report my results, should I report the r-squared, f, and t provided by rlm(), should bootstrapping influence r-squared, f and t values? In other words, when I'm constructing the regression results table, I will take the values of standard error from the  results of bootstrap, but can I take the t, f, and coefficients values returned by rlm()? SPSS provides two tables of coefficients, both with and without the bootstrap; the bootstrap coefficient table does not have a column for t, so, again, when reporting the results, do I report the t in coefficient table without the bootstrap, f in the ANOVA and r-squared in the model summary table. Field in his book, has provided a bootstrap example, but the process is shown graphically, and since I’m totally blind, I cannot see the image to determine whether his book has an answer to my dumbest question.
121836
If in standard OLS regressions two assumptions are violated (normal distribution of errors, homoscedasticity), is bootstrapping standard errors and confidence intervals an appropriate alternative to arrive at meaningful results with respect to the significance of regressor coefficients?  Do significance tests with bootstrapped standard errors and confidence intervals still "work" with heteroscedasticity?  If yes, what would be applicable confidence intervals that can be used in this scenario (percentile, BC, BCA)?  Finally, if bootstrapping is appropriate in this scenario, what would be the relevant literature that needs to be read and cited to arrive at this conclusion? Any hint would be greatly appreciated!
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
75641
I'm really confused about the difference in meaning regarding the context of linear regression of the following terms:   F statistic R squared Residual standard error   I found  which gave me great insight in the different terms involved in linear regression, however the terms mentioned above look a like quite a lot (as far as I understand). I will cite what I read and what confused me:      Residual Standard Error is measure of the quality of a linear   regression fit.......The Residual Standard Error is the average amount that the response (dist) will deviate from the true regression line.   1. This is thus actually the average distance of the observed values from the lm line?     The R-squared statistic provides a measure of how well the   model is fitting the actual data.   2. Now I'm getting confused because if RSE tells us how far our observed points deviate from the regression line a low RSE is actually telling us "your model is fitting well based on the observed data points" --> thus how good our models fits, so what is the difference between R squared and RSE?      F-statistic is a good indicator of whether there is a relationship   between our predictor and the response variables.   3. Is it true that we can have a F value indicating a strong relationship that is NON LINEAR so that our RSE is high and our R squared is low
122058
I want to predict future returns over a 20 days horizon using an ARMA-GARCH model fitted to my data. The goal is to estimate different risk measures like VaR or CVar. In particular say I use AR(1) GARCH(1,1). The smaple I use for estimation has 500 observation of daily logreturns. That's what I usually do:   Estimate AR, ARCH and GARCH coefficients Calculate standardized residual by dividing residuals by estimated conditional variances The standardized residuals constitute my INVARIANTS that is the i.i.d. series from which I extract bootstrap samples to generate scenarios.   The bootstrap samples are extracted by simulating a uniform between 1 and the sample size (500 in this case) and then taking the value corresponding to that position in the vector of the standardized residuals.  The problem is that I have only 500 standardized residuals and I think 500 is the maximum size of bootstrap samples I can extract.   My colleague instead extracts 100000 observations out of the original sample of $N=500$ observations.  I feel that this is somehow incorrect conceptually. Simulating only one step forward would produce exactly the same scenario as the initial one, but with repeted values that add no information.  My colleague claims that if he wants to project over a longer period, e.g. 20 days horizon, the 100000 extractions from the original sample of N=500 obs. would produce many different scenarios at the final horizon, providing a CDF that is smooth.  Actually thit is true because, although the values are simply repeated in the first step, after that they can sum up in many different ways.  That being said I don't feel this is right. I proposed an alternative that is:   From the standardized residuals create a smoothed empirical CDF, say kernel Exctract uniforms between 0 and 1 and feed it to the empirical smoothed CDF, i.e. inverse transform.   This way I feel more confortable to say that I can generate a bootstrab sample of size greater than the original one, but still I am not sure.  I am studying bootstrap theory on a book from Efron      Efron, Tibshirani - An Introduction to the Bootstrap - Springer US (1993)   but there are many concepts that I  don't understand yet.  My question are:   Would you give your opinion on the problem I just showed below? Would you suggest any valid matherial for studying bootstrap other than the book I mentioned? I think that this application of bootstrap is somehow different to the one explained in Efron's bookm, that is evaluating confidence intervals for estimated paramenters. What do you think about it?   Any comment would be much appreciated  I apologize for the length of the post but I tried to be as much concise as I could. Thank you
25743
I performed with R an ANOVA and I got significant differences. However when checking which pairs were significantly different using the Tukey's procedure I did not get any of them. How can this be possible?  Here is the code:  fit5_snow&lt;- lm(Response ~ Stimulus, data=audio_snow) anova(fit5_snow)  &gt; anova(fit5_snow) Analysis of Variance Table  Response: Response           Df Sum Sq Mean Sq F value  Pr(&gt;F)   Stimulus   5  73.79 14.7578  2.6308 0.02929 * Residuals 84 471.20  5.6095                   --- Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1   df&lt;-df.residual(fit5_snow) MSerror&lt;-deviance(fit5_snow)/df  comparison &lt;-  HSD.test(audio_snow$Response, audio_snow$Stimulus, df, MSerror, group=FALSE)  &gt; comparison &lt;-  HSD.test(audio_snow$Response, audio_snow$Stimulus, df, MSerror, group=FALSE)  Study:  HSD Test for audio_snow$Response   Mean Square Error:  5.609524   audio_snow$Stimulus,  means                  audio_snow.Response   std.err replication snow_dry_leaves            4.933333 0.6208034          15 snow_gravel                6.866667 0.5679258          15 snow_metal                 6.333333 0.5662463          15 snow_sand                  6.733333 0.5114561          15 snow_snow                  7.333333 0.5989409          15 snow_wood                  5.066667 0.7713110          15  alpha: 0.05 ; Df Error: 84  Critical Value of Studentized Range: 4.124617   Comparison between treatments means                                Difference   pvalue sig        LCL      UCL snow_gravel - snow_dry_leaves  1.9333333 0.232848     -0.5889913 4.455658 snow_metal - snow_dry_leaves   1.4000000 0.588616     -1.1223246 3.922325 snow_sand - snow_dry_leaves    1.8000000 0.307012     -0.7223246 4.322325 snow_snow - snow_dry_leaves    2.4000000 0.071587   . -0.1223246 4.922325 snow_wood - snow_dry_leaves    0.1333333 0.999987     -2.3889913 2.655658 snow_gravel - snow_metal       0.5333333 0.989528     -1.9889913 3.055658 snow_gravel - snow_sand        0.1333333 0.999987     -2.3889913 2.655658 snow_snow - snow_gravel        0.4666667 0.994348     -2.0556579 2.988991 snow_gravel - snow_wood        1.8000000 0.307012     -0.7223246 4.322325 snow_sand - snow_metal         0.4000000 0.997266     -2.1223246 2.922325 snow_snow - snow_metal         1.0000000 0.855987     -1.5223246 3.522325 snow_metal - snow_wood         1.2666667 0.687424     -1.2556579 3.788991 snow_snow - snow_sand          0.6000000 0.982179     -1.9223246 3.122325 snow_sand - snow_wood          1.6666667 0.393171     -0.8556579 4.188991 snow_snow - snow_wood          2.2666667 0.103505     -0.2556579 4.788991
102491
Can the group recommend a good introduction text/resource to applied resampling techniques? Specifically, I am interested in alternatives to classical parametric tests (e.g. t tests, ANOVA, ANCOVA) for comparing groups when assumptions such as normality are clearly violated.   An example problem type I would like to educate myself as to a better way to solve may involve something such as:  I) 2 Groups: Treatment and Control    Dependent Var: Change in account balance dollars after intervention  Covariate: Pre intervention account balance dollars.  Issue with applying ANCOVA: Many subjects will not have any change (many zeros).   II) 2 Groups: Treatment and Control  Dependent Var: new accounts added  Covariate: Pre intervention number of accounts.  *Many subjects will not have any added account (many zeros).   Can I use a bootstrap? A permutation test? This is the type of analysis I would like to apply nonparametric resampling methods to.
49931
I was reading about regression. I understand one of the how the variance of each $x_i$ needs to be normal distributed, but why doesn't that related to the $s_b$? It seems like it should. I was reading the AP College Board paper about Inference.      Two more assumptions are that the errors are normally distributed and that they have   the same standard deviation for all x values. Although it is not at all obvious, this is   not the same as saying that the residuals are normally distributed and all have the same   standard deviation. The difference is that the “errors” are the deviations between the   response variables and the “actual” underlying (invisible) linear phenomenon, whereas   the “residuals” are the deviations between the response variables and the least-squares   regression line, which is only an estimate of the “actual” line.   I understand (I think) that just because the error of the regression is normal doesn't mean that the variance of $x_i$ is normal because the sample could just randomly follow normality.  (don't know if I said that statistically correct), but if all the $x_i$ are normal and the variance are equal, then why doesn't that variance equal the variance that one would get for the $\beta$  I have gotten some great answers here that are very clear and I am hoping that I could get another answer to help clear up my confusion.    The kernel of my question is What exactly is the relationship between the variance of the residuals and the $x_i$ and is there any connection to the error, variance of $\beta$
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.
64658
I am using the epitools in R for calculating the confidence interval of relative risk.    There are three methods inside for calculations: namely Wald, Small and Boot.  I want to find some article describing the three methods, but I can't find any, can anyone help? Thanks!
33700
There are several threads on this site discussing  to determine if the OLS residuals are  normally distributed. Another way to evaluate the normality of the residuals with R code is provided in this excellent . This is another  on the practical difference between standardized and observed residuals.  But let's say the residuals are definitely not normally distributed, like in this . Here we have several thousand observations and clearly we must reject the normally-distributed-residuals assumption. One way to address the problem is to employ some form of robust estimator as explained in the answer. However I am not limited to OLS and in facts I would like to understand the benefits of other glm or non-linear methodologies.  What is the most efficient way to model data violating the OLS normality of residuals assumption? Or at least what should be the first step to develop a sound regression analysis methodology?
4346
Following on from : Imagine that you want to test for differences in central tendency between two groups (e.g., males and females) on a 5-point Likert item (e.g., satisfaction with life: Dissatisfied to Satisfied). I think a t-test would be sufficiently accurate for most purposes,  but that a bootstrap test of differences between group means would often provide more accurate estimate of confidence intervals. What statistical test would you use?
106557
I have searched far and wide for a non-parametric test for the 2x2 between subjects MANOVA I would like to do. So far I have found non-parametric equivalents for factorial ANOVA with corresponding R packages (). Does anyone have any ideas on how to do the same analysis, but with multiple DVs?
8582
I have a question on how a statistician would normally interpret an anova output. Say I have anova output from R.  &gt; summary(fitted_data)  Call: lm(formula = V1 ~ V2)  Residuals:      Min       1Q   Median       3Q      Max  -2.74004 -0.33827  0.04062  0.44064  1.22737   Coefficients:             Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)  2.11405    0.32089   6.588  1.3e-09 *** V2           0.03883    0.01277   3.040  0.00292 **  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 0.6231 on 118 degrees of freedom Multiple R-squared: 0.07262,    Adjusted R-squared: 0.06476  F-statistic:  9.24 on 1 and 118 DF,  p-value: 0.002917   &gt; anova(fit) Analysis of Variance Table  Response: V1            Df Sum Sq Mean Sq F value   Pr(&gt;F)    V2          1  3.588  3.5878  9.2402 0.002917 ** Residuals 118 45.818  0.3883                     --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1    From the above, I guess the most important value is Pr(>F), right? So this Pr, is less than 0.05 (95% level). How should my "explain" this? Do I explain it in "association", ie, V2 and V1 are associated (or not) ? or in terms of "significance"? I always felt that I couldn't understand when people say "This value is significant....". So what is "significant"? Is there a more intuitive form of explanation? like "I am 95% confident that ...." .   Also, is the Pr value the only important piece of information? or can i also look at residuals and the rest of the output to "explain" the result? thanks
104835
I have tried to reproduce some research (using PCA) from SPSS in R. In my experience, principal()  from package psych was the only function that came close (or if my memory serves me right, dead on) to match the output. To match the same results as in SPSS, I had to use parameter principal(..., rotate = "varimax"). I have seen papers talk about how they did PCA, but based on the output of SPSS and use of rotation, it sounds more like Factor analysis.  Question: Is PCA, even after rotation (using ), still PCA? I was under the impression that this might be in fact Factor analysis... In case it's not, what details am I missing?
5647
Let's have some linear model, for example just simple ANOVA:  # data generation set.seed(1.234)                       Ng &lt;- c(41, 37, 42)                     data &lt;- rnorm(sum(Ng), mean = rep(c(-1, 0, 1), Ng), sd = 1)       fact &lt;- as.factor(rep(LETTERS[1:3], Ng))   m1 = lm(data ~ 0 + fact) summary(m1)   Result is as follows:  Call: lm(formula = data ~ 0 + fact)  Residuals:      Min       1Q   Median       3Q      Max  -2.30047 -0.60414 -0.04078  0.54316  2.25323   Coefficients:       Estimate Std. Error t value Pr(&gt;|t|)     factA  -0.9142     0.1388  -6.588 1.34e-09 *** factB   0.1484     0.1461   1.016    0.312     factC   1.0990     0.1371   8.015 9.25e-13 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 0.8886 on 117 degrees of freedom Multiple R-squared: 0.4816,     Adjusted R-squared: 0.4683  F-statistic: 36.23 on 3 and 117 DF,  p-value: &lt; 2.2e-16    Now I try two different methods to estimate confidence interval of these parameters  c = coef(summary(m1))  # 1st method: CI limits from SE, assuming normal distribution cbind(low = c[,1] - qnorm(p = 0.975) * c[,2],      high = c[,1] + qnorm(p = 0.975) * c[,2])  # 2nd method confint(m1)   Questions:   What is the distribution of estimated linear regression coefficients? Normal or $t$? Why do both methods yield different results? Assuming normal distribution and correct SE, I'd expect both methods to have the same result.   Thank you very much!  data ~ 0 + fact  EDIT after an answer:  The answer is exact, this will give exactly the same result as confint(m1)!  # 3rd method cbind(low = c[,1] - qt(p = 0.975, df = sum(Ng) - 3) * c[,2],      high = c[,1] + qt(p = 0.975, df = sum(Ng) - 3) * c[,2])
49169
I am referring to the question and its answers:  by @Clark Chong and answers/comments by @Frank Harrell. and to the question  and the comments.    I have read the paper D.W. Hosmer, T. Hosmer, S. Le Cessie, S. Lemeshow, "A comparison of goodness-of-fit tests for the Logistic regression model",  Statistics in Medicine, Vol. 16, 965-980 (1997).    After reading I was confused because the question I referred to asks explicitly for "(probability) predictive ability", which is in my opinion not the same as what the goodness-of-fit-tests in the paper supra aim at:  As most of us know, the logistic regression assumes an S-shaped link between the explanatory variables and the probability of success, the functional form for the S-shape is   $P(y=1|_{x_i})=\frac{1}{1+e^{-(\beta_0+\sum_i \beta_i x_i)}}$  Without pretending that there are no shortcomings with the Hosmer&ndash;Lemeshow test, I think that we have to distinguish between tests for the (a) '(probability) predictive ability' and (b) 'goodness-of-fit'.    The former's goal is to test whether the probabilities are well predicted, while the goodness-of-fit tests test whether the S-shaped function above is the 'right' function.  More formally:    tests for 'probability predictive ability tests' have a $H_0$ stating that the success probabilities are well predicted by the model; while for goodness-of-fit tests $H_0$ is (see Hosmer et. al. ) that the S-shaped functional form supra is the correct one.  Hosmer et al. perform simulations where they find the power to detect two types of deviations from the null namely that the link function is wrong or that the exponent in the denominator is not linear.    Obviously, if the above function has the 'right' functional form (so if the tests concludes that we can accept $H_0$ for the goodness-of-fit test), then the predicted probabilities will be fine, ...  First remark  ...however, accepting the $H_0$ is a weak conclusion as explained in .   First question  The most important question/remark that I have is that if the goodness-of-fit $H_0$ is rejected, then the conclusion of the test is that the functional form was not the 'right' one, however, does this imply that the probabilities are not well predicted ?  Second question  Furthermore, I want to point to the conclusions of Hosmer et. al; (I cite from the abstract):   ''An examination of the performance of the tests when the correct model has a quadratic term but a model containing only the linear term has been fit shows that the Pearson chi-square, the unweighted sum-of-squares,the Hosmer-Lemeshow decile of risk, the smoothed residual sum-of-squares and Stukel’s score test,have power exceeding 50 per cent to detect moderate departures from linearity when the sample size is 100 and have power over 90 per cent for these same alternatives for samples of size 500. All tests had no power when the correct model had an interaction between a dichotomous and continuous covariate but only the continuous covariate model was fit. Power to detect an incorrectly specifed link was poor for samples of size 100. For samples of size 500 Stukel's score test had the best power but it only exceeded 50 per cent to detect an asymmetric link function. The power of the unweighted sum-of-squares test to detect an incorrectly specifed link function was slightly less than Stukel’s score test ''  Can I conclude from this which test has more power or that Hosmer&ndash;Lemeshow has less power (to detect these specific anomalies) ?  Second remark  The paper by Hosmer et. al. that I referred to supra, compute (simulate) the power to detect specific anomalies (the power can only be computed if an $H_1$ is specified).  This does at my opinion not imply that these results can be generalised to ''all possible alternatives $H_1$'' ?
83552
Suppose I have a statistics (say a price index) and I want to obtain standard errors for it.  I have heard that blind application of bootstrap may not be a good practice. If true  1- What could go wrong if I just apply nonparametric bootstrap and obtain standard errors for my index?   2- If it is OK, are the rsulting standard errors robust to heteroskedasticity, serial correlation etc.?  Thank you