lm and glm function in R I was running a logistic regression in r using glm() as
glm(Y ~X1 + X2 +X3, data = mydata, family = binomial(link = "logit"))

By accident I ran the model using lm instead:
lm(Y ~X1 + X2 +X3, data = mydata, family = binomial(link = "logit"))

I noticed that the coefficients from the model using lm() were a very good approximation to the marginals on the model using glm() (a difference of $0.005$).
Is this by coincidence or can I use the lm() as I specify to estimate the marginals for logistic regressions?
 A: Linear regression (lm in R) does not have link function and assumes normal distribution. It is generalized linear model (glm in R) that generalizes linear model beyond what linear regression assumes and allows for such modifications. In your case, the family parameter was passed to the ... method and passed further to other methods that ignore the not used parameter. So basically, you've run linear regression on your data.
A: If you take a look at the R help documentation you will note that there is no family argument for the lm function.  By definition, lm models (ordinary linear regression) in R are fit using ordinary least squares regression (OLS) which assumes the error terms of your model are normally distributed (i.e. family = gaussian) with mean zero and a common variance.  You cannot run a lm model using other link functions (there are other functions to do that, though if you wanted--you just can't use lm).  In fact, when you try to run the lm code you've presented above, R will generate a warning like this:
> > Warning message: In lm.fit(x, y, offset = offset, singular.ok =
>     > singular.ok, ...) :   extra argument ‘family’ is disregarded.

When you fit your model using glm, on the other hand, you specified that the error terms in your model were binomial using a logit link function.  This essentially constrains your model so that it assumes no constant error variance and it assumes the error terms can only be 0 or 1 for each observation.  When you used lm you made no such assumptions, but instead, your fit model assumed your errors could take on any value on the real number line.  Put another way, lm is a special case of glm (one in which the error terms are assumed normal).  It's entirely possible that you get a good approximation using lm instead of glm but it may not be without problems.  For example, nothing in your lm model will prevent your predicted values from lying outside $y\in [0, 1]$.  So, how would you treat a predicted value of 1.05 for example (or maybe even trickier 0.5)?  There are a number of other reasons to usually select the model that best describes your data, rather than using a simple linear model, but rather than my re-hashing them here, you can read about them in past posts like this one, this one, or perhaps this one.
Of course, you can always use a linear model if you wanted to--it depends on how precise you need to be in your predictions and what the consequences are of using predictions or estimates that might have the drawbacks noted.
