y = cbind(AB=1997, Ab=906, aB=904, ab=32) # Data from Fisher (1925)
fit = vglm(y ~ 1, AB.Ab.aB.ab(link="identity", init.p=0.9), trace=TRUE)
fit = vglm(y ~ 1, AB.Ab.aB.ab, trace=TRUE)
rbind(y, sum(y)*fitted(fit))
Coef(fit) # Estimated p
p = sqrt(4*(fitted(fit)[,4]))
p*p
summary(fit)
Run the code above in your browser using DataLab