I am performing deviance goodness-of-fit test on my model (used negative binomial regression), and the R summary() table of my model gives the following:  Call: glm.nb(formula = topPagesCount ~ DB_LEGAL_STATUS_CODE_V2 + BCORP_INDUSTRY_DENSITY +  DENSE_MSA + DENSE_MSA * BCORP_INDUSTRY_DENSITY, data = greatDF,  control = glm.control(maxit = 100), init.theta = 0.2103980872,  link = log)  Deviance Residuals:  Min       1Q   Median       3Q      Max   -1.9154  -1.2815  -0.7786  -0.2994   2.3033    Coefficients:                                     Estimate Std. Error z value Pr(&gt;|z|)       ( Intercept)                      4.4117492  0.4566746   9.661  &lt; 2e-16 ***   DB_LEGAL_STATUS_CODE_V212        -0.6420125  0.3006310  -2.136  0.03272 *     BCORP_INDUSTRY_DENSITY            0.0067614  0.0092156   0.734  0.46314       DENSE_MSA                         0.0311327  0.0104588   2.977  0.00291 **    BCORP_INDUSTRY_DENSITY:DENSE_MSA -0.0006298  0.0002075  -3.035  0.00241 **     ---   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1    (Dispersion parameter for Negative Binomial(0.2104) family taken to be 1)   Null deviance: 393.12  on 304  degrees of freedom  Residual deviance: 362.31  on 300  degrees of freedom  AIC: 2794.2   Number of Fisher Scoring iterations: 1             Theta:  0.2104        Std. Err.:  0.0158     2 x log-likelihood:  -2782.1530    the dispersion parameter of 0.2104 is getting me little confused...is it supposed to be a good thing? or is this a bad thing?   thank you
122808
I conducted a glm.nb by   glm1&lt;-glm.nb(x~factor(group))   with group being a categorial and x being a metrical variable. When I try to get the summary of the results, I get slightly different results, depending on if I use summary() or summary.glm. summary(glm1) gives me       ... Coefficients:                     Estimate Std. Error z value Pr(&gt;|z|)       (Intercept)       0.1044     0.1519   0.687   0.4921       factor(gruppe)2   0.1580     0.2117   0.746   0.4555       factor(gruppe)3   0.3531     0.2085   1.693   0.0904 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      (Dispersion parameter for Negative Binomial(0.7109) family taken to be 1)   whereas summary.glm(glm1) gives me      ... Coefficients:                     Estimate Std. Error t value Pr(&gt;|t|)       (Intercept)       0.1044     0.1481   0.705   0.4817       factor(gruppe)2   0.1580     0.2065   0.765   0.4447       factor(gruppe)3   0.3531     0.2033   1.737   0.0835 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      (Dispersion parameter for Negative Binomial(0.7109) family taken to be 0.9509067)   I understand the meaning of the dispersion parameter, but not of the line   (Dispersion parameter for Negative Binomial(0.7109) family taken to be 0.9509067).   In the handbook it says, it would be the estimated dispersion, but it seems to be a bad estimate, as 0.95 is not close to 0.7109, or is the estimated dispersion something different than the estimated dispersion parameter? I guess, I have to set the dispersion in the summary.nb(x, dispersion=) to something, but I'm not sure, if I have to set the dispersion to 1 (which will yield the same result as summary() or if I should insert an estimate of the dispersion parameter, in this case leading to summary.nb(glm1, dispersion=0.7109) or something else? Or am I fine with just using the summary(glm1)?
75191
In the book An introduction to statistical learning with Applications in R, when testing the relationship between the Response and the Predictors using null hypothesis, it gives a formula: $$F = \frac{(TSS - RSS)/p}{RSS/(n-p-1)}$$ where $$RSS = \sum_{i=1}^{n} (y_i-\hat{y_i})^2$$ $$\hat{y}=\hat{\beta_0} + \hat{\beta_1}x_1 + \dotsb + \hat{\beta_p}x_p$$ $$TSS = \sum_{i=1}^{n} (y_i-\bar{y})^2$$  I googled the F test and learned that the degree of freedom in F statistic means number of values in the final calculation of a statistic that are free to vary.(, )   But, here I cannot get why in this formula the degree of freedom are p and (n-p-1)?
88765
I have a question about the interpretation of the coefficients of an interaction between continuous and categorical variable. here is my model:  model_glm3=glm(cog~lg_hag+race+pdg+sex+as.factor(educa)+(lg_hag:as.factor(educa)),                 data=base_708)  Coefficients:                          Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)               21.4836     2.0698  10.380  &lt; 2e-16 *** lg_hag                     8.5691     3.7688   2.274  0.02334 *   raceblack                 -8.4715     1.7482  -4.846 1.61e-06 *** racemexican               -3.0483     1.7073  -1.785  0.07469 .   racemulti/other           -4.6002     2.3098  -1.992  0.04687 *   pdg                        2.8038     0.4268   6.570 1.10e-10 *** sexfemale                  4.5691     1.1203   4.078 5.15e-05 *** as.factor(educa)2         13.8266     2.6362   5.245 2.17e-07 *** as.factor(educa)3         21.7913     2.4424   8.922  &lt; 2e-16 *** as.factor(educa)4         19.0179     2.5219   7.541 1.74e-13 *** as.factor(educa)5         23.7470     2.7406   8.665  &lt; 2e-16 *** lg_hag:as.factor(educa)2 -21.2224     6.5904  -3.220  0.00135 **  lg_hag:as.factor(educa)3 -19.8083     6.1255  -3.234  0.00129 **  lg_hag:as.factor(educa)4  -8.5502     6.6018  -1.295  0.19577     lg_hag:as.factor(educa)5 -17.2230     6.3711  -2.703  0.00706 ***   lets say the equation of the model is:  E[cog] = a + b1(lg_hag) + b2(educa2*lg_hag) + b3(educa3*lg_hag) + b4(educa4*lg_hag)  + b5(pdg, centered) + other covars, where  b1 = difference in cog  with higher lg_hag among lowest education (coded as 1) b1 + b2 = difference in cog with higher lg_hag among middle education (coded as 2) b1 + b3 = difference in cog with higher lg_hag among high education (coded as 3) b1 + b3 = difference in cog with higher lg_hag among very high education (coded as 4) b5 = difference in cog with each unit increase in pdg   My question is: if my interpretation is right, how to construct confidence intervals for each effect estimate of interactions (e.g: b1+b2) from the confidence intervals of b1 and b2.
83034
I know I'm missing something in my understanding of logistic regression, and would really appreciate any help.  As far as I understand it, the logistic regression assumes that the probability of a '1' outcome given the inputs, is a linear combination of the inputs, passed through an inverse-logistic function. This is exemplified in the following R code:  #create data: x1 = rnorm(1000)           # some continuous variables  x2 = rnorm(1000) z = 1 + 2*x1 + 3*x2        # linear combination with a bias pr = 1/(1+exp(-z))         # pass through an inv-logit function y = pr &gt; 0.5               # take as '1' if probability &gt; 0.5  #now feed it to glm: df = data.frame(y=y,x1=x1,x2=x2) glm =glm( y~x1+x2,data=df,family="binomial")   and I get the following error message:      Warning messages:   1: glm.fit: algorithm did not converge    2: glm.fit: fitted probabilities numerically 0 or 1 occurred    I've worked with R for some time now; enough to know that probably I'm the one to blame.. what is happening here?
60723
I have very recently started learning about Generalised Linear Mixed Models and was using R to explore what difference it makes to treat group membership as either fixed or random effect. In particular, I am looking at the example dataset discussed here:      As outlined in this tutorial, the effect of Doctor ID is appreciable and I was expecting the mixed model with a random intercept to give better results. However, comparing AIC values for the two methods suggest that this model is worse:  &gt; require(lme4) ; hdp = read.csv("http://www.ats.ucla.edu/stat/data/hdp.csv") &gt; hdp$DID = factor(hdp$DID) ; hdp$Married = factor(hdp$Married) &gt; GLM = glm(remission~Age+Married+IL6+DID,data=hdp,family=binomial);summary(GLM)  Call: glm(formula = remission ~ Age + Married + IL6 + DID, family = binomial,  data = hdp)  Deviance Residuals:  Min       1Q   Median       3Q      Max   -2.5265  -0.6278  -0.2272   0.5492   2.7329    Coefficients:               Estimate Std. Error z value Pr(&gt;|z|)     (Intercept) -1.560e+01  1.219e+03  -0.013    0.990     Age         -5.869e-02  5.272e-03 -11.133  &lt; 2e-16 *** Married1     2.688e-01  6.646e-02   4.044 5.26e-05 *** IL6         -5.550e-02  1.153e-02  -4.815 1.47e-06 *** DID2         1.805e+01  1.219e+03   0.015    0.988     DID3         1.932e+01  1.219e+03   0.016    0.987     [...]  DID405       1.566e+01  1.219e+03   0.013    0.990     DID405       1.566e+01  1.219e+03   0.013    0.990     DID406      -2.885e-01  3.929e+03   0.000    1.000     DID407       2.012e+01  1.219e+03   0.017    0.987     --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)  Null deviance: 10353  on 8524  degrees of freedom Residual deviance:  6436  on 8115  degrees of freedom AIC: 7256  Number of Fisher Scoring iterations: 17   &gt; GLMM = glmer(remission~Age+Married+IL6+(1|DID),data=hdp,family=binomial) ; m  Generalized linear mixed model fit by the Laplace approximation  Formula: remission ~ Age + Married + IL6 + (1 | DID)  Data: hdp  AIC  BIC logLik deviance 7743 7778  -3867     7733 Random effects: Groups Name        Variance Std.Dev. DID    (Intercept) 3.8401   1.9596   Number of obs: 8525, groups: DID, 407  Fixed effects:              Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)  1.461438   0.272709   5.359 8.37e-08 *** Age         -0.055969   0.005038 -11.109  &lt; 2e-16 *** Married1     0.260065   0.063736   4.080 4.50e-05 *** IL6         -0.053288   0.011058  -4.819 1.44e-06 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Correlation of Fixed Effects:          (Intr) Age    Marrd1 Age      -0.898               Married1  0.070 -0.224        IL6      -0.162  0.012 -0.033   &gt; extractAIC(GLM) ; extractAIC(GLMM)  [1]  410.000 7255.962 [1]    5.000 7743.188   Thus, my questions are:  (1) Is it appropriate to compare the AIC values provided by the two functions? If so, why does the fixed effect model do better?  (2) What is the best way to identify if fixed or random effects are more important (ie to quantify that the variability due to the doctor is more important than patient characteristics?
19685
I'm quite new on this with binomial data tests, but needed to do one and now I´m not sure how to interpret the outcome. The y-variable, the response variable, is binomial and the explanatory factors are continuous. This is what I got when summarizing the outcome:  glm(formula = leaves.presence ~ Area, family = binomial, data = n)  Deviance Residuals:  Min      1Q  Median      3Q     Max   -1.213  -1.044  -1.023   1.312   1.344    Coefficients:                         Estimate Std. Error z value Pr(&gt;|z|)  (Intercept)           -0.3877697  0.0282178 -13.742  &lt; 2e-16 *** leaves.presence        0.0008166  0.0002472   3.303 0.000956 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)  Null deviance: 16662  on 12237  degrees of freedom Residual deviance: 16651  on 12236  degrees of freedom (314 observations deleted due to missingness) AIC: 16655 Number of Fisher Scoring iterations: 4   There's a number of things I don't get here, what does this really say:                          Estimate Std. Error z value Pr(&gt;|z|)  (Intercept)           -0.3877697  0.0282178 -13.742  &lt; 2e-16 *** leaves.presence        0.0008166  0.0002472   3.303 0.000956 ***   And what does AIC and Number of Fisher Scoring iterations mean?  &gt; fit Call:  glm(formula = Lövförekomst ~ Areal, family = binomial, data = n)  Coefficients: (Intercept)        Areal   -0.3877697    0.0008166    Degrees of Freedom: 12237 Total (i.e. Null);  12236 Residual (314 observations deleted due to missingness) Null Deviance:      16660  Residual Deviance: 16650        AIC: 16650   And here what does this mean:  Coefficients: (Intercept)        Areal   -0.3877697    0.0008166
23908
I have a problem with the normality of some data I have:  I've done a Kolmogorov test which says it isn't normal with p=.0000,  I don't understand: the skewness of my distribution =-.497, and the kurtosis =-0,024  Here is the plot of my distribution which looks very much normal ...  (I have three scores, and each one of this scores isn't normal with a significant p-value for the Kolmogorov test ... I really don't understand )
101070
I read in several sources that the Quasi Poisson model and the Negative Binomial, should produce (on average) the same results. I tried a simple example and, although very close to each other, the numbers are not quite the same.  Does this depends on the numerical algorithm or is there something else I am not considering? Thank you phi = 3 a = 1/50 b = 3 x = 1:100  mu = exp(a*x+b)  y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator  glmQP = glm(y~x, family=quasipoisson) glmNB = glm.nb(y~x)  glmQP glmNB  glmQP$fitted.values glmNB$fitted.values  plot(x, y) lines(x, glmQP$fitted.values, col = "red") lines(x, glmNB$fitted.values, col = &quot;blue&quot;)
119109
If you split traffic on a website in order to conduct an AB test, optimizing/ monitoring conversion, and you split users by evens/ odds, and have the following information before the test is conducted: group     total_recs    n_success   conversion_rate control    1091653       162424      0.148787 treatment  1098513       153955      0.140149  where the control group has a 6.16% relative improved conversion rate over the treatment group for &gt; 1,000,000 samples (even though this is relatively unlikely, since there's no reason odd id users would convert more than evens etc..) If we then run an AB test for 2 weeks and obtain the following results: group     total_recs   n_success    conversion_rate control    46812        7198         0.153764 treatment  47128        7294         0.154770  where the treatment group shows a relative improvment of 0.65%, intuitively, this suggests that the treatment is an improvement. Assume that the null cannot be rejected, and that for whatever reason, no more data can be collected. Is there some way to adjust for this? I.e. normalize current results based on historic evidence? And if so, how can you do this? Additionally; is this something that actually makes sense to do? Instead, if you were to use a Bayesian approach for model comparison, and incorporate MCMC sampling to generate posterior probabilities, how can the Bayes Factors be calculated from the posterior probabilities? Cheers, Appreciate responses!
999
Assuming I need to find the ML estimator for p, p being the chance of success in a Binomial experiment $Bin(N,p)$, I would expect my density function to be:  $$ f(y) = {{N}\choose{y}} p^y(1-p)^{N-y} $$  And so my likelihood function should be:  $$L(p) = \prod_i^n(f(y_i)) =  \prod_i^n({{N}\choose{y_i}})p^{\sum_i^n y_i} (1-p)^{nN-\sum_i^n yi}$$  However, the last exponent, $nN-\sum yi$ seems to be wrong because when I derive the log likelihood and isolate $p$, I get $p = \frac{\sum_i^ny_i}{nN}$, while it should be $p = \frac{\sum_i^ny_i}{N}$. Indeed, when checking online, I find different exponents:   they have $N-\sum yi$, for no obvious reason to me (last time I checked, $\sum_i^n(N-y_i) = nN-\sum y_i$ !)   as well, but there they start from the likelihood function of a Bernoulli experiment. It makes sense that $\sum_i^n(1-y_i) = n-\sum y_i$, but what is more obscure to me is why they take the likelihood function of a Bernoulli while the problem is clearly about a Binomial. I am aware of the link between the two, but not enough to see why their likelihood functions seem to be substitutable to estimate p, especially since it doesn't give me the same result.
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
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
84305
Is there any way to compare the fit of quasi-Poission and negative binomial models? If so, can it be done in R?
46882
I was trying to fit my data into various models and figured out that the fitdistr function from library MASS of R gives me Negative Binomial as the best-fit. Now from the  page, the definition is given as:     NegBin(r,p) distribution describes the probability of k failures and r   successes in k+r Bernoulli(p) trials   with success on the last trial.   Using R to perform model fitting gives me two parameters mean and dispersion parameter. I am not understanding how to interpret these because I cannot see these parameters on the wiki page. All I can see is the following formula:    where k is the number of observations and r=0...n. Now how do I relate these with the parameters given by R? The help file does not provide much information either.   Also, just to say a few words about my experiment: In a social experiment that I was conducting, I was trying to count the number of people each user contacted in a period of 10 days. The population size was 100 for the experiment.   Now, if the model fits the Negative Binomial, I can blindly say that it follows that distribution but I really want to understand the intuitive meaning behind this. What does it mean to say that the number of people contacted by my test subjects follows a negative binomial distribution? Can someone please help clarify this?
47771
In logistic regression with an intercept term and with at least one dependent variable which is categorical, is there a closed form for the variance of the sum of the intercept and the coefficient of the categorical variable, or do you have to sample from a multivariate distribution with the means and variances of the intercept and the coefficient to get a reliable measure of the variance of this sum?  $P(y) = \frac{exp(y)}{1 + exp(y)}$  $y = \beta_{0} + \beta_{1}x + \epsilon$  where x is categorical.  Would the formula for the variance of a sum (of two random variables) be applicable here?  $Var(\beta_{0} + \beta_{1}) = Var(\beta_{0}) + Var(\beta_{1}) + 2Cov(\beta_{0}, \beta_{1})$  The reason I ask, is that in  I got the impression that no closed form existed for the variance in question, and the advice was to sample from a multivariate distribution with the means and variances of $\beta_{0}$ and $\beta_{1}$  set.seed(1) dependent.var &lt;- sample(c(TRUE, FALSE), 100, replace = TRUE, prob = c(0.3, 0.7)) independent.var &lt;- ifelse(dependent.var, sample(c("Red", "Blue"), replace = TRUE,   size = 10, prob = c(0.8, 0.2)), sample(c("Red", "Blue"), size = 10,   replace = TRUE, prob = c(0.4, 0.6))) table(dependent.var, independent.var) ##               independent.var ## dependent.var Blue Red ##         FALSE   42  26 ##         TRUE     7  25 my.fit &lt;- glm(dependent.var ~ 1 + independent.var, family = binomial(logit)) coef(summary(my.fit))                     Estimate Std. Error   z value     Pr(&gt;|z|) (Intercept)        -1.791759  0.4082482 -4.388897 0.0000113927 independent.varRed  1.752539  0.4951042  3.539737 0.0004005255 &gt; vcov(my.fit)                    (Intercept) independent.varRed (Intercept)          0.1666666         -0.1666666 independent.varRed  -0.1666666          0.2451282   The logit of TRUE for a "Red" case is $ \beta_{0}+\beta_{1} \approx -0.039$. Is the variance for this estimate exactly  vcov(my.fit)[1,1] + vcov(my.fit)[2,2] + 2 * vcov(my.fit)[1,2] [1] 0.07846154 ?   Or is this only an approximation, and a more accurate measure is to be found by sampling, e.g.  library(MASS) var(rowSums(mvrnorm(n = 1E7, mu = coef(my.fit), Sigma = vcov(my.fit)))) [1] 0.07842985 ?   In this simple example, the sampling method does not seem to provide more accurate estimates of the variance (using 1E7 samples).   it is stated that "There is a correspondance between the covariance matrix of the fit parameters and Δχ2 confidence regions only for the case of Gaussian uncertainties on the input measurements.". Is that a reason against relying on the closed form above, or is there perhaps another reason for the advice to sample instead of deriving the variance analytically in cases like this?  EDIT: (In response to the answer given by StasK). The advice I originally got was to simulate from the full model, not from the vcov(), so here is the code to simulate from the full model:  library(arm) sim.i &lt;- sim(my.fit, 100000) logit.for.TRUE.red &lt;- sim.i@coef[,1] + sim.i@coef[,2] var(logit.for.TRUE.red) [1] 0.07781206
99864
I have created a mixed distribution model comprising 80% $H_0$ plus 20% $H_1$ to illustrate the link between the expected proportions of true and false positives and negatives in the PDF, CDF and p-value distributions.  When the effect size is small (25% increase in expected correlations vs $H_0$), everything looks consistent:    However, when I increase the effect size to 75% by separating $H_1$ from $H_0$, things don't seem to add up because I get apparent false negatives (grey area) in the p-value distribution, which are no longer there in the PDF and CDF distributions:    Conversely, when I reduce the effect size to 0% (so that $H_1$ = $H_0$), the grey area disappears altogether from the p-value distribution, whereas it is even more prominent in the PDF and CDF distributions:    Now, I'm sure I have done all the maths right and the actual distribution curves are correct and consistent, but perhaps I am misinterpreting the areas between these curves?  What do these four different coloured areas mean in each chart, and how do they relate to the expected proportions of true and false positives and negatives?  How can I show the link between these proportions in the three distributions, and why don't they correspond with each other as shown?  Please help me understand the apparent discrepancy, thanks!
105635
I ran a linear regression of acceptance into college against SAT scores and family / ethnic background. The data are fictional. This is a follow-up on a prior question, already answered. The question focuses in the gathering and interpretation of odds ratios when leaving the SAT scores aside for simplicity.   The variables are Accepted (0 or 1) and Background ("red" or "blue"). I set up the data so that people of "red" background were more likely to get in:  fit &lt;- glm(Accepted~Background, data=dat, family="binomial") exp(cbind(Odds_Ratio_RedvBlue=coef(fit), confint(fit)))                          Odds_Ratio_RedvBlue             2.5 %       97.5 % (Intercept)             0.7088608                     0.5553459   0.9017961 Backgroundred           2.4480042                     1.7397640   3.4595454   Questions:   Is 0.7 the odd ratio of a person of "blue" background being accepted? I'm asking this because I also get 0.7 for "Backgroundblue" if instead I run the following code:  fit &lt;- glm(Accepted~Background-1, data=dat, family="binomial") exp(cbind(OR=coef(fit), confint(fit)))  Shouldn't the odds ratio of "red" being accepted ($\rm Accepted/Red:Accepted/Blue$) just the reciprocal: ($\rm OddsBlue = 1 / OddsRed$)?