### Supporting R code for the linear model elicitation example in Daneshkhah and Oakley (2010) ### Implements the method in ### Kadane, J. B., Dickey, J. M., Winkler, R. L., Smith, W. S. and Peters, S. C. (1980). ### Interactive elicitation of opinion for a normal linear model. ### Journal of the American Statistical Association, 75: 845-854. ### Specify design matrix X x1<-rep(c(20000,70000),4) x2<-rep(c(250,250,500,500),2) x3<-c(10,10,10,10,20,20,20,20) x4<-c(600,200,200,600,600,200,200,600) x5<-c(4,8,8,4,8,4,4,8) X<-cbind(rep(1,8),x1,x2,x3,x4,x5) # Can obtain fractional factorial design using the BHH2 library: # library(BHH2) # ffDesMatrix(5,gen=list(c(4,1,2),c(5,3,4))) ### Initial judgements y5<-matrix(c(385,200,575,260,700,390,765,575),8,1) y75<-matrix(c(410,235,605,295,735,430,800,620),8,1) y9<-matrix(c(430,265,630,330,770,470,840,665),8,1) ### Obtain b b<-solve(t(X)%*%X)%*%t(X)%*%y5 ### Compare fitted medians with elicited medians fitted.medians<-X%*%b plot(y5,fitted.medians) ### Choose d a<-(y9-y5)/(y75-y5) #tail ratios. Should inspect for outliers d.fit<-function(delta,mean.a){ (qt(0.9,delta)/qt(0.75,delta)-mean.a)^2 } d<-optimise(d.fit,interval=c(0.1,100),mean.a=mean(a))\$minimum ### hypothetical data y0<-matrix(c(417,259,627,306,775,471),6,1) ### Conditional medians ### Element i,j is median for Y_j conditional on y^0_1,...y^0_i C.C<-matrix(0,6,6) C.C[1,]<-c(390,205,580,265,705,395) C.C[2,]<-c(0,210,585,270,710,400) C.C[3,]<-c(0,0,590,275,715,405) C.C[4,]<-c(0,0,0,280,720,410) C.C[5,]<-c(0,0,0,0,730,415) C.C[6,]<-c(0,0,0,0,0,425) ### Conditional 75th percentiles ### Element i+1 is median for Y_{i+1} conditional on y^0_1,...y^0_i new.75<-c(0,240,615,310,755,455) ### Estimating R and omega S.C<-rep(0,6) U.i<-matrix(((y75[1]-y5[1])/qt(0.75,d))^2,1,1) # U_1 S.C[1]<-U.i[1,1] omega.hat<-rep(0,6) ### First estimate of omega cospread<-(C.C[1,1]-y5[1])*U.i[1,1]/(y0[1]-y5[1]) omega.hat[1]<-U.i[1,1]-cospread for(i in 1:5){ #S(Y_{i+1} | y_1^0,...,y_i^0) S.C[i+1]<-((new.75[i+1]-C.C[i,(i+1)])/qt(0.75,i+d))^2 #L_{i+1} L.i.plus.1<-matrix(0,i,1) for(j in 1:i){ L.i.plus.1[j,1]<-C.C[j,(i+1)]-y5[(i+1),1] } # Matrix used to calulate T_{i+1} M.i.plus.1<-matrix(-y5[1:i],i,i,byrow=T) if(i>1){ for(j in 1:(i-1)){ for(k in (j+1):i){ M.i.plus.1[j,k]<-M.i.plus.1[j,k] + C.C[j,k] } } } for(k in 1:i){ for(j in k:i){ M.i.plus.1[j,k]<- M.i.plus.1[j,k]+y0[k] } } #T_{i+1} T.i.plus.1<-solve(M.i.plus.1)%*%L.i.plus.1 # Element (i+1,i+1) of U_{i+1} S.i.plus.1<-S.C[i+1]*(1+i/d) / (1 + t(y0[1:i]-y5[1:i])%*%solve(U.i)%*%(y0[1:i]-y5[1:i])/d) +t(T.i.plus.1)%*%U.i%*%T.i.plus.1 # Construct U_{i+1} U.i.plus.1<-matrix(0,i+1,i+1) U.i.plus.1[1:i,1:i]<-U.i U.i.plus.1[1:i,(i+1)]<-U.i%*%T.i.plus.1 U.i.plus.1[(i+1),1:i]<-t(U.i.plus.1[1:i,(i+1)]) U.i.plus.1[(i+1),(i+1)]<-S.i.plus.1 U.i<-U.i.plus.1 # Get estimate of omega #F_{i+1} F.i.plus.1<-1/(1+ t(y0[1:i]-X[1:i,]%*%b)%*%solve(U.i.plus.1[1:i,1:i])%*%(y0[1:i]-X[1:i,]%*%b)/d) #H_{i+1} H.i.plus.1<-(C.C[(i+1),(i+1)]-C.C[i,(i+1)])/(y0[i+1]-C.C[i,(i+1)])*S.C[i+1] # Estimate of omega omega.hat[i+1]<-(S.C[i+1] - H.i.plus.1 )*(d+i)/d*F.i.plus.1 } # Smallest eigenvalue of U_6 min.eig<-min(eigen(U.i)\$values) # Either estimate omega from #omega<-mean(omega.hat[(omega.hat