ymatrix = cbind(A = 725, B = 258, AB = 72, O = 1073) # Order matters, not the name
fit = vglm(ymatrix ~ 1, ABO(link = identity), trace = TRUE, cri = "coef")
coef(fit, matrix = TRUE)
Coef(fit) # Estimated pA and pB
rbind(ymatrix, sum(ymatrix) * fitted(fit))
sqrt(diag(vcov(fit)))
Run the code above in your browser using DataLab