##################### TLNise.mv1 README FILE  #######################
####### Copyright 2000, Phil Everson, Swarthmore College. ###########

TLNise stands for "Two-Level Normal independent sampling estimation". 
This program enables inferences for the following two-level model:

Level-1:  Yj|thetaj,Vj ~ Np(thetaj, Vj), j=1,..,J

Level-2: thetaj|gamma,A ~ Np(Wjgamma, A)

With tlnise.mv1, responses Yj may be scalars or p-dimensional vectors,
in which case the Vj's are pxp covariance matrices. See Everson and
Morris (2000) JRSSB 62 (399-412) for details. The Vj's are assumed 
known for this analysis. The next release of TLNise will allow for 
unknown level-1 variances in the univariate case.

INSTALLATION:

1) Save tlnisemv1.f with that name, and compile the Fortran code in 
   your environment; e.g., type "f77 -c tlnisemv.f"  (or "g77") to
   create a Fortran object file (tlnise.o).

2) Edit TLNise.mv1.src to point to your copy of tlnise.o:

#
# 			!!!IMPORTANT!!! 
# you must change the path in the following line to specify the directory in
# which you stored the Fortran object file tlnise.o.
 dyn.load("/home/everson/TLNise/tlnisemv1.o")
#          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -> change to agree with your path
#

3) Source TLNise.mv1.src into Splus ( "source("/..path../TLNise.mv1.src") to
   load the Splus functions.


EXECUTION:

The minimal inputs to tlnise are Y and V. If p=1 (univariate) then both
are vectors. If p > 1 (multivariate) then Y is a Jxp (or pxJ) matrix, and
V is a (pxpxJ) array of level-1 covariance matrices.


# INPUTS:
#        Y: Jxp (or pxJ) matrix of p-dimensional Normal outcomes
#        V: pxpxJ array of pxp Level-1 covariances (assumed known)
#        w: Jxq (or qxJ) covariate matrix (adds column of 1's if 
#	     not included and intercept = T)
#       V0: "typical" Vj (default is average of Vj's)
#    prior: Prior parameter: p(B0) = |B0|^{(prior-p-1)/2}
#      prior = -(p+1): uniform on level-2 covariance matrix A (default) 
#      prior =      0: Jeffreys' prior
#      prior =  (p+1): uniform on shrinkage matrix B0
#        N: number of Constrained Wishart draws for inference
#      Tol: tolerance for determining modal convergence
#  maxiter: maximum number of EM iterations for finding mode
#intercept: if True, an intercept term is included in the regression
#   labelw: optional names vector for covariates
#   labelY: optional names vector for the J observations
#  labelYj: optional names vector for the p elements of Yj
#   digits: number of significant digits for reporting results
#    brief: level of output, from 0 (minimum) to 2 (maximum)
#     prnt: Controls printing during execution
#
# OUTPUTS:
#  *** brief=0 ***
#    gamma:  matrix of posterior mean and sd estimates of Gamma,
#             and their ratios.
#    theta:  pxJ matrix of posterior mean estimates for thetaj's.
#  SDtheta:  pxJ matrix of posterior SD estimates for thetaj's.
#        A:  pxp estimated posterior mean of variance matrix A.
#      rtA:  p-vector of between group SD estimates.
#
#  *** brief=1 *** 
#   brief=0 outputs, plus
#  Dgamma: rxr estimated posterior covariance matrix for Gamma.
#  Vtheta: pxpxJ array of estimated covariances for thetaj's.
#      B0: pxpxN array of simulated B0 values.
#      lr: N-vector of log density ratios for each B0 value.
#  
#  *** brief=2 *** 
#   brief=1 outputs, plus
#      lf: N-vector of log f(B0|Y) evaluated at each B0.
#     lf0: N-vector of log f0(B0|Y) evaluated at each B0
#           (f0 is the CWish envelope density for f).
#      df: degrees of freedom for f0.
#   Sigma: scale matrix for f0.
#    nvec: number of matrices begun, diagonal and off-diagonal
#           elements simulated to get N CWish matrices.
#    nrej: number of rejections that occurred at each step 1,..,p.
###########



EXAMPLES:

1) UNIVARIATE:

J_10
n_c(5,5,10,10,15,15,25,25,30,30)
V_150/n
# eg, Yj's are averages of nj N(thetaj,150) observations, so Vj=150/nj
w_c(1,2,3,4,5,5,4,3,2,1)
mu_5+2*w    # gamma = (5,2)'
A_10        # (level-2 standard deviation: sqrt(A) = 3.162278)
set.seed(5)
theta_sqrt(A)*rnorm(10)+mu
#> theta   (simulated theta values)
# [1]  3.036296  6.681080 14.467340  9.643961 20.387069 19.355127 13.179958
# [8]  9.668313  7.811781  9.129335
Y_sqrt(V)*rnorm(10)+theta
#> Y	   (simulated Y values)
# [1] 12.574769  8.714413 19.244139  8.835750 22.274406 15.603354 11.504766
# [8]  9.500900  9.150601  8.482405
#

## Inference via TLNise.mv1:

source("/home/everson/TLNise/TLNise.mv1.src")
# (specify path to appropriate directory)

out_tlnise(Y,V,w,seed=123456789)


[1] ******** Prior Parameter = -2

[1] ******** Locating Posterior Mode ************ 

