x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
n <- c(59, 60, 62, 56, 63, 59, 62, 61)

p=y/n; z=log(p/(1-p))
m1=lm(z~x,weights=n*p*(1-p))
print(summary(m1))
print("*******************")

m2=glm(cbind(y,n-y)~x,family=binomial)
print(summary(m2))
print("*******************")

negloglik=function(b) sum( -y*(b[1]+b[2]*x) +n*log(1+exp(b[1]+b[2]*x)) )
out=nlm(negloglik,c(-50,20),hessian=TRUE)
print(out$estim)
print(sqrt(diag(solve(out$hess)))) #SE
print(2*(out$min-sum(log(choose(n,y))))+4) #AIC

a=m1$coeff; b=m2$coeff
curve( 1/(1+exp(-(a[1]+a[2]*x))),1.6,2.0)
curve( 1/(1+exp(-(b[1]+b[2]*x))),1.6,2.0,add=TRUE,col=2)
points(x,p)