I am running a logistic regression in R, such as: glm(cbind(successes,tries-successes) ~ A + B, family=&quot;binomial&quot;)  A is numeric, B is a factor. A only contains 4 different unique values; 2 of these never lead to 'successes', one value almost never does and the fourth value leads to a lot of successes. All 4 values are on a clearly interpretable linear scale, but it seems there is sort of a threshold value, above which successes can be reached. The factor B has two levels, both of which can lead to successes. When I run the GLM and look at my estimators for the different levels of B, the confidence intervals are huge. If I run the same GLM only including B, the confidence intervals for the estimators are reasonably narrow. If I run the full model, but treating both, A and B as factors, I also get reasonable confidence intervals. This leads me to the conclusion that including A as a numeric fixed effect 'ruins' the model. However, turning it into a factor also doesn't seem right, as the order (and possibly magnitude) does matter. So my questions are:  Would it make sense to treat A as an ordered factor? This looses information on the magnitude of differences on the numeric scale, but at least keeps the ordering. And it seems to be the only way of building a functioning model including both, A and B. Why do I see this weird behaviour that A as a numeric variable 'destroys' the whole model? I find this especially intriguing, as there is no interaction term included. Do you think there is another way of solving the problem? Would it help including B maybe as a numeric column, but within a quadratic term, to acknowledge more this kind of 'threshold' behaviour?
2225
If you have a variable which perfectly separates zeroes and ones in target variable, R will yield the following "perfect or quasi perfect separation" warning message:  Warning message: glm.fit: fitted probabilities numerically 0 or 1 occurred    We still get the model but the coefficient estimates are inflated.   How do you deal with this in practice?
119147
I am fitting a binomial family glm in R, and I have a whole troupe of explanatory variables, and I need to find the best (R-squared as a measure is fine). Short of writing a script to loop through random different combinations of the explanatory variables and then recording which performs the best, I really don't know what to do. And the leaps function from package  does not seem to do logistic regression.  Any help or suggestions would be greatly appreciated.
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$)?
94302
I'm performing a GLM, the response variable is number of individuals, and response variables are    habitat (4 levels) and  season (4 levels).    I need some help since I know the summary() shows p-values but not for the first (alphabetical) level of factor. I mean, I don't know how to interpret that the Intercept has a significant p-value. I can't reach a biological explanation for this model. Hope you can help me ...  Call: glm(formula = individuals ~ habitat + season, family = poisson(link = "log"),  data = data)  Deviance Residuals:   Min        1Q    Median        3Q       Max   -2.85859  -0.51541  -0.08508   0.36497   2.29058    Coefficients:                  Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)            1.6774     0.1624  10.331  &lt; 2e-16 *** habitatDeciduous       0.1340     0.1696   0.790  0.42950     habitatSemiDeciduous   0.2102     0.1675   1.255  0.20933     habitatWetland        -0.1861     0.2039  -0.912  0.36151     seasondry2018         -0.1138     0.1510  -0.753  0.45123     seasonwet2016         -0.2699     0.1576  -1.713  0.08677 .   seasonwet2017         -0.4383     0.1656  -2.647  0.00813 **   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for poisson family taken to be 1)  Null deviance: 78.487  on 63  degrees of freedom Residual deviance: 64.556  on 57  degrees of freedom AIC: 287.77
3015
I have a dataset of around 120000 (120K) unique individuals. I am fitting a binary logistic regression, where I have around 150 variables to choose from. For the categorical variables, some are very imbalanced. For example, there is a variable that can take on Yes (=1) or No (=0). There are 119K individuals with Yes and the remaining 1K with No.    What are the consequences of having such a variable in the logistic regression?  Can it make my results unstable if I have many of these (say 10-20 similar variables), when using stepwise (forward) or backward automatic selection procedures?  What if the predictor / variable is not binary, but has for example 5 categories where 1 category is largely over or underrepresented?
103290
Ok, I have a logistic regression and have used the predict() function to develop a probability curve based on my estimates.   ## LOGIT MODEL: library(car) mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))  ## PROBABILITY CURVE: all.x &lt;- expand.grid(won=unique(won), bid=unique(bid)) y.hat.new &lt;- predict(mod1, newdata=all.x, type="response") plot(bid&lt;-000:1000,predict(mod1,newdata=data.frame(bid&lt;-c(000:1000)),type="response"), lwd=5, col="blue", type="l")   This is great but I'm curious about plotting the confidence intervals for the probabilities. I've tried plot.ci() but had no luck. Can anyone point me to some ways to get this done, preferably with the car package or base R.
3273
I just fit a model in lme4, and I'm wondering what the heck I fit...  I have individuals id, and each is measured pass/fail on items that can be described using two factors, f1 and f2. My theory says that id and f2 can interact. So I want to compare these models:  # f2 doesn't matter m1 &lt;- lmer(pass ~ f1 + (1|id), family="binomial")  # f2 is a fixed effect, doesn't interact with individuals m2 &lt;- lmer(pass ~ f1 + f2 + (1|id), family="binomial")  # f2 is a nested random effect, interacting with individuals m3 &lt;- lmer(pass ~ f1 + (1|f2/id), family="binomial")   First, am I right that there's a random effect of each level of f2, and each individual's level of f2 is partially pooled, depending on the number of observations for that individual?  Second, does it make sense to fit f2 as a fixed effect with random individuals nested within each level of f2? How does one write that?  Third, what's this, with f2 and id reversed in the formula? I typed it by accident, but now I want to understand it.  m3weird &lt;- lmer(pass ~ f1 + (1|id/f2), family="binomial")
71408
Using R, I have fitted a linear model for a single response variable from a mix of continuous and discrete predictors. This is uber-basic, but I'm having trouble grasping how a coefficient for a discrete factor works.   Concept: Obviously, the coefficient of the continuous variable 'x' is applied in the form y = coefx(varx) + intercept but how does that work for a factor z if the factor is non-numeric? y = coefx(varx) + coefz(factorz???) + intercept   Specific: I have fitted a model in R as lm(log(c) ~ log(d) + h + a + f + h:a) where h and f are discrete, non-numeric factors. The coefficients are:  Coefficients:               Estimate  (Intercept)  -0.679695  log(d)        1.791294  h1            0.870735   h2           -0.447570   h3            0.542033    a             0.037362   f1           -0.588362   f2            0.816825  f3            0.534440 h1:a         -0.085658 h2:a         -0.034970  h3:a         -0.040637   How do I use these to create the predictive equation:  log(c) =  1.791294(log(d)) + 0.037362(a) + h??? + f???? + h:a???? + -0.679695   Or am I doing it wrong?   I THINK that that concept is if the subject falls in category h1 and f2, the equation becomes:   log(c) =  1.791294(log(d)) + 0.037362(a) +  0.870735  + 0.816825  + h:a???? + -0.679695   But I'm really not clear on how the h:a interactive term gets parsed. Thanks for going easy on me.
28205
In answering  question John Christie suggested that the fit of logistic regression models should be assessed by evaluating the residuals.  I'm familiar with how to interpret residuals in OLS, they are in the same scale as the DV and very clearly the difference between y and the y predicted by the model.  However for logistic regression, in the past I've typically just examined estimates of model fit, e.g. AIC, because I wasn't sure what a residual would mean for a logistic regression.  After looking into  a little bit I see that in R there are five types of glm residuals available, c("deviance", "pearson", "working","response", "partial").  The help file refers to:     Davison, A. C. and Snell, E. J. (1991) , eds. Hinkley, D. V., Reid, N. and Snell, E. J., Chapman &amp; Hall.     I do not have a copy of that.  Is there a short way to describe how to interpret each of these types?  In a logistic context will sum of squared residuals provide a meaningful measure of model fit or is one better off with an Information Criterion?
62405
I want to perform logistic regression with the following binomial response and with $X_1$ and $X_2$ as my predictors.      I can present the same data as Bernoulli responses in the following format.    The logistic regression outputs for these 2 data sets are mostly the same.  The deviance residuals and AIC are different.  (The difference between the null deviance and the residual deviance is the same in both cases - 0.228.)  The following are the regression outputs from R. The data sets are called binom.data and bern.data.  Here is the binomial output.  Call: glm(formula = cbind(Successes, Trials - Successes) ~ X1 + X2,      family = binomial, data = binom.data)  Deviance Residuals:  [1]  0  0  0  Coefficients:             Estimate Std. Error z value Pr(&gt;|z|) (Intercept)  -2.9649    21.6072  -0.137    0.891 X1Yes        -0.1897     2.5290  -0.075    0.940 X2            0.3596     1.9094   0.188    0.851  (Dispersion parameter for binomial family taken to be 1)  Null deviance:  2.2846e-01  on 2  degrees of freedom Residual deviance: -4.9328e-32  on 0  degrees of freedom AIC: 11.473  Number of Fisher Scoring iterations: 4   Here is the Bernoulli output.  Call: glm(formula = Success ~ X1 + X2, family = binomial, data = bern.data)  Deviance Residuals:      Min       1Q   Median       3Q      Max   -1.6651  -1.3537   0.7585   0.9281   1.0108    Coefficients:             Estimate Std. Error z value Pr(&gt;|z|) (Intercept)  -2.9649    21.6072  -0.137    0.891 X1Yes        -0.1897     2.5290  -0.075    0.940 X2            0.3596     1.9094   0.188    0.851  (Dispersion parameter for binomial family taken to be 1)  Null deviance: 15.276  on 11  degrees of freedom Residual deviance: 15.048  on  9  degrees of freedom AIC: 21.048  Number of Fisher Scoring iterations: 4   My questions:   I can see that the point estimates and standard errors between the 2 approaches are equivalent in this particular case.  Is this equivalence true in general? How can the answer for Question #1 be justified mathematically?   Why are the deviance residuals and AIC different?
12043
I have built a binomial regression model using the equation: logit_1 &lt;- glm(past_adl_binary~., family = binomial(link = &quot;cloglog&quot;),                 data = newdata)  The outcome of the model is as follows:  Can anyone help me interpret the outcome of this model? I understand that all of the coefficients are significant apart from sunangle. Can I use this equation (100*exp(parameter estimate)-1) to calculate the increase/ decrease in percent probability of my y-variable (1 or 0) changing based on my predictor variables?  For example, for habitat_categoryONDEEP: 100*(exp(-0.173)-1) = -15.886.  Would this suggest that a habitat category of ONDEEP decreases the probability of my y-variable being 1 by 15.9%, compared to the reference category?
98436
First off, I'll say I am a biologist and new to the statistics side of things so excuse my ignorance  I have a data set that consists of a binary outcome and then a bunch of trinary explanatory variables that looks something like this:  head()  Category block21_hap1 block21_hap2 block21_hap3 block21_check 1        1            1            1            0             2 2        1            2            0            0             2 3        1            1            0            1             2 4        1            1            0            1             2 5        1            1            1            0             2 6        1            1            1            0             2   A quick summary of the data  summary() Category block21_hap1 block21_hap2 block21_hap3 block21_check  1:718    0:293        0:777        0:1026       2:1467         0:749    1:709        1:577        1: 390                               2:465        2:113        2:  51     and another summary grouped by outcome levels  by(hap.ped.final, hap.ped.final$Category, summary) hap.ped.final$Category: 1  block21_hap1 block21_hap2 block21_hap3 block21_check  0:146        0:374        0:518        2:718          1:336        1:286        1:174                       2:236        2: 58        2: 26                      ----------------------------------------------------------------------------  hap.ped.final$Category: 0  block21_hap1 block21_hap2 block21_hap3 block21_check  0:147        0:403        0:508        2:749          1:373        1:291        1:216                       2:229        2: 55        2: 25             So I am trying to run logistic regression on this data. When I do this:  fit = glm(Category~ block21_hap1 + block21_hap2 + block21_hap3, data = hap.ped.final ,family = "binomial") summary(fit)  Deviance Residuals:     Min      1Q  Median      3Q     Max   -1.301  -1.177   1.059   1.177   1.200    Coefficients: (1 not defined because of singularities)                              Estimate Std. Error z value Pr(&gt;|z|) (Intercept)                 -0.039221   0.280110  -0.140    0.889 hap.ped.final$block21_hap11  0.123555   0.183087   0.675    0.500 hap.ped.final$block21_hap12  0.009111   0.295069   0.031    0.975 hap.ped.final$block21_hap21 -0.084334   0.183087  -0.461    0.645 hap.ped.final$block21_hap22 -0.013889   0.337468  -0.041    0.967 hap.ped.final$block21_hap31  0.201113   0.183087   1.098    0.272 hap.ped.final$block21_hap32        NA         NA      NA       NA  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 2033  on 1466  degrees of freedom Residual deviance: 2028  on 1461  degrees of freedom AIC: 2040  Number of Fisher Scoring iterations: 3   So I don't really know what a singularity is or what's going wrong here that is throwing up NA's as a result of my analysis. Is it my data, or what I'm doing to it. I tried googling the warning (or whatever you might call it) and I got some pages talking about collinearity and multilinearity, which I do not understand at all.  Again, sorry for lack of knowledge here. I wish I had done more maths in undergrad.
93158
I've conducted a psychological experiment on the same subject, under two different condition. For each condition I've collected the number of correct and wrong answer for each stimulus (number of trials per stimulus = 10; number of stimulus = 15). From the data collected I fitted a psychometric curve (probit or logit curve), and now I want to compare the results.  Consider subject 1: His answers are collected in the following matrices  **COND_1**     correct cnt   Y  N   1   0 10   2   0 10   3   0 10   4   0 10   5   0 10   6   0 10   7   1  9   8   4  6   9   5  5   10  7  3   11  8  2   12  7  3   13 10  0   14 10  0   15 10  0  **COND_2**     correct cnt   Y  N   1   0 10   2   2  8   3   2  8   4   3  7   5   2  8   6   4  6   7   8  2   8  10  0   9   8  2   10 10  0   11  9  1   12 10  0   13 10  0   14 10  0   15 10  0  COND_1.1 &lt;- matrix(c(0, 2, 2, 3, 2, 4, 8, 10, 8, 10, 9, 10, 10, 10, 10, 10, 8, 8, 7, 8, 6, 2, 0, 2, 0, 1, 0, 0, 0, 0), byrow=F, ncol=2)  COND_1.2 &lt;- matrix(c(0, 0, 0, 0, 0, 0, 1, 4, 5, 7, 8, 7, 10, 10, 10, 10, 10, 10, 10, 10, 10, 9, 6, 5, 3, 2, 3, 0, 0, 0), byrow=F, ncol=2)  cnt &lt;- seq(from=0, to=1.4, by=0.1)  ddprob.1.1 &lt;- glm(COND_1.1 ~ cnt, family = binomial(link = "probit"))  ddprob.2.1 &lt;- glm(COND_2.1 ~ cnt, family = binomial(link = "probit"))     Observing the plot, you can surely say that performance under condition 2 is better then the performance under condition 1. In fact if you compute some parameters like AUC, p25, p50, p75 (the value of the stimulus at which the number of correct answers are 25%, 50% and 75% each) [threshold parameters], you'll see that:  COND_1      COND_2 AUC    &lt;    AUC p25    &gt;    p25 p50    &gt;    p50 p75    &gt;    p75     Now consider subject 2:  COND_1     correct cnt   Y  N   1   0 10   2   0 10   3   0 10   4   0 10   5   3  7   6   3  7   7   4  6   8   5  5   9   9  1   10  9  1   11 10  0   12 10  0   13 10  0   14 10  0   15 10  0  COND_2     correct cnt   Y  N   1   0 10   2   1  9   3   0 10   4   4  6   5   2  8   6   6  4   7   4  6   8   7  3   9   4  6   10  7  3   11  7  3   12 10  0   13  9  1   14 10  0   15 10  0   COND_1.2 &lt;- matrix(c(0, 0, 0, 0, 3, 3, 4, 5, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 7, 7, 6, 5, 1, 1, 0, 0, 0, 0, 0), byrow=F, ncol=2) COND_2.2 &lt;- matrix(c(0, 1, 0, 4, 2, 6, 4, 7, 4, 7, 7, 10, 9, 10, 10, 10, 9, 10, 6, 8, 4, 6, 3, 6, 3, 3, 0, 1, 0, 0), byrow=F, ncol=2) cnt &lt;- seq(from=0, to=1.4, by=0.1)  ddprob.1.1 &lt;- glm(COND_1.2 ~ cnt, family = binomial(link = "probit")) ddprob.2.1 &lt;- glm(COND_2.2 ~ cnt, family = binomial(link = "probit"))   The plotted curves are shown in the image below.    Can you say what performance is better? Under condition 1 or under condition 2? Extrapolating AUC, p25, p50, p75, you have:  COND_1      COND_2 AUC     &gt;   AUC p25     &gt;   p25 p50     =   p50 p75     &lt;   p75   So my question is: Is there a method to say that a generic curve (logit, probit or whatelse) is HIGHER then another? Is there a method to compare (a single number of "performance") that describe the differences of the two curves? My example shows that AUC, p25, p50, p75 are not good parameters. I'd like to compute a single numeric parameter for each curves, to make a simple t-paired test, extrapolated from the distribution seen.
17662
If I have a dataset with a very rare positive class, and I down-sample the negative class, then perform a logistic regression, do I need to adjust the regression coefficients to reflect the fact that I changed the prevalence of the positive class?  For example, let's say I have a dataset with 4 variables: Y, A, B and C.  Y, A, and B are binary, C is continuous.  For 11,100 observations Y=0, and for 900 Y=1:  set.seed(42) n &lt;- 12000 r &lt;- 1/12 A &lt;- sample(0:1, n, replace=TRUE) B &lt;- sample(0:1, n, replace=TRUE) C &lt;- rnorm(n) Y &lt;- ifelse(10 * A + 0.5 * B + 5 * C + rnorm(n)/10 &gt; -5, 0, 1)   I fit a logistic regression to predict Y, given A, B and C.   dat1 &lt;- data.frame(Y, A, B, C) mod1 &lt;- glm(Y~., dat1, family=binomial)   However, to save time I could remove 10,200 non-Y observations, giving 900 Y=0, and 900 Y=1:  require('caret') dat2 &lt;- downSample(data.frame(A, B, C), factor(Y), list=FALSE) mod2 &lt;- glm(Class~., dat2, family=binomial)   The regression coefficients from the 2 models look very similar:  &gt; coef(summary(mod1))               Estimate Std. Error   z value     Pr(&gt;|z|) (Intercept) -127.67782  20.619858 -6.191983 5.941186e-10 A           -257.20668  41.650386 -6.175373 6.600728e-10 B            -13.20966   2.231606 -5.919353 3.232109e-09 C           -127.73597  20.630541 -6.191596 5.955818e-10 &gt; coef(summary(mod2))               Estimate  Std. Error     z value    Pr(&gt;|z|) (Intercept) -167.90178   59.126511 -2.83970391 0.004515542 A           -246.59975 4059.733845 -0.06074284 0.951564016 B            -16.93093    5.861286 -2.88860377 0.003869563 C           -170.18735   59.516021 -2.85952165 0.004242805   Which leads me to believe that the down-sampling did not affect the coefficients.  However, this is a single, contrived example, and I'd rather know for sure.
58580
I am trying to produce a model for which I have a response variable which is a proportion between 0 and 1, this includes quite a few 0s and 1s but also many values in between. I am thinking about attempting a beta regression. The package I have found for R (betareg) only allows values  in between 0 and 1 but not including 0 or 1 them selves. I have read elsewhere that theoretically the beta distribution should be able to handle values of 0 or 1 but I do not know how to handle this in R.I have seen some people add 0.001 to the zeros and take 0.001 from the ones, but I am not sure this is a good idea?  Alternatively I could logit transform the response variable and use linear regression. In this case I have the same problem with the 0 and 1's which cannot be log transformed.
94380
I've compared the logistic regression models on R (glm) and on Spark (LogisticRegressionWithLBFGS) on a dataset of 390 obs. of 14 variables.  The results are completely different in the intercept and the weights. How to explain this?  Here is the results of Spark (LogisticRegressionWithLBFGS) :  model.intercept  :   1.119830027739959 model.weights :  GEST    0.30798496002530473  DILATE  0.28121771009716895  EFFACE  0.01780105068588628  CONSIS -0.22782058111362183  CONTR  -0.8094592237248102  MEMBRAN-1.788173534959893  AGE    -0.05285751197750732  STRAT  -1.6650305527536942  GRAVID  0.38324952943210994  PARIT  -0.9463956993328745  DIAB   0.18151162744507293  TRANSF -0.7413500749909346  GEMEL  1.5953124037323745   Here is the result of R :               Estimate Std. Error z value Pr(&gt;|z|)    (Intercept)  3.0682091  3.3944407   0.904 0.366052     GEST         0.0086545  0.1494487   0.058 0.953821     DILATE       0.4898586  0.2049361   2.390 0.016835 *   EFFACE       0.0131834  0.0059331   2.222 0.026283 *   CONSIS       0.1598426  0.2332670   0.685 0.493196     CONTR        0.0008504  0.5788959   0.001 0.998828     MEMBRAN     -1.5497870  0.4215416  -3.676 0.000236 ***    AGE         -0.0420145  0.0326184  -1.288 0.197725     STRAT       -0.3781365  0.5860476  -0.645 0.518777     GRAVID       0.1866430  0.1522925   1.226 0.220366     PARIT       -0.6493312  0.2357530  -2.754 0.005882 **   DIAB         0.0335458  0.2163165   0.155 0.876760     TRANSF      -0.6239330  0.3396592  -1.837 0.066219 .   GEMEL        2.2767331  1.0995245   2.071 0.038391 *   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
32113
We have comprehensive coverage on the topic of perfect separation in logistic regression. When it happens in R we usually see two warnings:     Warning messages:   1: glm.fit: algorithm did not converge   2: glm.fit: fitted probabilities numerically 0 or 1 occurred    Is it possible that the algorithm did converge and fitted probabilities numerically 0 or 1 occurred nonetheless?