ClearAll["Global`*"]
expit[x_]:=1/(1+Exp[-x])
f1[X_?NumericQ]:=expit[b+b1+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
f2[X_?NumericQ]:=expit[b+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
E1=Integrate[f1[X],{X,-10,10}]
E0=Integrate[f2[X],{X,-10,10}]
b=1
b1=1
Plot[E1,{b2,-5,5}]
Plot[E0,{b2,-5,5}]
Plot[E1-E0,{b2,-5,5}]
Plot[(E1/(1-E1))/(E0/(1-E0)),{b2,-5,5}]

Clear[expit,X,b,b1,b2,f1,f2,E0,E1,Biais]


expit[x_]:=1/(1+Exp[-x])
f1[X_?NumericQ]:=expit[b+b1+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
f2[X_?NumericQ]:=expit[b+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]

(*As association with X changes*)
(*l'effet conditionnel est fixe (b1) mais la distribution des issues change (b2)*)
E1=Integrate[f1[X],{X,-10,10}]
E0=Integrate[f2[X],{X,-10,10}]
b=1
b1=1
Biais=(Exp[b1]-((E1/(1-E1))/(E0/(1-E0))))/((E1/(1-E1))/(E0/(1-E0)))*100
MOR=(E1/(1-E1))/(E0/(1-E0))
AbsBias = Exp[b1]-MOR
Plot[Biais,{b2,-5,5},AxesLabel->{"b2","% error"}]
Plot[((E1/(1-E1))/(E0/(1-E0))),{b2,-5,5},AxesLabel->{"b2","marginal odds ratio"}]
Plot[AbsBias,{b2,-5,5},AxesLabel->{"b2","Absolute bias"}]

Clear[E1,E0,b1,b2,b,Biais,f1,f2]
(*As conditional effect changes*)
b=0
b2=2
(*reparametrizing by the conditional odds rather than log-odds b1*)
f1[X_?NumericQ]:=Expb1*Exp[b+b2*X]/(1+Expb1*Exp[b+b2*X])*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
f2[X_?NumericQ]:=expit[b+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
E1=Integrate[f1[X],{X,-10,10}]
E0=Integrate[f2[X],{X,-10,10}]
Biais=(Expb1-((E1/(1-E1))/(E0/(1-E0))))/((E1/(1-E1))/(E0/(1-E0)))*100
Plot[Biais,{Expb1,0.01,5},AxesLabel->{"exp(b1)","% error"}]
Plot[((E1/(1-E1))/(E0/(1-E0))),{Expb1,0.01,5}]

Clear[E1,E0,b1,b2,b,Biais,f1,f2]
(*As baseline risk changes*)
b1=1
b2=2
E1=Integrate[f1[X],{X,-10,10}]
E0=Integrate[f2[X],{X,-10,10}]
f1[X_?NumericQ]:=expit[b+b1+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
f2[X_?NumericQ]:=expit[b+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
Biais=(Exp[b1]-((E1/(1-E1))/(E0/(1-E0))))/((E1/(1-E1))/(E0/(1-E0)))*100
(*Plot[((E1/(1-E1))/(E0/(1-E0))),{b,-5,5},AxesLabel->{b,"marginal odds ratio"}]*)
Plot[Biais,{b,-10,10},AxesLabel->{b,"% error"}]

ClearAll["Global`*"]
expit[x_]:=1/(1+Exp[-x])
(*reparametrizing by the conditional odds rather than log-odds b1*)
f1[X_?NumericQ]:=Expb1*Exp[b+b2*X]/(1+Expb1*Exp[b+b2*X])*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
f2[X_?NumericQ]:=expit[b+b2*X]*1/Sqrt[2*Pi]*Exp[-1/2*X^2]
E1=Integrate[f1[X],{X,-10,10}]
E0=Integrate[f2[X],{X,-10,10}]
b=1
Biais=(Expb1-((E1/(1-E1))/(E0/(1-E0))))/((E1/(1-E1))/(E0/(1-E0)))*100
Plot3D[Biais,{Expb1,0.01,5},{b2,-5,5},AxesLabel->{"exp(b1)","b2","% error"}]
Plot3D[((E1/(1-E1))/(E0/(1-E0))),{Expb1,0.01,5},{b2,-5,5}]