#################### TLNise README FILE ############################
####### Copyright 1999, Phil Everson, Swarthmore College. ###########
#
#This program should be cited if used in published research.
#This is a test version. The results are not guaranteed.
#The author will appreciate notification of any errors by
#   emailing peverson@swarthmore.edu
#

TLNise (Two-Level Normal independent sampling estimation) enables
inference about the parameters {theta1,...,thetaJ,gamma,A} of the 
model:

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

Level-2: thetaj|gamma,A ~ N(Wj'gamma, A), indep., j=1,...,J

####### TLNise.u.v2 HANDLES UNIVARIATE OUTCOMES ONLY ###### 

For version TLNise.u.v2, the outcomes Yj must be univariate, and the 
scalar level-1 variances V={V1,...,VJ} are assumed known, and must be 
input. 

##################### INSTALLATION ################

(1) Compile Fortran code in tlnise.f, creating an object file 
     called tlnise.o.

(2) In the Splus code contained in the file TLNise.u.v2.txt, locate 
     the dynamic load call:

# 			!!!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/Frtrn/Nise/TLNise/tlnise.o")
#          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -> change to agree with your path
#
     Change this path to point to your copy of tlnise.o.

(3) Source the Splus code from within Splus using the command: 

      source("/...path.../TLNise.u.v2.txt") 

     changing "...path..." to point to the directory where TLNise.u.v2.txt
     is stored.

#####################################################

#######  DETAILS AND EXAMPLES  ################

# Inputs to tlnise (from within Splus) are as follows:

tlnise_function(Y,V,W="NA",V0="NA",Vs="NA",P=1,N=1000, alpha=.05,
 Astart="NA",ds="NA",nuo="NA",Tol=1e-6,maxiter=1000,intercept=T, 
 labelW="NA",labelY="NA",digits=3,srt=F,brief=0){
#
# inputs:
#        Y: J-vector of Normal outcomes
#        V: J-vector of Level-1 variances
#        W: Jxr covariate matrix (adds column of 1's if not
#	      included and intercept = T)
#       V0: "typical" Vj (default is average of Vj's)
#       Vs: V* for defining prior on B* = V*/(V*+A)
#        P: prior specification; 
#         P=0: uniform on B0; 
#	  P=1: uniform on A;
#         P=2: uniform on B* = V*/(V*+A) (default V* = V0)
#	  P=3: p(B*) = B*^(ds/2)         (default ds = -2)
#        N: number of truncated Gamma draws for inference
#    alpha: specifies (1-apha)100% confidence level (default=95%)
#   Astart: starting value for locating mode of B0
#	ds: required for prior option P=3
#      nu0: alternate value of df - default is J-r (nu0 <= J-r)
#      Tol: tolerance for determining modal convergence
#  maxiter: maximum number of EM iterations for finding mode
#intercept: if T an intercept term is included in the regression
#   labelW: optional names vector for covariates
#   labelY: optional names vector for observations
#   digits: number of digits beyond decimal for reporting results
#      srt: if True, sorts observations on V
#    brief: level of output, from 0 (minimum) to 2 (maximum)
#
The program output are controlled by the input brief. 
					      =====
brief=0
========
For brief=0 (default), the program outputs two matrices, labeled
"gamma" and "theta". The matrix gamma contains posterior estimates of
the rx1 level-2 regression coefficient for predicting thetaj from the
covariates in wj (of length r, including a possible constant term).  It
includes estimates of the posterior mean and sd and their ratio, as
well as lower and upper alpha/2 cutoffs for the marginal posterior
distributions of each component of theta, where alpha is an input with
default 0.05. Note that the lower and upper quantiles are not exactly
symmetric about the posterior mean estimate because the posterior
distribution is not Normal after integrating out A.

The matrix theta contains posterior mean (theta*) and sd (rtVtheta)
estimates for theta, as well as estimates of mu=Wj'gamma and
B=Vj/(Vj+A). It also produces lower and upper alpha/2 cutoffs of the
posterior distribution for each thetaj.

brief=1
=======
Outputs gamma and theta, as above plus estimated quantiles of the
posterior distribution of A (A), the posterior mode of B0=V0/(V0+A)
(modeB0), the parameters of the enveloping gamma density used to
approximate f(B0|Y) (param=c(2*alpha,2*lambda), and the N approximate
draws made from f(A|Y) (Adraws) and their associated scaled importance
ratios (ratio).

brief=2
=======
Following the brief=1 output are thetam and gammam, matrices of the
posterior means of theta and of gamma, conditional on each sampled
value of A (note that these are JxN and rxN matrices).

If you set brief=3, you will get a glut of output, most likely only
interpretable by me.


Univariate 2-Level Normal Example:
----------------------------------
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
V0_mean(V)  # = 13.2
W_c(1,2,3,4,5,5,4,3,2,1)
mu_5+2*W
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.u.v2:
############################
#
source("/home/everson/Frtrn/Nise/TLNise/TLNise.u.v2.txt")  
# specify path to appropriate directory 
#
out_tlnise(Y,V,W)  
[1] flat prior distribution assumed for A
[1] ********* Locating posterior mode of B0 ************
[1] converged to within tolerance in 23 EM iterations
[1] posterior mode for B0 = 0.382 ( A = 21.32 )

[1] ******** Drawing truncated Gammas *********** 

[1] ************ Processing draws *************** 
[1] Average Importance Weight = 0.946
[1] Estimated quantiles of f(A|Y):
    5%   10%   25% Median    75%    90%    95% 
 1.567 2.648 6.879 14.522 29.034 53.275 72.863
[1] sqrt(A):
    5%   10%   25% Median   75%   90%   95% 
 1.252 1.627 2.623  3.811 5.388 7.299 8.536
[1] Estimated quantiles of f(B0|Y):
[1] B0 =  13.2 /( 13.2 + A ) 
    5%   10%   25% Median   75%   90%   95% 
 0.153 0.199 0.313  0.476 0.657 0.833 0.894
-----
$gamma:
    est    sd ratio 
0 6.346 4.480 1.417
1 1.982 1.333 1.487

$theta:
        Y rt(V)     mu     B theta* rt(Vtheta) 
 1 12.575 5.477  8.328 0.674  9.747      3.834
 2  8.714 5.477 10.310 0.674  9.716      3.473
 3 19.244 3.873 12.292 0.508 15.603      3.254
 4  8.836 3.873 14.274 0.508 11.639      3.085
 5 22.274 3.162 16.256 0.408 19.624      2.981
 6 15.603 3.162 16.256 0.408 15.883      2.632
 7 11.505 2.449 14.274 0.292 12.438      2.157
 8  9.501 2.449 12.292 0.292 10.429      2.122
 9  9.151 2.236 10.310 0.256  9.471      1.952
10  8.482 2.236  8.328 0.256  8.380      2.039



compare squared error loss:
---------------------------

sum((Y-theta)^2/V)/10    		# observed values
[1] 0.7413655
sum((out$theta[,"mu"]-theta)^2/V)/10	# regression means
[1] 0.8507454
sum((out$theta[,"theta"]-theta)^2/V)/10 # shrinkage estimates
[1] 0.4274357