Error in R-code for evaluation of optimal design

#1
For my project I have to choose an optimal design for a given problem and I have to evaluate this design. The problem follows a linear regression model (cubic relationship), with a variance of 23. Out of the restrictions that were given, I know the 4 coefficients of the parameters.

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<-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
 

rogojel

TS Contributor
#2
Hi,
I may not understand the code fully, but it seems to me that you are calculating N=10000 times the exactly same model, with no random variations whatsoever. Do I miss something?

regards
 
#3
Hi,

Thanks for your answer! Today I also concluded that this lus is always calculating the same model. But I'm still stuck, I don't know with what I have to compare this model with. With random Y-calculations? But then I get a strong biased beta1 estimator.
I have to do an evaluation of the chosen design (G-optimal design) and the point estimators. A cubic relationship is assumed. So what I have done now, is calculating the Y-values with the cubic relationship, the optimal designpoints, and the calculated coefficients (bètas). But indeed, it's not smart to do estimations of beta with this model.. But at the moment, I have no idea of how to evaluate this.
 
#5
We are not so good in understanding Dutch. Your code explanations are cryptic.
(Some of us are not so good in English either but we understand the elementary part.)

"Today I also concluded that this lus is always...."
It seems like "lus" in Dutch means "loop".

Make a fully reproducible example and explain what the problem is and highlight the code and click on the # symbol.