Error in R-code for evaluation of optimal design

Evelien_V

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

Evelien_V

New Member
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.

GretaGarbo

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