el2.test.wtm: Computes maximum-likelihood probability jumps for multiple...
In emplik2: Empirical Likelihood Ratio Test for Two Samples with Censored Data

Description Usage Arguments Details Value Author(s) References Examples

This function computes the maximum-likelihood probability jumps for multiple mean-type hypotheses, based on two samples that are independent, uncensored, and weighted. The target function for the maximization is the constrained log(empirical likelihood) which can be expressed as,

∑_{dx_i=1} wx_i \log{μ_i} + ∑_{dy_j=1} wy_j \log{ν_j} - η ( 1 -∑_{dx_i=1} μ_i ) - δ ( 1 -∑_{dy_j=1} ν_j ) - λ ( μ^T H_1 ν, … , μ^T H_p ν )^T

where the variables are defined as follows:

x is a vector of uncensored data for the first sample

y is a vector of uncensored data for the second sample

wx is a vector of estimated weights for the first sample

wy is a vector of estimated weights for the second sample

μ is a vector of estimated probability jumps for the first sample

ν is a vector of estimated probability jumps for the second sample

H_k = [ g_k(x_i, y_j) - mean_k ], k=1, … , p, where g_k(x,y) is a user-chosen function

H = [H_1, ... , H_p] (used as argument in el.cen.EMm function, which calls el2.test.wtm)

mean is a vector of length p of hypothesized means, such that mean_k is the hypothesized value of E(g_k(x,y))

E indicates “expected value”

1	el2.test.wtm(xd1,yd1,wxd1new, wyd1new, muvec, nuvec, Hu, Hmu, Hnu, p, mean, maxit=15)

`xd1`	a vector of uncensored data for the first sample
`yd1`	a vector of uncensored data for the second sample
`wxd1new`	a vector of estimated weights for `xd1`
`wyd1new`	a vector of estimated weights for `yd1`
`muvec`	a vector of estimated probability jumps for `xd1`
`nuvec`	a vector of estimated probability jumps for `yd1`
`Hu`	H_u = [ H_1 - [mean_1], … , H_p - [mean_p] ], dx_i=1, dy_j=1
`Hmu`	a matrix, whose calculation is shown in the example below
`Hnu`	a matrix, whose calculation is shown in the example below
`p`	the number of hypotheses
`mean`	a vector of hypothesized values of E(g_k(u,v)), k=1, …,p
`maxit`	a positive integer used to control the maximum number of iterations in the Newton-Raphson algorithm; default is 10

This function is called by el2.cen.EMm. It is listed here because the user may find it useful elsewhere.

The value of mean_k should be chosen between the maximum and minimum values of g_k(xd1_i,yd1_j); otherwise there may be no distributions for xd1 and yd1 that will satisfy the the mean-type hypothesis. If mean_k is inside this interval, but the convergence is still not satisfactory, then the value of mean_k should be moved closer to the NPMLE for E(g(xd1,yd1)). (The NPMLE itself should always be a feasible value for mean_k.) The calculations for this function are derived in Owen (2001).

el2.test.wtm returns a list of values as follows:

`constmat`	a matrix of row vectors, where the kth row vector holds successive values of μ^T H_k ν , k=1, …, p, generated by the Newton-Raphson algorith m
`lam`	the vector of Lagrangian mulipliers used in the calculations
`muvec1`	the vector of probability jumps for the first sample that maximize the weighted empirical likelihood
`nuvec1`	the vector of probability jumps for the second sample that maximize the weighted empirical likelihood
`mean`	the vector of hypothesized means

William H. Barton <bbarton@lexmark.com>

Owen, A.B. (2001). Empirical Likelihood. Chapman and Hall/CRC, Boca Raton, pp.223-227.

 
#Ho1: P(X>Y) = 0.5
#Ho2: P(X>1060) = 0.5
#g1(x) = I[x > y]
#g2(y) = I[x > 1060]

mean<-c(0.5,0.5)
p<-2

xd1<-c(10,85,209,273,279,324,391,566,852,881,895,954,1101,1393,1669,1823,1941)
nx1=length(xd1)
yd1<-c(21,38,39,51,77,185,524,610,612,677,798,899,946,1010,1074,1147,1154,1329,1484,1602,1952)
ny1=length(yd1)

wxd1new<-c(2.267983, 1.123600, 1.121683, 1.121683, 1.121683, 1.121683, 1.121683,
 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.261740, 2.912753,
 2.912753, 2.912753)
muvec<-c(0.08835785, 0.04075290, 0.04012084, 0.04012084, 0.04012084, 0.04012084,
 0.04012084, 0.03538020, 0.03389263, 0.03389263, 0.03389263, 0.03322693, 0.04901516,
 0.05902008, 0.13065491, 0.13065491, 0.13065491)

wyd1new<-c(1.431653, 1.431653, 1.431653, 1.431653, 1.431653, 1.438453, 1.079955, 1.080832,
 1.080832, 1.080832, 1.080832, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000,
 1.222883, 1.227865, 1.739636, 5.809616)
nuvec<-c(0.04249966, 0.04249966, 0.04249966, 0.04249966, 0.04249966, 0.04316922, 0.03425722,
 0.03463312, 0.03463312, 0.03463312, 0.03463312, 0.03300598, 0.03300598, 0.03333333,
 0.03333333, 0.03382827, 0.03382827, 0.04136800, 0.04229270, 0.05992020, 0.22762676)


H1u<-matrix(NA,nrow=nx1,ncol=ny1)
H2u<-matrix(NA,nrow=nx1,ncol=ny1)
for (i in 1:nx1) {
   for (j in 1:ny1) {
        H1u[i,j]<-(xd1[i]>yd1[j])
        H2u[i,j]<-(xd1[i]>1060) } }
Hu=matrix(c(H1u,H2u),nrow=nx1,ncol=p*ny1)
for (k in 1:p) {
     M1 <- matrix(mean[k], nrow=nx1, ncol=ny1)
     Hu[,((k-1)*ny1+1):(k*ny1)] <- Hu[,((k-1)*ny1+1):(k*ny1)] - M1}
Hmu <- matrix(NA,nrow=p, ncol=ny1*nx1)
Hnu <- matrix(NA,nrow=p, ncol=ny1*nx1) 
for (i in 1:p) {
   for (k in 1:nx1) {
        Hmu[i, ((k-1)*ny1+1):(k*ny1)] <-  Hu[k,((i-1)*ny1+1):(i*ny1)] } }
for (i in 1:p) {
   for (k in 1:ny1) {
        Hnu[i,((k-1)*nx1+1):(k*nx1)] <- Hu[(1:nx1),(i-1)*ny1+k]} }

el2.test.wtm(xd1,yd1,wxd1new, wyd1new, muvec, nuvec, Hu, Hmu,
  Hnu, p, mean, maxit=10)

#muvec1
# [1] 0.08835789 0.04075290 0.04012083 0.04012083 0.04012083 0.04012083 0.04012083
# [8] 0.03538021 0.03389264 0.03389264 0.03389264 0.03322693 0.04901513 0.05902002
# [15] 0.13065495 0.13065495 0.13065495

#nuvec1
# [1] 0.04249967 0.04249967 0.04249967 0.04249967 0.04249967 0.04316920 0.03425722
# [8] 0.03463310 0.03463310 0.03463310 0.03463310 0.03300597 0.03300597 0.03333333
# [15] 0.03333333 0.03382827 0.03382827 0.04136801 0.04229269 0.05992018 0.22762677

#  $lam
#        [,1]      [,2]
#  [1,] 8.9549 -10.29119