I used the function SSfpl (in R) for finding out the correct initial parameters. Why aren't they right? This is my data:

I used SSfpl for finding out the correct initial parameters
 nls(y~ SSfpl (-x, b1,b2,b3,b4), data=estatus) 

And with these initial parameters, I used 
nls(y~b1 + ((b2-b1)/(1 + exp(b3*(x-b4)))), data=estatus, start=emp)
I can't resolve why I have the next error:
Error in nlsModel(formula, mf, start, wts) :   singular gradient matrix at initial parameter estimates

 A: Your problem seems to stem from a typo/mistake in your equation.
Setup:
estatus <- data.frame(x = c(19, 25,38,47,53,69),
                     y = c(65, 61, 56, 28, 12, 10))
plot(y~x,data=estatus)

Let's see what the initial parameter values look like:
library(nlme)
(i1 <- getInitial(y ~ SSfpl(x, A, B, xmid, scal), data = estatus))
##         A         B      xmid      scal 
## 62.923635  9.107685 44.738601  3.363543
## compute curve for initial fit ...
xvec <- 20:70
y1 <- with(as.list(i1),SSfpl(xvec, A, B, xmid, scal))
lines(xvec,y1)


This looks OK.
The easiest and most reliable way to fit the model is to use SSfpl() again: among other things, this also provides gradients for the fitting procedure.
nls(y ~ SSfpl(x, A, B, xmid, scal), data = estatus)
##      A      B   xmid   scal 
## 62.924  9.108 44.739  3.364 

Your equation was b1 + ((b2-b1)/(1 + exp(b3*(x-b4)))): this makes two mistakes (1) it switches the third (xmid) and fourth (scal) parameters; (2) it multiplies rather than dividing by the scal parameter (it also switches the first and second parameters, lower and upper asymptotes, but that seems relatively harmless).  This works fine:
nls(y~B + ((A-B)/(1 + exp((x-xmid)/scal))), data=estatus, start=i1)

