Can I run a linear probability model or logistic regression if my binary dependent variable only has 4 1's and 150 0's? I am trying to estimate a linear probability model (LPM) / logit model with about 3 covariates. The main dependent variable is an extremely imbalanced binary variable (there are 4 1's and 150 0's in the whole data). Why would this regression not work?

There doesn't have to be a problem. For instance, consider
n = 154
set.seed(11) # I choose this seed by hand to get 3 positives.

x1 = rnorm(n)
x2 = rnorm(n)
x3 = rnorm(n)

a0 = -2.5
a1 = a2 = a3 = 1

z = a0 + a1*x1 + a2*x2 + a3*x3 + rnorm(n)
y = z > 1

sum(y) # 3

mod = glm(y ~ x1 + x2 + x3, binomial(link = "probit"))
summary(mod)

The summary doesn't look too bad:
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -3.7928     1.1936  -3.178  0.00149 **
x1            1.3886     0.6931   2.003  0.04514 * 
x2            1.1549     0.5122   2.255  0.02416 * 
x3            0.5848     0.5449   1.073  0.28314   

