SO.mc.est {CorrBin}

Order-restricted MLE assuming marginal compatibility

Description

SO.mc.est computes the nonparametric maximum likelihood estimate of the distribution of the number of responses in a cluster P(R=r|n) under a stochastic ordering constraint. Umbrella ordering can be specified using the turn parameter.

Usage

SO.mc.est(cbdata, turn = 1, control = soControl())

Arguments

cbdata

an object of class CBData.

turn

integer specifying the peak of the umbrella ordering (see Details). The default corresponds to a non-decreasing order.

control

an optional list of control settings, usually a call to soControl. See there for the names of the settable control values and their effect.

Details

Two different algorithms: EM and ISDM are implemented. In general, ISDM (the default) should be faster, though its performance depends on the tuning parameter max.directions: values that are too low or too high slow the algorithm down.

SO.mc.est allows extension to an umbrella ordering: D_1 \geq^{st} \cdots \geq^{st} D_k \leq^{st} \cdots \leq^{st} D_n by specifying the value of k as the turn parameter. This is an experimental feature, and at this point none of the other functions can handle umbrella orderings.

Value

A list with components:

Components Q and D are unlikely to be needed by the user.

MLest

data frame with the maximum likelihood estimates of P(R_i=r|n)

Q

numeric matrix; estimated weights for the mixing distribution

D

numeric matrix; directional derivative of the log-likelihood

loglik

the achieved value of the log-likelihood

converge

a 2-element vector with the achieved relative error and the performed number of iterations

Author(s)

Aniko Szabo, [e-mail]

References

Szabo A, George EO. (2010) On the Use of Stochastic Ordering to Test for Trend with Clustered Binary Data. Biometrika 97(1), 95-108.

See Also

soControl

Examples

 data(shelltox)
 ml <- SO.mc.est(shelltox, control=soControl(eps=0.01, method="ISDM"))
 attr(ml, "converge")
 
 require(lattice)
 panel.cumsum <- function(x,y,...){
   x.ord <- order(x)
   panel.xyplot(x[x.ord], cumsum(y[x.ord]), ...)}

 xyplot(Prob~NResp|factor(ClusterSize), groups=Trt, data=ml, type="s",
      panel=panel.superpose, panel.groups=panel.cumsum,
      as.table=TRUE, auto.key=list(columns=4, lines=TRUE, points=FALSE),
      xlab="Number of responses", ylab="Cumulative Probability R(R>=r|N=n)",
      ylim=c(0,1.1), main="Stochastically ordered estimates\n with marginal compatibility")