#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)
}