# vcbd.r
#
# Variance Components for Triply-Nested Binary Responses
#
# Call: vcbd(i, m, y)
#   i : ID of practice (top-level unit)
#   m : Number of patients seen by a provider (2nd level unit)
#   y : Number of responses
#
# Each provider (2nd level) appears as one row of data
#
#
# Model:
#
# At the top (first) level alpha_i has mean theta and variance
#                       theta*(1-theta)*lambda_1.
# Given alpha_i, at the second level, beta_ij has mean alpha_i
# and variance
#                    alpha_i*(1-alpha_i)*lambda_2.
# Given alpha_i and beta_ij, at the third level, the responses Y_ijk are
# independent Bernoulli with mean beta_ij.
#
# Estimation:
#
# Compute sums of squares for a nested study of subjects within providers
# within practices.  Equate expected and observed sums of squares and
# solve for theta, lambda_1 and lambda_2.
#
# Output:
#
# Returns a list containing parameter estimates and the following derived
# quantities: a breakdown of the total variance into the three components
# given as percentages, the two correlations and two odds ratios.
# The list contains also:
# Number of top level units   k
# Number of 2nd level units   n
# Number of 3rd level units   m
# Mean response               ybar
#
# The auxiliary function vcbd.print generates a friendly printout.
#
# Example:
#
# source('vcbd.r')
#
# i = c ( 1, 1, 2, 3, 3, 3, 4, 4, 5, 5, 5, 5 ) # practice ID
# m = c ( 4, 4, 4, 3, 3, 4, 4, 3, 6, 4, 5, 3 ) # patients
# y = c ( 1, 3, 3, 2, 3, 3, 4, 3, 1, 4, 1, 0 ) # patients with response=1
# vc = vcbd(i, m, y)
# vcbd.print(vc)
# vc = vcbd(i, m, y, print_flag = T)  # prints right away
#
#************************************************************************;

vcbd = function (i, mij, yij, print_flag=F) {

smry = function (x, g, f) {
  # summarize variable x by group g using function f
  u = sort(unique(g))
  k = length(u)
  y = double(k)
  for (i in 1:k) y[i] = f(x[g == u[i]])
  list (group=u, smry=y)
}

ssq = function(x) {sum(x^2)}

corr2or = function(u, r) {
  # odds ratio (Y1,Y2) if E[Y1]=E[Y2]=u, corr(Y1,Y2)=r
  u12 = u*u + r*u*(1-u)
  odds_ratio = u12*(1-u-u+u12)/(u-u12)^2
  return(odds_ratio)
}

  a  = sum(yij ^ 2 / mij)
  yi = smry(yij,  i, sum)$smry
  y  = sum(yi)
  mi = smry(mij,  i, sum)$smry
  m  = sum(mi)
  h  = sum(mi^2)
  b  = sum(yi ^ 2  / mi)
  qi = smry(mij,  i, ssq)$smry
  q  = sum(qi)
  r  = sum(qi/mi)
  n  = length(yij)
  k  = length(yi)

  ss1 = y - a;
  ss2 = a - b;
  ss3 = b - y ^ 2 / m;
  x_ = ss1/(m-n);
  y_ = (ss2 - (n-m-k+r)*x_)/(m-r);
  z_ = (ss3 - (k-1-r+q/m)*x_ - (r-q/m-m+h/m)*y_)/(m-h/m)

  ybar = y/m
  if (4*z_ > 1)   th = ybar
  else {
    th1 = (1 + sqrt(1-4*z_)) / 2
    th2 = (1 - sqrt(1-4*z_)) / 2
    if (abs(ybar-th1) < abs(ybar-th2))   th = th1
    else th = th2       # theta is the value closest to ybar
  }

  la1 = 1 - y_ / z_ #  lambda1
  la2 = 1 - x_ / y_ #  lambda2
  if (la1 < 0)  la1 = 0 else if (la1 > 1) la1 = 1
  if (la2 < 0)  la2 = 0 else if (la2 > 1) la2 = 1

  comp1 = la1 * 100
  comp2 = (1-la1)*la2 * 100
  comp3 = (1-la1)*(1-la2) * 100
  corr1 = comp1 / 100
  corr2 = (comp1 + comp2) / 100
  or1 = corr2or(th, corr1)
  or2 = corr2or(th, corr2)

  vc =
  list(k = k, n = n, m = m, ybar = ybar,
       theta = th, lambda1 = la1, lambda2 = la2,
       comp1 = comp1, comp2 = comp2, comp3 = comp3,
       corr1 = corr1, corr2 = corr2, or1 = or1, or2 = or2 )

  if (print_flag)  vcbd.print(vc)

  return(vc)
}

vcbd.print = function (vc) {

  print("Variance Components for Triply-Nested Binary Responses (Version 1.0)", quote = F)
  print("(c) Bahjat Qaqish (2005)", quote = F)
  print("Department of Biostatistics, CB 7420", quote = F)
  print("Chapel Hill,  NC 27599-7420", quote = F)
  print("", quote = F)

  print(paste("Number of practices:  ", vc$k), quote = F)
  print(paste("Number of physicians: ", vc$n), quote = F)
  print(paste("Number of patients:   ", vc$m), quote = F)
  print(paste("Mean response         ", vc$ybar), quote = F)

  print("Estimates:", quote = F)
  print(paste("    Theta:    ", vc$theta ), quote = F)
  print(paste("    Lambda_1: ", vc$lambda1), quote = F)
  print(paste("    Lambda_2: ", vc$lambda2), quote = F)

  print("Variance components:", quote = F)
  print(paste("    Between practices:                 ", vc$comp1, "%"), quote = F)
  print(paste("    Between providers within practice: ", vc$comp2, "%"), quote = F)
  print(paste("    Between patients  within provider: ", vc$comp3, "%"), quote = F)

  print("Correlation:", quote = F)
  print(paste("    Between providers within practice: ", vc$corr1), quote = F)
  print(paste("    Between patients  within provider: ", vc$corr2), quote = F)

  print("Odds ratio:", quote = F)
  print(paste("    Between providers within practice: ", vc$or1), quote = F)
  print(paste("    Between patients  within provider: ", vc$or2), quote = F)
}