(All of the following is done in R, code to reproduce the dataset is given at the end of this post.)  I have a simulated data set, generated in the following way:   Make 10 categories and label them 1-10. Assign a probability value to each of the categories, so that the first two have probability 0, and the rest have probability drawn uniformly from the interval $[0, 1]$. For each category, draw 50 samples from the Bernoulli distribution with success probability equal to the category's assigned probability.   Thus, the dataset has 500 observations, a sample is shown below:  &gt; dt      names     probs response   1:     1 0.0000000        0   2:     1 0.0000000        0   3:     1 0.0000000        0   4:     1 0.0000000        0   5:     1 0.0000000        0  ---                          496:    10 0.9446753        0 497:    10 0.9446753        1 498:    10 0.9446753        1 499:    10 0.9446753        1 500:    10 0.9446753        1   Now, I fit a logistic regression with formula response ~ name.  My understanding of this is that the predicted value assigned to each category is just the mean of the response for that category. This holds true.  However, the regression coefficients are all insignificant:  Call: glm(formula = response ~ names, family = "binomial", data = dt)  Deviance Residuals:       Min        1Q    Median        3Q       Max   -2.53727  -0.45904  -0.00013   0.54922   2.14597    Coefficients:               Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -1.857e+01  9.224e+02  -0.020    0.984 names2       8.728e-11  1.305e+03   0.000    1.000 names3       2.174e+01  9.224e+02   0.024    0.981 names4       1.762e+01  9.224e+02   0.019    0.985 names5       1.762e+01  9.224e+02   0.019    0.985 names6       1.881e+01  9.224e+02   0.020    0.984 names7       2.022e+01  9.224e+02   0.022    0.983 names8       1.637e+01  9.224e+02   0.018    0.986 names9       2.038e+01  9.224e+02   0.022    0.982 names10      2.132e+01  9.224e+02   0.023    0.982  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 692.50  on 499  degrees of freedom Residual deviance: 343.65  on 490  degrees of freedom AIC: 363.65   Why does this happen and what does it mean? I have seen other questions dealing with some insignificant categories, where the advice then is that the significance of addition of the categorical variable should be tested using likelihood ratio or a chi-square test. In this case the chi-square test shows that adding names to the null model improves the fit.    Code to generate data and models in R:  library(dplyr) library(data.table)   df &lt;- data.frame(names = factor(1:10)) set.seed(0) df$probs &lt;- c(0, 0, runif(8, 0, 1))     df$response = lapply(df$probs, function(i){   rbinom(50, 1, i)   })    dt &lt;- data.table(df)  dt &lt;- dt[, list(response = unlist(response)), by = c('names', 'probs')]   lm0 &lt;- glm(data = dt, formula = response ~ 1, family = 'binomial') summary(lm0) lm1 &lt;- glm(data = dt, formula = response ~ names, family = 'binomial') summary(lm1)  anova(lm0, lm1, test = 'Chisq')
17738
I'm using R to run some logistic regression. My variables were continuous, but I used cut to bucket the data. Some particular buckets for these variables always result in dependent variable being equal to 1. As expcted, the coefficient estimate for this bucket is very high, but the p-value is also high. There are about ~90 observations in either these buckets, and around 800 total observations, so I don't think it's a problem of sample size. Also, this variable should not be related to other variables, which would naturally reduce their p-values.  Are there any other plausible explanations for the high p-value?  Example:  myData &lt;- read.csv("application.csv", header = TRUE) myData$FICO &lt;- cut(myData$FICO, c(0, 660, 680, 700, 720, 740, 780, Inf), right = FALSE) myData$CLTV &lt;- cut(myData$CLTV, c(0, 70, 80, 90, 95, 100, 125, Inf), right = FALSE) fit &lt;- glm(Denied ~ CLTV + FICO, data = myData, family=binomial())   Results are something like this:  Deviance Residuals:       Min        1Q    Median        3Q       Max   -1.53831  -0.77944  -0.62487   0.00027   2.09771    Coefficients:                 Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)     -1.33630    0.23250  -5.747 9.06e-09 *** CLTV(70,80]     -0.54961    0.34864  -1.576 0.114930     CLTV(80,90]     -0.51413    0.31230  -1.646 0.099715 .   CLTV(90,95]     -0.74648    0.37221  -2.006 0.044904 *   CLTV(95,100]     0.38370    0.37709   1.018 0.308906     CLTV(100,125]   -0.01554    0.25187  -0.062 0.950792     CLTV(125,Inf]   18.49557  443.55550   0.042 0.966739     FICO[0,660)     19.64884 3956.18034   0.005 0.996037     FICO[660,680)    1.77008    0.47653   3.715 0.000204 *** FICO[680,700)    0.98575    0.30859   3.194 0.001402 **  FICO[700,720)    1.31767    0.27166   4.850 1.23e-06 *** FICO[720,740)    0.62720    0.29819   2.103 0.035434 *   FICO[740,780)    0.31605    0.23369   1.352 0.176236     --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 1037.43  on 810  degrees of freedom Residual deviance:  803.88  on 798  degrees of freedom AIC: 829.88  Number of Fisher Scoring iterations: 16   FICO in the range [0, 660) and CLTV in the range (125, Inf] indeed always results in Denial = 1, so their coefficients are very large, but why are they also "insignificant"?
20499
We're working with some logistic regressions and we have realized that the average estimated probability always equals the proportion of ones in the sample; that is, the average of fitted values equals the average of the sample.  Can anybody explain me the reason or give me a reference where I can find this demonstration?
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.
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)?
53093
I'm working on a multiple logistic regression in R using glm. The predictor variables are continuous and categorical. An extract of the summary of the model shows the following:  Coefficients:                Estimate Std. Error z value Pr(&gt;|z|) (Intercept)   2.451e+00  2.439e+00   1.005   0.3150 Age           5.747e-02  3.466e-02   1.658   0.0973 . BMI          -7.750e-02  7.090e-02  -1.093   0.2743 ... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Confidence intervals:                    2.5 %       97.5 % (Intercept)  0.10969506 1.863217e+03 Age          0.99565783 1.142627e+00 BMI          0.80089276 1.064256e+00 ...   Odd ratios:                   Estimate Std. Error   z value Pr(&gt;|z|) (Intercept)  1.159642e+01  11.464683 2.7310435 1.370327 Age          1.059155e+00   1.035269 5.2491658 1.102195 B            9.254228e-01   1.073477 0.3351730 1.315670 ...   The first output shows that $Age$ is significant. However, the confidence interval for $Age$ includes the value 1 and the odds ratio for $Age$ is very close to 1. What does the significant p-value from the first output mean? Is $Age$ a predictor of the outcome or not?
56845
I have a data set with a sample size over three million numeric values. Close to 20% are either 0 or 1, with the maximum being nearly 18500. So the data is clearly quite heavily positively skewed.  I am trying to categorize some of this data by putting it into bins of equal width, so I decided to try and find the optimal number of bins. Using the Freedman-Diaconis rule it gave me a value of 126044.0262335108, this is clearly a ridiculously large number of bins for the data.   Breaking the set into the Inter-decile range also proved fruitless giving me [0, 1, 1, 2, 3, 5, 8, 17, 47]  Reading elsewhere the square root of the sample size was suggested, this gave 1732.05081 which is more reasonable. However the method is quite crude.  I also looking into Doane's formula given . But reading up on this method it seems to have been based on an incorrect .  How should I deal with this level of skew in the data?  What is the best way to categorize this data?
66847
I was reading . They are performing a logistic regression (with R). At some point they calculate the probability of getting a residual deviance higher than the one they got on a $\chi^2$ distribution with degrees of freedom equal to the degrees of freedom of the model. Copying-pasting from their website...   &gt; glm( cbind(using,notUsing) ~ age + hiEduc + noMore, family=binomial)  Call:  glm(formula = cbind(using, notUsing) ~ age + hiEduc + noMore,            family = binomial)   Coefficients: (Intercept)     age25-29     age30-39     age40-49       hiEduc       noMore       -1.9662       0.3894       0.9086       1.1892       0.3250       0.8330    Degrees of Freedom: 15 Total (i.e. Null);  10 Residual Null Deviance:      165.8  Residual Deviance: 29.92        AIC: 113.4        The residual deviance of 29.92 on 10 d.f. is highly significant:  &gt; 1-pchisq(29.92,10) [1] 0.0008828339       so we need a better model     Why does it make sense to compute 1-pchisq(29.92,10) and why does a low probability indicate that something is going wrong with their model?
17052
I have a data that where the outcome is the proportion of a species observed in an area by a machine on 2 separate days. Since the outcome is a proportion and does not include 0 or 1 I used a beta regression to fit the model. Temperature is used as  an independent variable. Here is some toy R code:  set.seed(1234) library(betareg) d &lt;- data.frame(   DAY = c(1,1,1,1,2,2,2,2),   Proportion = c(.4,.1,.25, .25, .5,.3,.1,.1),   MACHINE = c("A","B","C","D","H","G","K","L"),   TEMPERATURE = c(rnorm(8)*100) ) b &lt;- betareg(Proportion ~ TEMPERATURE,   data= d, link = "logit", link.phi = NULL, type = "ML") summary(b) ## Call: ## betareg(formula = Proportion ~ TEMPERATURE, data = d, link = "logit", link.phi = NULL, type = "ML") ##  ## Standardized weighted residuals 2: ##     Min      1Q  Median      3Q     Max  ## -1.2803 -1.2012  0.3034  0.6819  1.6494  ##  ## Coefficients (mean model with logit link): ##               Estimate Std. Error z value Pr(&gt;|z|)     ## (Intercept) -1.0881982  0.2620518  -4.153 3.29e-05 *** ## TEMPERATURE  0.0003469  0.0023677   0.147    0.884     ##  ## Phi coefficients (precision model with identity link): ##       Estimate Std. Error z value Pr(&gt;|z|)   ## (phi)    9.305      4.505   2.066   0.0389 * ## --- ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1    Above you can see the TEMPERATURE coefficient is .0003469. Exponentiating, exp(.0003469) = 1.000347  Update incorporating replies and comments:  You can see here how increasing temperature by 1 unit from -10 to 10 increase the proportion  nd &lt;- data.frame(TEMPERATURE = seq(-10, 10, by = 1)) nd$Proportion &lt;- predict(b, newdata = nd) nd$proportion_ratio &lt;- nd$Proportion/(1 - nd$Proportion) plot(Proportion ~ TEMPERATURE, data = nd, type = "b")   The interpretation is: A 1-unit change in TEMPERATURE leads to a relative change of 1.000347 ≈0.04%  in the Proportion: $$\frac{\mathrm{E}(\mathtt{Proportion})}{1-\mathrm{E}(\mathtt{Proportion})}$$   The key word there is relative change so when you compare exp(coef(b))[2] to nd$proportion_ratio[2] / nd$proportion_ratio[1]  you'll see they are the same  ## ratio of proportion nd$proportion_ratio[2] / nd$proportion_ratio[1]  exp(coef(b))[2] nd$proportion_ratio[-1] / nd$proportion_ratio[-20]
117961
I know this is a fairly specific R question, but I may be thinking about proportion variance explained, $R^2$, incorrectly. Here goes.  I'm trying to use the R package randomForest. I have some training data and testing data. When I fit a random forest model, the randomForest function allows you to input new testing data to test. It then tells you the percentage of variance explained in this new data. When I look at this, I get one number.  When I use the predict() function to predict the outcome value of the testing data based on the model fit from the training data, and I take the squared correlation coefficient between these values and the actual outcome values for the testing data, I get a different number. These values don't match up.   Here's some R code to demonstrate the problem.  # use the built in iris data data(iris)  #load the randomForest library library(randomForest)  # split the data into training and testing sets index &lt;- 1:nrow(iris) trainindex &lt;- sample(index, trunc(length(index)/2)) trainset &lt;- iris[trainindex, ] testset &lt;- iris[-trainindex, ]  # fit a model to the training set (column 1, Sepal.Length, will be the outcome) set.seed(42) model &lt;- randomForest(x=trainset[ ,-1],y=trainset[ ,1])  # predict values for the testing set (the first column is the outcome, leave it out) predicted &lt;- predict(model, testset[ ,-1])  # what's the squared correlation coefficient between predicted and actual values? cor(predicted, testset[, 1])^2  # now, refit the model using built-in x.test and y.test set.seed(42) randomForest(x=trainset[ ,-1], y=trainset[ ,1], xtest=testset[ ,-1], ytest=testset[ ,1])
109942
Is there i.i.d. assumption on the response variable of logistic regression?  For example, suppose we have $1000$ data points. It seems the response $Y_i$ is coming from a Bernoulli distribution with $p_i=\text{logit}^{-1}(\beta_0+\beta_1 x_i)$. Therefore, we should have $1000$ Bernoulli distributions, with different parameter $p$.  So, they are "independent", but are not "identical".  Am I right?    PS. I learned logistic regression from "machine learning" literature, where we optimize the objective function and check if it is good in testing data, without talking too much about assumptions.   My question started with this post  Where I try to learn more on statistical assumptions.
47334
This question is sort of general and long-winded, but please bear with me.  In my application, I have many datasets, each consisting of ~20,000 datapoints with ~50 features and a single dependent binary variable. I am attempting to model the datasets using regularized logistic regression (R package )  As part of my analysis, I've created residual plots as follows. For each feature, I sort the datapoints according to the value of that feature, divide the datapoints into 100 buckets, and then compute the average output value and the average prediction value within each bucket. I plot these differences.   Here is an example residual plot:    In the above plot, the feature has a range of [0,1] (with a heavy concentration at 1). As you can see, when the feature value is low, the model appears to be biased towards overestimating the likelihood of a 1-output. For example, in the leftmost bucket, the model overestimates the probability by about 9%.  Armed with this information, I'd like to alter the feature definition in a straightforward manner to roughly correct for this bias. Alterations like replacing  $x \rightarrow \sqrt{x}$  or  $x \rightarrow f_a(x) = \cases{a &amp; if $x&lt;a$ \cr x &amp; else}$  How can I do this? I'm looking for a general methodology so that a human could quickly scroll through all ~50 plots and make alterations, and do this for all datasets and repeat often to keep models up-to-date as the data evolves over time.  As a general question, is this even the right approach? Google searches for "logistic regression residual analysis" don't return many results with good practical advice. They seem to be fixated on answering the question, "Is this model a good fit?" and offer various tests like Hosmer-Lemeshow to answer. But I don't care about whether my model is good, I want to know how to make it better!
37401
If I run the R code below, it will generate two independent vectors then test them to see if they are related in some way (i.e. p-value &lt; 0.05).  If I repeat this 1,000 times, then 50 of these (5%) will have a false positive, with a p-value &lt; 0.05. This is a type 1 error.  If we increase sampleSize to 1000, or even 100,000, the results are the same (5% false positives).  I am struggling to understand this, as I would have thought that if we collected enough samples, the chance of a false positive would drop towards 0 (just like correlation).  So, I guess my question is: "How is it possible that the number of false positives, based on p-values generated from comparing two independent datasets, be independent of sample size?".  # R code to demonstrate that with a large dataset, we can still  # get significant p-values, purely by chance.  # Change this to 1000, and we still get the same number of false positives (50, or 5%). sampleSize = 20  cat("Sample size:",sampleSize,"\n")  set.seed(1010093) n=1000 pValues &lt;- rep(NA,n) for(i in 1:n){     y &lt;- rnorm(sampleSize)     x &lt;- rnorm(sampleSize)     pValues[i] &lt;- summary(lm(y ~ x))$coeff[2,4] } # Controls false positive rate fp = sum(pValues &lt; 0.05) cat("Out of ",n," tests, ",fp," had a p-value &lt; 0.05, purely by by chance.\n",sep="")  ----output---- Running "true positives - none.R" ... Sample size: 1000  Out of 1000 tests, 52 had a p-value &lt; 0.05, purely by by chance.
26061
Suppose I have a $2 \times 2$ table that looks like:             Disease       No Disease Treatment         55                67 Control           42                34  I would like to do a logistic regression in R on this table. I understand that the standard way is to use the glm function with a cbind function in the response. In other words, the code looks like: glm(formula = cbind(c(55,67),c(42,34)) ~ as.factor(c(1, 0)), family = binomial())  I am wondering why R requires us to use the cbind function and why simply using proportions is not sufficient. Is there a way to write this out explicitly as a formula? What would it look in the form of: $$ log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1X $$ where $X = 1$ if we have treatment and $X=0$ for control? Right now it seems like I am regressing on a matrix for the dependent value.
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;)
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.
31320
In R I have a categorical variable that I performed logistic regression on and got the following result:  glm(formula = mortality ~ SMOKE, family = binomial, data = c.data)  Deviance Residuals:      Min       1Q   Median       3Q      Max   -0.2155  -0.2155  -0.2155  -0.1860   2.8515    Coefficients:             Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)  -4.0483     0.3189 -12.694   &lt;2e-16 *** SMOKEN        0.2968     0.3559   0.834    0.404     --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 492.45  on 2369  degrees of freedom Residual deviance: 491.72  on 2368  degrees of freedom AIC: 495.72  Number of Fisher Scoring iterations: 6   Is the value for the intercept the same as SMOKEY (has a history of smoking)?
29229
I'm modeling (using GLM) one binary explanatory variable against a binary response variable.  The output looks ok to me, but I'm a biologist with a basic knowledge of stats, I just wanted to ask if there's any rule/caveat/problem for modeling a binary versus a binary?