#Copyright Sept 1993, Phil Everson and Carl Morris, #Dept. of Statistics, Harvard University. #Last Update 8/31/94 #NOTE: The previous version of this program gave incorrect #estimates for beta in the case when Tau^2 was estimated to #be zero. In this case the beta estimates should agree with #Weighted Least Squares results with weights equal to 1/Vi. #This program should be cited if used in published research. #This is a test version. The results are not guaranteed. #The authors will appreciate notification of any errors by # emailing everson@stat.harvard.edu #Permission from the authors is necessary for further distribution # and use of this program. We will notify users of new versions. #Reference: # Parametric Empirical Bayes Inference: Theory and Applications # Carl N. Morris, JASA, March 1983 #estimate parameters for normal/normal hierarchical model. #Allows matrix of covariates, X, for estimating prior mean. #Model: # Y.i|theta.i ~ N(theta.i,V.i) independently (V.i assumed known) # theta.i ~ N(X.i*beta,Tau^2) independently #Inputs: Y = vector of observed values (unbiased estimates for # true parameters) # V = vector of variances for each observation # X = matrix of covariates (may include column of 1's # if a constant term is to be included in the model) #Intercept = T if intercept to be included in model. A column # of 1's will be added to X if it wasn't included. # title = title for output # labelx = labels for the regression coefficients # labelobs = labels for observations # srt = if true, will sort output on V # tol = tolorance level for convergence # dig = # digits for output #Output: Convergence Iterations for Tau^2 # Summary: summary statistics for Y, rtV and X # Hyperparameter: beta^ estimates, standard errors and ratio # Tau2: estimate of prior variance and standard error # Tau: sqrt(Tau^2) # aveB: average of estimated shrinkage factors # output: matrix of # Y = observed outcomes, # sqrtV = standard deviations # mu = estimated prior means, # B = estimated shrinkage factors # se(B) = estimated standard errors for B # theta^ = estimated posterior means # se(thet^) = estimated standard errors for theta^ norm.hm=function(Y,V,X=NA,intercept=T,title="",labelx=NULL, labelobs=NULL,srt=F,tol=1e-5,dig=3){ k=length(Y) rtV=sqrt(V) if(missing(X)){ X=matrix(rep(1,k),ncol=1) r=1} else{ X=as.matrix(X) if(intercept){ if(max(abs(X[,1]-rep(1,k)))!=0) X=cbind(rep(1,k),X) } r=dim(as.matrix(X))[[2]] } if(missing(labelx)){ if(intercept) labelx=c(0:(r-1)) else labelx=c(1:r) } if(missing(labelobs)) labelobs=1:k dimnames(X)=list(labelobs,labelx) dat=cbind(Y,rtV,X) summary.mat=apply(dat,2,summary) summary.mat=apply(summary.mat*10^dig,2,round)/10^dig tau2=mean(V) #initial guess at tau2 Tol=1 print(paste("**********", title, "**********"),quote=F) print(" Hierarchical Normal Regression Model",quote=F) print(" copyright Sept. 1993 P. J. Everson & C. N. Morris",quote=F) cat("\n\n") cat("**********CONVERGENCE ITERATIONS FOR TAU^2**********", "\n\n") while(Tol>tol){ #iterate until within tolorance W=1/(V+tau2) D=diag(W) var.beta=solve(t(X)%*%D%*%X) beta=var.beta%*%t(X)%*%D%*%Y mu=X%*%beta tau22=max((sum(W*((k/(k-r))*(Y-mu)^2-V))/sum(W)),0) Vtau22=2/((k-r)*mean(W)^2) Tol=abs(tau2-tau22)/sqrt(Vtau22) tau2=tau22 cat(signif(tau2,5),",") } W=1/(V+tau2) D=diag(W) var.beta=solve(t(X)%*%D%*%X) beta=var.beta%*%t(X)%*%D%*%Y mu=X%*%beta Vtau2=2/((k-r)*mean(W)^2) cat("\n\n") sd.beta=sqrt(diag(var.beta)) ratio=round(beta/sd.beta,dig) sd.beta=round(sd.beta,dig) beta=round(beta,dig) beta.mat=cbind(beta,sd.beta,ratio) tau=round(sqrt(tau2),dig) B=((k-r-2)/(k-r))*V/(V+tau2) rr=k*W*diag(X%*%solve(t(X)%*%D%*%X)%*%t(X)) Vbar=sum(W*V)/sum(W) v=(2/(k-r-2))*B^2*(Vbar+tau2)/(V+tau2) rtv=sqrt(v) s2=V*(1-((k-rr)/k)*B)+v*(Y-mu)^2 theta=(1-B)*Y+B*mu s=sqrt(s2) aveB=mean(B) dimnames(beta.mat)=list(paste("b^",labelx),c("beta^","se","beta^/se")) tau2.vec=c(tau2,signif(sqrt(Vtau2),dig-1)) names(tau2.vec)=c("Tau^2","SE(Tau^2)") out=cbind(Y,sqrt(V),mu,B,rtv,theta,s) out=apply((10^dig)*out,2,round)/10^dig dimnames(out)=list(labelobs,c("Y","sqrtV","mu","B","se(B)", "theta^", "se(thet^)")) if(srt) out=out[order(V),] list(summary=summary.mat,hyperparameter=beta.mat,Tau2=tau2.vec, Tau=tau,aveB=round(aveB,dig),output=out) }