# Why is my quasibinomial GLM estimator biased - Monte Carlo simulation

I'm playing with some Monte Carlo simulations to get an idea of the properties of some linear and non-linear models. The linear OLS model in my case is specified as:

$$Y_t = \beta_0 + \beta_1x+ \epsilon,$$

with $$\beta_0=\beta_1=0.5$$

My point of interest here is to figure out if the estimates $$\hat{\beta_0}$$ and $$\hat{\beta_1}$$ are unbiased and roughly equal to the true parameters $$\beta_0$$ and $$\beta_1$$.

Below, you can see my R code for doing a Monte Carlo simulation with both the linear OLS model and a non-linear quasibinomial probit GLM version of it.

#Monte Carlo simulation

N <- 100 #number of observations
mean <- 0 #mean
sd <- 1 #standard deviation
x <- rnorm(N,mean = mean,sd=sd) #simulate random exogeneous variables

n_sim <- 1000 #number of simulations

b0 <- 0.5 #true parameter beta 0
b1 <- 0.5 #true parameter beta 1

b0_lin <- numeric(n_sim) #stores the estimates of b0
b1_lin <- numeric(n_sim) #stores the estimates of b1
b0_probit <- numeric(n_sim) #stores the estimates of b0
b1_probit <- numeric(n_sim) #stores the estimates of b1

#apply the Monte Carlo simulation
for (i in 1:n_sim){
set.seed(1+i)
#simulate dependent variable
y_lin <- b0+b1*x+rnorm(N, mean=mean, sd=sd)
#linear model
mod <- lm(y_lin~x)
#save vector of estimated parameters
b0_lin[i] <- coef(mod)[]
b1_lin[i] <- coef(mod)[]

#transform dependend variable with inverse probit link funktion, thus it lies between zero and 1
y_prob <- pnorm(y_lin)
#quasibinomial probit GLM model
mod <- glm(y_prob~x, family = quasibinomial("probit"))
#save vector of estimated parameters
b0_probit[i] <- coef(mod)[]
b1_probit[i] <- coef(mod)[]
}

#Save means of the parameter vectors
b0_lin_mean <- mean(b0_lin)
b1_lin_mean <- mean(b1_lin)

b0_probit_mean <- mean(b0_probit)
b1_probit_mean <- mean(b1_probit)

Montecarlo_comparison <- data.frame("b0_mean"=c("True parameters"=b0,
"Linear"=b0_lin_mean,
"Probit"=b0_probit_mean),
"b1_mean"=c(b1,
b1_lin_mean,
b1_probit_mean)
)
print(Montecarlo_comparison)


As you can see below, the estimator from a linear OLS model is not biased. However, the quasibinomial probit GLM estimator appears to be severely biased. My questions: Why is the quasibinomial probit GLM estimator biased in my Monte Carlo simulation? Is the process of my Monte Carlo simulation technically correct or did I maybe make a mistake in the model specifications? If I did it correct, what is the explanation for the fact that the quasibinomial probit GLM estimator is biased here?