As to the first question, regardless of what regression model you choose, logistic, probit, ANOVA, the predicted means of the response on the probability scale will be the exact same values since your single predictor is a grouping variable. So all models will yield identical fit to the response variable in terms of prediction. This changes if you add covariates in an ANCOVA type situation. Here's some sample R code to show predicted probabilities are the same:
set.seed(12345)
# create three category grouping variable
x <- sample(1:3, 100, TRUE)
y <- rep(NA, 100)
# grouping 1 average probability = .2, group 2 = .4, group 3 = .6
y[x == 1] <- rbinom(length(y[x == 1]), 1, .2)
y[x == 2] <- rbinom(length(y[x == 2]), 1, .4)
y[x == 3] <- rbinom(length(y[x == 3]), 1, .6)
x <- factor(x) # make x categorical
# Obtain predicted probabilities
unique(fitted(glm(y ~ x))) # anova
[1] 0.5263158 0.6000000 0.2187500
unique(fitted(glm(y ~ x, binomial))) # logistic
[1] 0.5263158 0.6000000 0.2187500
unique(fitted(glm(y ~ x, binomial(link = "probit")))) # probit
[1] 0.5263158 0.6000000 0.2187500
The next question is about statistical inference. With a binary outcome, your errors and residuals if you check them will neither be normally distributed nor will they have constant variance, so you violate some of the classical assumptions. In practice though, it does not matter. There is a 1972 paper by Glass, Peckham and Sanders that talks about this after a review of the literature. It is also easy to check this by simulation. Your p-values behave like they do when you meet classical assumptions. Moreover, power is not affected, relative to using a logistic regression model. It is useful to know that there is a weighted least squares alternative to ANOVA that is technically correct, but in this situation, it does not matter.
Here's a simple simulation to check whether ANOVA maintains the nominal error rate when the null is true. It's a three group design, balanced, with total sample size of 50.
set.seed(12345)
p <- replicate(500, {
x <- factor(sample(1:3, 50, TRUE))
y <- rbinom(50, 1, .5)
# x and y are uncorrelated
# conduct ANOVA and save p values
summary(aov(y ~ x))[[1]]$`Pr(>F)`[1]
})
# what is the statistics power of the test?
# it should be close to 5% if the null is true
mean(p < .05)
[1] 0.042