My evaluation (R-code below) shows biased results and there would be no correct coverage of my 95% CI. Where did I go wrong in my script?

# de model parameters (data generend model)

beta0<-10

beta1<-69.8333333

beta2<- -15.7500000

beta3<-0.9166667

# effect dat moet geschat worden

effect<-beta1

#steekproefgroottes

n<-8 # steekproefgrootte is 8 voor een G-optimal design

#initialisatie van vectoren

beta1.N<-c()

s2.beta1.N<-c()

cnt<-0

#dataframe voor component

efficientie<-data.frame(component=c(10,0,2.45,3.12,7.25,10,0,7.25))

#aantal simulaties

N<-10000

for(i in 1:N){

X1<-c(10,0,2.45,3.12,7.25,10,0,7.25)

Y<-beta0+beta1*X1+beta2*X1^2+beta3*X1^3

#steek alle data in de data matrix

efficientie$component<-c(10,0,2.45,3.12,7.25,10,0,7.25)

efficientie$component2<-X1^2

efficientie$component3<-X1^3

efficientie$score<-Y

#fit het werkingsmodel aan de data

model.fit<-lm(score~component+component2+component3,data=efficientie)

#extraheer de parameterschatting beta1 uit het gefitte model

beta1.N

*<-summary(model.fit)$coeff[2,1]*

#extraheer de geschatte variantie op beta1_hat

s2.beta1.N

#extraheer de geschatte variantie op beta1_hat

s2.beta1.N

*<-summary(model.fit)$coeff[2,2]^2*

#95% betrouwbaarheidsinterval

CI<-confint(model.fit)

#Nagaan van de 95% bedekking

if((effect>CI[2,1])&(effect<CI[2,2]))cnt<-cnt+1

}

mean(beta1.N)

var(beta1.N)

mean(s2.beta1.N)

cnt/N#95% betrouwbaarheidsinterval

CI<-confint(model.fit)

#Nagaan van de 95% bedekking

if((effect>CI[2,1])&(effect<CI[2,2]))cnt<-cnt+1

}

mean(beta1.N)

var(beta1.N)

mean(s2.beta1.N)

cnt/N