Dear Gavin,

glm reported exactly what it noticed, giving a warning that some very  
small fitted probabilities have been found.
However, your data are **not** quasi-separated. The maximum likelihood  
estimates are really those reported by glm.

A first elementary way is to change the tolerance and maximum number  
of iterations in glm and see if you get the same result:
#
 > mod1

Call:  glm(formula = analogs ~ Dij, family = binomial, data = dat,  
control = glm.control(epsilon = 1e-16,      maxit = 1000))

Coefficients:
(Intercept)          Dij
       4.191      -29.388

Degrees of Freedom: 4033 Total (i.e. Null);  4032 Residual
Null Deviance:	    1929
Residual Deviance: 613.5 	AIC: 617.5
#
This is exactly the same fit as the one you have. If separation  
occured the effects ususally diverge as we allow more iterations to  
glm and at some point.

**************
Secondly an inspection of the estimated asymptotic standard errors,  
reveals nothing to worry for.
#
 > summary(mod1)

Call:
glm(formula = analogs ~ Dij, family = binomial, data = dat, control =  
glm.control(epsilon = 1e-16,
     maxit = 1000))

Deviance Residuals:
        Min          1Q      Median          3Q         Max
-1.676e+00  -1.319e-02  -1.250e-04  -1.958e-06   4.104e+00

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   4.1912     0.3248   12.90   <2e-16 ***
Dij         -29.3875     1.9345  -15.19   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

     Null deviance: 1928.62  on 4033  degrees of freedom
Residual deviance:  613.53  on 4032  degrees of freedom
AIC: 617.53

Number of Fisher Scoring iterations: 11
#
If separation occurred the estimated asymptotic standard errors would  
be unnaturally large. This is because, in the case of separation  
(quasi or not) glm would calculate the standard errors taking the sqrt  
of the diagonal elements of minus the hessian of the log-likelihood,  
in a point where the log-likelihood appears to be flat for the given  
tolerance.

**************
To be certain, you could also try fitting with brglm, which is  
guaranteed to give finite estimates, that have bias of smaller order  
than the MLE and compare the results.
#
 > library(brglm)
 > mod.br <- brglm(analogs ~ Dij, data = dat, family = binomial)
 > mod.br

Call:  brglm(formula = analogs ~ Dij, family = binomial, data = dat)

Coefficients:
(Intercept)          Dij
       4.161      -29.188

Degrees of Freedom: 4033 Total (i.e. Null);  4032 Residual
Deviance:	    613.5448
Penalized Deviance: 610.2794 	AIC: 617.5448
#
The estimates are similar a bit shrunk towards the origin which is  
natural for bias removal. If separation occurred, and given the  
previous discussion, the bias-reduced estimates would be considerably  
different than the estimates that glm reports.

**************
Lastly, the more certain way to check for separation is to inspect the  
profiles of the log-likelihood. Vito suggested this but the chosen  
limits for the xval are not appropriate. If separation would occur the  
estimate would be -Inf so that the profiling as done in his email  
should be done starting from example from -40 rather than -20. This  
would reveal that the profile deviance starts increasing again, while  
if separation occured there would be an asymptote on the left. Below I  
give the correct profiles, as reported by profileModel.
 > library(profileModel)
 > pp <- profileModel(mod1, quantile = qchisq(0.95, 1), objective =  
"ordinaryDeviance")
Preliminary iteration .. Done

Profiling for parameter (Intercept) ... Done
Profiling for parameter Dij ... Done
 > plot(pp)
The profiles are quite quadratic. In the case of separation you would  
have seen asymptotes on the left or on the right (see  
help(profileModel) for an example).

**************
It appears that the fitted logistic curve, while steep still has a  
finite gradient, for example, at the LD50 point
 > library(MASS)
 > dose.p(mod)
               Dose          SE
p = 0.5: 0.1426167 0.003646903
When separation occurs the LD50 point cannot be identified (computer  
software would return something with enormous estimated standard error).

In conclusion, if you get data sets that result in large estimated  
effects on the log-odds scale, the above checks can be used to  
convince you whether separation occurred or not. If there is  
separation (not the case in the current example) then, you could use  
an alternative to maximum likelihood for estimation ---such as  
penalized maximum likelihood in brglm--- which always return  finite  
estimates. Though in that case, I suggest you incorporate the  
uncertainty on how large the estimated effects are in having  
confidence intervals with one infinite endpoint, for example  
confidence intervals as in help(profile.brglm).

Hope this helps,

Best wishes,

Ioannis


On 15 Dec 2008, at 18:03, Gavin Simpson wrote:

> Dear List,
>
> Apologies for this off-topic post but it is R-related in the sense  
> that
> I am trying to understand what R is telling me with the data to hand.
>
> ROC curves have recently been used to determine a dissimilarity
> threshold for identifying whether two samples are from the same "type"
> or not. Given the bashing that ROC curves get whenever anyone asks  
> about
> them on this list (and having implemented the ROC methodology in my
> analogue package) I wanted to try directly modelling the probability
> that two sites are analogues for one another for given dissimilarity
> using glm().
>
> The data I have then are a logical vector ('analogs') indicating  
> whether
> the two sites come from the same vegetation and a vector of the
> dissimilarity between the two sites ('Dij'). These are in a csv file
> currently in my university web space. Each 'row' in this file
> corresponds to single comparison between 2 sites.
>
> When I analyse these data using glm() I get the familiar "fitted
> probabilities numerically 0 or 1 occurred" warning. The data do not  
> look
> linearly separable when plotted (code for which is below). I have read
> Venables and Ripley's discussion of this in MASS4 and other sources  
> that
> discuss this warning and R (Faraway's Extending the Linear Model  
> with R
> and John Fox's new Applied Regression, Generalized Linear Models, and
> Related Methods, 2nd Ed) as well as some of the literature on Firth's
> bias reduction method. But I am still somewhat unsure what
> (quasi-)separation is and if this is the reason for the warnings in  
> this
> case.
>
> My question then is, is this a separation issue with my data, or is it
> quasi-separation that I have read a bit about whilst researching this
> problem? Or is this something completely different?
>
> Code to reproduce my problem with the actual data is given below. I'd
> appreciate any comments or thoughts on this.
>
> #### Begin code snippet  
> ################################################
>
> ## note data file is ~93Kb in size
> dat <- read.csv(url("http://www.homepages.ucl.ac.uk/~ucfagls/ 
> dat.csv"))
> head(dat)
> ## fit model --- produces warning
> mod <- glm(analogs ~ Dij, data = dat, family = binomial)
> ## plot the data
> plot(analogs ~ Dij, data = dat)
> fit.mod <- fitted(mod)
> ord <- with(dat, order(Dij))
> with(dat, lines(Dij[ord], fit.mod[ord], col = "red", lwd = 2))
>
> #### End code snippet  
> ##################################################
>
> Thanks in advance
>
> Gavin
> -- 
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
> Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
> ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
> Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
> Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
> UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
> %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

------------------------------------------
Ioannis Kosmidis
D0.05, Dept. of Statistics,
University of Warwick,
Coventry, CV4 7AL, UK

Webpage: http://go.warwick.ac.uk/kosmidis
Voice: +44(0)2476 150778
Fax:  +44(0)2476 524532
------------------------------------------


	[[alternative HTML version deleted]]

