I'm a newbie to R and am hoping someone can show me the error of my ways. I'm getting what appears to be a logistic regression and plot but I wasn't expecting the slope around 25.2 to be so steep and was expecting a more typical logistic curve.  Also, the glm function threw a couple warnings but I don't think it's necessarily indicative of a problem.  Warning messages: 1: glm.fit: algorithm did not converge  2: glm.fit: fitted probabilities numerically 0 or 1 occurred   p1ShapeResponse &lt;- structure(list(Response = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L), Temperature = c(23.74959742, 23.77777778,23.88043478,24.00322061,24.0394525,24.06360709,24.19041868,24.25885668,24.30112721,24.50040258,24.5326087,24.76811594,24.78019324,24.82045089,24.87077295,24.94927536,25.0821256,25.41022544,25.65982287,25.82085346)), .Names = c("Response", "Temperature"), class = "data.frame", row.names = c(NA, -27L)) p1Shape.glm &lt;- glm(Response ~ Temperature, data=p1ShapeResponse, family=binomial(link=logit))  plot(p1ShapeResponse$Temperature, p1ShapeResponse$Response, xlab="Mean Temperature", ylab="") curve(predict(p1Shape.glm, data.frame(Temperature=x), type="resp"), add=TRUE, col="black") mtext("p1 Shape", side=2, line=2, col="blue") axis(side=1,col="black",col.axis="black",las=1) points(p1ShapeResponse$Temperature,fitted(p1Shape.glm),pch=20)   p2ShapeResponse &lt;- structure(list(Response = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L), Temperature = c(23.74577573, 23.77265745, 23.88018433, 23.99923195, 24.00691244, 24.03763441, 24.05683564, 24.19124424, 24.26036866, 24.30261137, 24.4984639, 24.53302611, 24.76728111, 24.7749616, 24.8172043, 24.87096774, 25.4124424, 25.65821813, 25.82334869)), .Names = c("Response", "Temperature"), class = "data.frame", row.names = c(NA, -27L)) p2Shape.glm &lt;- glm(Response ~ Temperature, data=p2ShapeResponse, family=binomial(link=logit))  ## Allow a second plot on the same graph par(new = TRUE)  ## Plot the second plot and put axis scale on right plot(p2ShapeResponse$Temperature, p2ShapeResponse$Response, axes=FALSE, bty="n", pch="+", col="blue", xlab="", ylab="") curve(predict(p2Shape.glm, data.frame(Temperature=x), type="resp"), add=TRUE, col="blue")  ## a little farther out (line=4) to make room for labels mtext("p2 Shape", side=4, line=2, col="blue") axis(side=4, col="blue",col.axis="blue",las=1)
98301
Why does logistic regression not converge for a linearly separable data set?   For linear separable data sets the model parameters go to infinity when mimizing the error function (according to Bishop2006, Pattern recognition and machine learning, section 4.3.2). I don't understand why. It must be connected to setting the score function, i.e. the derivative of the error function, to zero.  The error function is given by $E(\mathbf{w})$ as shown in the figure below. Where $\mathbf{t}$ is the target variable vector of the training set and $\mathbf{w}$ the parameter / coefficient vector. The posterior probabiltiy for class 1 is given by  $$y(n) = p(C_1|\mathbf{x}_n) = \sigma(\mathbf{w}'\mathbf{x}_n).$$ Where sigma is the logistic sigmoid function and $\mathbf{x}_n$ is the feature vector for instance $n$ including the predictors.     My second question to logistic regression is if the to posterior probabilities $p(C_1|\mathbf{x}_n)$ and $p(C_2|\mathbf{x}_n)$ sum to one?
78696
I always have a hard time explaining statistical techniques to audience with no statistical background. If I wanted to explain what GLM is to such audience (without throwing out statistical jargon), what would be the best or most effective way?  I usually explain GLM with three parts -- (1) the random component which is response variable, (2) the systematic component which is linear predictors, and (3) the link function which is the "key" to connecting (1) and (2). Then I would give an example of linear or logistic regression and explain how the link function is selected based on the response variable. Hence it acts as the key connecting two components.
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
114123
I'd like to solve the heteroskedasticity in logistic regression. In my problem, I have two numeric and 23 dummies variables. I tried to transform the two numerical variables using log, min-max normalization and standard normal transformation but the model continues presenting this phenomenon. How to solve this problem?  My R output  Call: glm(formula = TURMA_PROFICIENTE ~ ., family = "binomial", data = treinamento3,      model = T)  Deviance Residuals:      Min       1Q   Median       3Q      Max   -1.5633  -0.6633  -0.4702  -0.2725   3.2180    Coefficients:                                 Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)                   -11.468260   0.234033 -49.003  &lt; 2e-16 *** MODA_ID_DEPENDENCIA_ADM_TURMA   0.207687   0.029116   7.133 9.82e-13 *** TAMANHO_TURMA                   0.025761   0.002113  12.191  &lt; 2e-16 *** PERC_ALUNOS_GOSTAM_MT           0.855038   0.092606   9.233  &lt; 2e-16 *** TX_RESP_Q001B                   0.294212   0.029333  10.030  &lt; 2e-16 *** TX_RESP_Q004S_EM                0.204347   0.087208   2.343 0.019119 *   TX_RESP_Q005                    0.139776   0.012944  10.798  &lt; 2e-16 *** TX_RESP_Q008                    0.073287   0.014984   4.891 1.00e-06 *** TX_RESP_Q010                    0.032345   0.006231   5.191 2.09e-07 *** TX_RESP_Q018                    0.057162   0.020725   2.758 0.005815 **  TX_RESP_Q020                    0.042434   0.017486   2.427 0.015233 *   TX_RESP_Q022C                   0.133927   0.031147   4.300 1.71e-05 *** TX_RESP_Q028                    0.026202   0.014779   1.773 0.076234 .   TX_RESP_Q048                    0.188193   0.022012   8.549  &lt; 2e-16 *** TX_RESP_Q052                    0.239548   0.015695  15.263  &lt; 2e-16 *** TX_RESP_Q054                    0.031970   0.011816   2.706 0.006814 **  TX_RESP_Q060                    0.036555   0.016207   2.255 0.024106 *   TX_RESP_Q074                    0.166943   0.032754   5.097 3.45e-07 *** TX_RESP_Q075                    0.121384   0.033159   3.661 0.000252 *** TX_RESP_Q095                    0.206870   0.023490   8.807  &lt; 2e-16 *** TX_RESP_Q096                    0.328982   0.016370  20.097  &lt; 2e-16 *** TX_RESP_Q098                    0.117467   0.033336   3.524 0.000426 *** TX_RESP_Q099                    0.203174   0.013005  15.622  &lt; 2e-16 *** TX_RESP_Q106                    0.469938   0.022099  21.265  &lt; 2e-16 *** TX_RESP_Q108                    0.047157   0.015743   2.995 0.002740 **  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 39156  on 42108  degrees of freedom Residual deviance: 34932  on 42084  degrees of freedom AIC: 34982  Number of Fisher Scoring iterations: 5   Breush-Pagan test   bptest(fit3)      studentized Breusch-Pagan test  data:  fit3 BP = 3559.6, df = 24, p-value &lt; 2.2e-16   My plot of the fitted_values vs residuals
43807
I am running a binomial logistic regression with a logit link function in R. My response is factorial [0/1] and I have two multilevel factorial predictors - let's call them $a$ and $b$ where $a$ has 4 factor levels  $(a_1,a_2,a_3,a_4)$ and $b$ has 9 factor levels $(b_1,b_2,\dotsc,b_9)$. Therefore:  mod &lt;- glm(y ~ a + b, family = binomial(logit), data=pretend) summary(mod)   The model output would then show all the information about the model as well as the coefficients.  There is a factor level for both a and b missing (a1 and b1) from the summary output. I understand that it is constrained in the "intercept" of the model. I have read that if I want to remove the intercept term and see the estimates for those factor levels I can just add -1 to the model formula, i.e.:  mod2 &lt;- glm(y ~ a + b - 1, family=binomial(logit), data=pretend) summary(mod2)   In the new model (mod2) the intercept term is then gone and variable a's factor-level a1 is given amongst the list of coefficients. But, variable b's factor-level b1 is still missing and given that there is no intercept term anymore, how can I interpret the odds-ratio for that factor level then?  Could someone please explain to me how to get the coefficient for b1 too and why this is happening?
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?
29900
Specifically, I want to know if there is a difference between lm(y ~ x1 + x2) and glm(y ~ x1 + x2, family=gaussian). I think that this particular case of glm is equal to lm. Am I wrong?
99401
As far as I know, the difference between logistic model and fractional response model (frm) is that the dependent variable (Y) in which frm is [0,1], but logistic is {0, 1}. Further, frm uses the quasi-likelihood estimator to determine its parameters. Normally, we can use glm to obtain the logistic models by glm(y ~ x1+x2, data = dat, family = binomial(logit)). For frm, we change family = binomial(logit) to family = quasibinomial(logit). I noticed we can also use family = binomial(logit) to obtain frm's parameter since it gives the same estimated values. See the following example library(foreign) mydata &lt;- read.dta(&quot;k401.dta&quot;)   glm.bin &lt;- glm(prate ~ mrate + age + sole + totemp,                 data = mydata, family = binomial('logit')) summary(glm.bin)  return,  Call: glm(formula = prate ~ mrate + age + sole + totemp,      family = binomial(&quot;logit&quot;),      data = mydata)  Deviance Residuals:      Min       1Q   Median       3Q      Max   -3.1214  -0.1979   0.2059   0.4486   0.9146    Coefficients:               Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)  1.074e+00  8.869e-02  12.110  &lt; 2e-16 *** mrate        5.734e-01  9.011e-02   6.364 1.97e-10 *** age          3.089e-02  5.832e-03   5.297 1.17e-07 *** sole         3.636e-01  9.491e-02   3.831 0.000128 *** totemp      -5.780e-06  2.207e-06  -2.619 0.008814 **  --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for binomial family taken to be 1)      Null deviance: 1166.6  on 4733  degrees of freedom Residual deviance: 1023.7  on 4729  degrees of freedom AIC: 1997.6  Number of Fisher Scoring iterations: 6  And for family = quasibinomial('logit'), glm.quasi &lt;- glm(prate ~ mrate + age + sole + totemp,   data = mydata ,family = quasibinomial('logit')) summary(glm.quasi)  return,  Call: glm(formula = prate ~ mrate + age + sole + totemp,      family = quasibinomial(&quot;logit&quot;),      data = mydata)  Deviance Residuals:      Min       1Q   Median       3Q      Max   -3.1214  -0.1979   0.2059   0.4486   0.9146    Coefficients:               Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)  1.074e+00  4.788e-02  22.435  &lt; 2e-16 *** mrate        5.734e-01  4.864e-02  11.789  &lt; 2e-16 *** age          3.089e-02  3.148e-03   9.814  &lt; 2e-16 *** sole         3.636e-01  5.123e-02   7.097 1.46e-12 *** totemp      -5.780e-06  1.191e-06  -4.852 1.26e-06 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  (Dispersion parameter for quasibinomial family taken to be 0.2913876)      Null deviance: 1166.6  on 4733  degrees of freedom Residual deviance: 1023.7  on 4729  degrees of freedom AIC: NA  Number of Fisher Scoring iterations: 6  The estimated Beta from both `family` are the same, but the difference is the SE values.  However, to obtain the correct SE, we have to use `library(sandwich)` as in this [post][1].  Now, my questions:  1. What is the difference between these two codes? 2. Is frm about to obtain robust SE?  If my understanding is not correct, please give some suggestions.      [1]: https://stackoverflow.com/questions/37584715/fractional-response-regression-in-r
97082
I have fitted random coefficient Poisson analysis in R. I have obtained the following results:   Generalized linear mixed model fit by maximum likelihood ['glmerMod'] Family: poisson ( log )  Formula: frequency ~ 1 + cc + ageveh + make + (1 | AREA)      Data: x        AIC       BIC    logLik  deviance  1359.1477 1389.7370 -672.5739 1345.1477   Random effects:  Groups Name        Variance Std.Dev.  AREA   (Intercept) 1.323    1.15   Number of obs: 584, groups: AREA, 8  Fixed effects:             Estimate Std. Error z value Pr(&gt;|z|)   (Intercept) -0.12902    0.44432  -0.290   0.7715  ccL          0.05656    0.12371   0.457   0.6475 agevehO      0.02136    0.09264   0.231   0.8177 make2       -0.45454    0.20632  -2.203   0.0276 * make3       -0.31799    0.21422  -1.484   0.1377  make4       -0.29708    0.14469  -2.053   0.0401 * --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Correlation of Fixed Effects:         (Intr) ccL    agevhO make2  make3 ccL      0.052          agevehO -0.179 -0.232                make2   -0.171 -0.007 -0.001  make3   -0.156  0.022 -0.078  0.366      make4   -0.300 -0.235  0.167  0.544  0.522   However I am unable to interpret the results.
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
7359
I am using glmnet where my dependent variable is binary (class 0, class 1). I want to report percentage accuracy of the model. So I use the predict function for my test dataset. However, the values returned are in decimals instead of being 0 and 1. So I set a threshold of 0.5, meaning if the predicted value > 0.5, I consider it as 1 and if predicted value &lt;= 0.5, I consider it as 0. Next I create a confusion matrix by comparing the predicted and actual values of my test data. From this I find the accuracy. I have pasted my sample code below. I am not sure if this is the right approach for reporting the accuracy percentage for a glmnet model predicting a binary dependent variable.  data &lt;- read.csv('datafile', header=T) mat  &lt;- as.matrix(data) X    &lt;- mat[, c(1:ncol(mat)-1)] y    &lt;- mat[, ncol(mat)]  fit  &lt;- cv.glmnet(X, y, family="binomial", type.measure="class", alpha=0.1)  t                             &lt;- 0.2*nrow(mat) #20% of data t                             &lt;- as.integer(t)  testX                         &lt;- mat[1:t, 1:ncol(mat)-1] predicted_y                   &lt;- predict(fit, s=0.01, testX, type='response') predicted_y[predicted_y&gt;0.5]  &lt;- 1 predicted_y[predicted_y&lt;=0.5] &lt;- 0 Yactual                       &lt;- mat[1:t, ncol(mat)] confusion_matrix              &lt;- ftable(Yactual, predicted_y) accuracy                      &lt;- 100* (sum(diag(confusion_matrix)) / length(predicted_y))
20940
I used the inbuilt dataset stackloss here. I used following R code for creating the multiple regression model.  &gt; head(stackloss)   Air.Flow Water.Temp Acid.Conc. stack.loss   1        80         27         89         42   2        80         27         88         37   3        75         25         90         37   4        62         24         87         28   5        62         22         87         18   6        62         23         87         18  &gt; fit &lt;- lm(stack.loss ~ .,data=stackloss) &gt; plot(fit)   And when I plot the model, the four following graphs are displayed.      How to interpret these graphs? What are the significance of these graphs? How to use these graphs?
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.
73628
Let's say that there exists some "true" relationship between $y$ and $x$ such that $y = ax + b + \epsilon$, where $a$ and $b$ are constants and $\epsilon$ is i.i.d normal noise. When I randomly generate data from that R code: x &lt;- 1:100; y &lt;- ax + b + rnorm(length(x)) and then fit a model like y ~ x, I obviously get reasonably good estimates for $a$ and $b$.  If I switch the role of the variables as in (x ~ y), however, and then rewrite the result for $y$ to be a function of $x$, the resulting slope is always steeper (either more negative or more positive) than that estimated by the y ~ x regression. I'm trying to understand exactly why that is and would appreciate it if anyone could give me an intuition as to what's going on there.
4511
When I graph a density plot in R, and all the numbers are slightly greater than 0, I get essentially a vertical line at x = 0.  But when all the numbers are exactly equal to 0, I get some sort of bell curve.  Why is that?  It seems counterintuitive.  The command I used to plot the curves was  for (i in 0:60) {     cur_data = subset(data, time == i)     p &lt;- ggplot(cur_data, aes(x=error)) +                 geom_density() +                 theme_bw() +                 xlab(paste("Error distribution (minute ", i, ")", sep="")) +                 xlim(0, 1)     ggsave(...) }   At i = 60, cur_data should be entirely populated by the values 0.0.      (Originally posted on Stack Overflow; was told to post here.)
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
15933
I've been thinking about this problem.  The usual logistic function for modeling binary data is: $$ \log\left(\frac{p}{1-p}\right)=\beta_0+\beta_1X_1+\beta_2X_2+\ldots $$ However is the logit function, which is an S-shaped curve, always the best for modeling the data?  Maybe you have reason to believe your data does not follow the normal S-shaped curve but a different type of curve with domain $(0,1)$.  Is there any research into this?  Maybe you can model it as a probit function or something similar, but what if it is something else entirely?  Could this lead to better estimation of the effects?  Just a thought I had, and I wonder if there is any research into this.