[1] Converged in 50 EM iterations.
[1] Posterior mode of B0:

      [,1] 
[1,] 0.382

[1] lf(modeB0) = -4.054 ; lf0(modeB0) = -3.924 ; adj = 0.1295

[1] ******** Drawing Constrained Wisharts ******** 


[1] CWish acceptance rate = 1000 / 1000 = 1

[1] ******** Processing Draws ************** 

[1] Average scaled importance weight = 0.9415

[1] Posterior mean estimate of A:
      [,1] 
[1,] 25.34

[1] Between-group SD estimate:
      [,1] 
[1,] 5.034


## Parameter inference:

out$gamma

       est       se   est/se 
0 6.361285 4.550758 1.397852
1 1.980062 1.354638 1.461690

out$theta
        1        2        3        4        5        6       7        8 
 9.796958 9.710735 15.69414 11.57806 19.70825 15.87692 12.4019 10.39493
       9       10 
 9.46121 8.391604

out$SDtheta
        1        2        3        4        5        6        7        8 
 3.863312 3.512786 3.248165 3.096908 2.966085 2.649943 2.164782 2.132232
        9       10 
 1.964715 2.046164



2) MULTIVARIATE:

J_20
p_3
# Set an A matrix to generate data.
A_1.5*diag(3)
A[2,1]_A[1,2]_0.3
A[3,2]_A[2,3]_0.2
rtA_t(chol(A))
# generate thetas ~ N3(0, A)  (no covariates):
theta_matrix(0,ncol=J,nrow=p)
set.seed(5555555)
for(j in 1:J){
 thetaj_matrix(rnorm(p))
 theta[,j]_rtA%*%thetaj
 }
nj_c(rep(10,5),rep(20,5),rep(30,5),rep(40,5))
Y_matrix(0,ncol=J,nrow=p)
V_array(0,c(3,3,J))
# Generate p=3 dimensional N(thetaj,50*I) samples of size nj, j=1,..,J.
V0_0*A
for(j in 1:J){
 Yj_matrix(rnorm(nj[j]*p,0,sqrt(50)),ncol=nj[j],nrow=p)+theta[,j]
 Y[,j]_apply(Yj,1,mean)
 Zj_Yj-Y[,j]
 V0_V0+Zj%*%t(Zj)
 }
V0_V0/sum(nj) # estimate of within covariance (=50*I)

V0

#[1,] 50.7150198   1.741149   0.5318248
#[2,]  1.7411489  54.992984  -1.1571739
#[3,]  0.5318248  -1.157174  46.8230153

# Compute Vj = V0/nj, j=1,..,J:

for(j in 1:J){
 V[,,j]_V0/nj[j]
 }

# Inference via tlnise:

out_tlnise(Y,V,seed=987654321)
[1] ******** Prior Parameter = -4

[1] ******** Locating Posterior Mode ************ 

[1] Converged in 54 EM iterations.
[1] Posterior mode of B0:

         [,1]     [,2]    [,3] 
[1,] 0.505273 0.006369 0.10653
[2,] 0.006369 0.393624 0.08319
[3,] 0.106530 0.083195 0.44365

[1] lf(modeB0) = -32.54 ; lf0(modeB0) = -30.37 ; adj = 2.171

[1] ******** Drawing Constrained Wisharts ******** 


[1] CWish acceptance rate = 1000 / 1304 = 0.7669

[1] ******** Processing Draws ************** 

[1] Average scaled importance weight = 0.5776

[1] Posterior mean estimate of A:
       [,1]   [,2]   [,3] 
[1,]  3.542  0.378 -1.205
[2,]  0.378  5.227 -1.398
[3,] -1.205 -1.398  4.270

[1] Between-group SD estimate:
[1] 1.882 2.286 2.066

# estimates of the means of each thetaij (p=3):
# (true means are all 0).
out$gamma
        est        se    est/se 
0 0.5185656 0.5400422 0.9602316
1 0.2885743 0.6231494 0.4630901
2 0.5374147 0.5664382 0.9487613


# estimates of 1st 5 thetas:

out$theta[,1:5]
            1          2          3          4          5 
1 -0.25071812  0.6075118 -0.2415335  1.3768032 -0.3271834
2 -0.07639761 -0.7983987  0.1679034  0.8021381  1.3590052
3  1.47776465  1.3437457  1.2525367 -1.2017726  0.5979628

out$SDtheta[,1:5]
         1        2        3        4        5 
1 1.362928 1.353548 1.358573 1.398708 1.397352
2 1.556334 1.569360 1.551943 1.589873 1.601760
3 1.417908 1.415371 1.406529 1.495524 1.424882


# compare to true thetas:
theta[,1:5]
           [,1]       [,2]       [,3]       [,4]       [,5] 
[1,] -1.1305050 -1.0780958 -1.3360882 -0.3071049 -0.2415728
[2,] -1.9076087 -0.1059074 -0.8573725 -0.5588133  0.4234559
[3,]  0.5093163  0.9792807 -0.4555741 -1.0687344 -0.5667646

summary((out$theta-theta)/out$SDtheta)

   Min. 1st Qu. Median   Mean 3rd Qu.  Max. 
 -1.358 -0.2756  0.209 0.2776  0.8766 1.651

# deviations are all small in terms of standard deviations.