## When using this code, please cite this article: de Bekker-Grob, E.W., Donkers, B., Jonker, M.F. et al. Patient (2015) 8: 373. https://doi.org/10.1007/s40271-015-0118-z##
############################################################################################################################################################################
############################################################################################################################################################################


test_alpha=0.05
z_one_minus_alpha<-qnorm(1-test_alpha)


test_beta=0.2
z_one_minus_beta<-qnorm(1-test_beta)


parameters<-c(1.23 , -0.31 , -0.21 , -0.44 , 0.028 , -1.10 , -0.04 , -0.0015)


ncoefficients=8
nalts=3
nchoices=16


# load the design information
design<-as.matrix(read.table("D:\\Design_matrix_Illustration_DCE_Osteoporosis_treatment.txt",header=FALSE));


#compute the information matrix
# initialize a matrix of size ncoefficients by ncoefficients filled with zeros.
info_mat=matrix(rep(0,ncoefficients*ncoefficients), ncoefficients, ncoefficients) 
# compute exp(design matrix times initial parameter values) 
exputilities=exp(design%*%parameters)
# loop over all choice sets
for (k_set in 1:nchoices) {
# select alternatives in the choice set
alternatives=((k_set-1)*nalts+1) : (k_set*nalts)
# obtain vector of choice shares within the choice set
p_set=exputilities[alternatives]/sum(exputilities[alternatives])
# also put these probabilities on the diagonal of a matrix that only contains zeros
p_diag=diag(p_set)
# compute middle term P-pp’
middle_term<-p_diag-p_set%o%p_set
# pre- and postmultiply with the Xs from the design matrix for the alternatives in this choice set
full_term<-t(design[alternatives,])%*%middle_term%*%design[alternatives,]
# Add contribution of this choice set to the information matrix
info_mat<-info_mat+full_term 
} # end of loop over choice sets
#get the inverse of the information matrix (i.e., gets the variance-covariance matrix)
sigma_beta<-solve(info_mat,diag(ncoefficients)) 


# Use the parameter values as effect size. Other values can be used here.
effectsize<-parameters
# formula for sample size calculation is n>[(z_(beta)+z_(1-alpha))*sqrt(S??)/delta]^2 
N<-((z_one_minus_beta + z_one_minus_alpha)*sqrt(diag(sigma_beta))/abs(effectsize))^2
# Display results (required sample size for each coefficient)
N