I've put two packages for R on my home page at http://www.xs4all.nl/~jhckx/R/. The "pcnt" package is for multiway percentage tables. I've posted a first effort called "ctab" on this group and a request for enhancing "ftable" with percentages on the wishlist. The "mcl" package is for estimating multinomial logistic models using conditional logistic regression. This gives greater flexibility in imposing restrictions on the dependent variable. One application is to estimate loglinear models for sqaure tables (e.g. quasi-independence, quasi-symmetry) with covariates. The mcl package therefore contains a number of functions for models for square tables as well. A caveat is that "clogit" in R doesn't produce the same estimates as "multilog", although the likelihood functions for both models are the same. The maximum absolute difference is 0.0034, the mean absolute difference is 0.00069. Stata's "clogit" and "mlogit" produce the same estimates and match those of "multilog" to at least 6 decimal points accuracy. See the notes in http://www.xs4all.nl/~jhckx/R/mcl.html Can anyone shed any light on this? John Hendrickx
mcl models, percentages
4 messages · John Hendrickx, Thomas Lumley, Göran Broström
On Wed, 14 May 2003, John Hendrickx wrote:
A caveat is that "clogit" in R doesn't produce the same estimates as "multilog", although the likelihood functions for both models are the same. The maximum absolute difference is 0.0034, the mean absolute difference is 0.00069. Stata's "clogit" and "mlogit" produce the same estimates and match those of "multilog" to at least 6 decimal points accuracy. See the notes in http://www.xs4all.nl/~jhckx/R/mcl.html Can anyone shed any light on this?
Two possibilities 1/ not having converged far enough: the convergence tolerance for coxph is by default only 1e-4 (yes, I should change it) 2/ Weights. In one of your examples you have frequency weights passed to clogit. This doesn't work (at least, it isn't equivalent to passing in the expanded data) because weighting doesn't make coxph compute the full set of permutations that you need for the likelihood in a large stratum. -thomas
--- Thomas Lumley <tlumley at u.washington.edu> wrote:
On Wed, 14 May 2003, John Hendrickx wrote:
A caveat is that "clogit" in R doesn't produce the same estimates
as
"multilog", although the likelihood functions for both models are
the
same. The maximum absolute difference is 0.0034, the mean
absolute
difference is 0.00069. Stata's "clogit" and "mlogit" produce the
same
estimates and match those of "multilog" to at least 6 decimal
points
accuracy. See the notes in http://www.xs4all.nl/~jhckx/R/mcl.html
Can
anyone shed any light on this?
Two possibilities 1/ not having converged far enough: the convergence tolerance for coxph is by default only 1e-4 (yes, I should change it) 2/ Weights. In one of your examples you have frequency weights passed to clogit. This doesn't work (at least, it isn't equivalent to passing in the expanded data) because weighting doesn't make coxph compute the full set of permutations that you need for the likelihood in a large stratum.
It appears to be the convergence tolerance. It's not just the weights, the problem occurs for the logan data as well. Specifying "eps=1e-9" in coxph produces the estimates found in other programs. Could you specify a lower default value for eps? It would also be useful to add options in "clogit" for eps and weights. John Hendrickx
On Wed, 14 May 2003, Thomas Lumley wrote:
On Wed, 14 May 2003, John Hendrickx wrote:
A caveat is that "clogit" in R doesn't produce the same estimates as "multilog", although the likelihood functions for both models are the same. The maximum absolute difference is 0.0034, the mean absolute difference is 0.00069. Stata's "clogit" and "mlogit" produce the same estimates and match those of "multilog" to at least 6 decimal points accuracy. See the notes in http://www.xs4all.nl/~jhckx/R/mcl.html Can anyone shed any light on this?
Two possibilities 1/ not having converged far enough: the convergence tolerance for coxph is by default only 1e-4 (yes, I should change it)
I tried to change 'eps' to 1.e-8, found yet another way of crashing R: ------------------------------------------------------------------------
library(survival) coxph(Surv(enter, exit, event) ~ x, data = tt,
control = list(eps = 1.e-8))
Program received signal SIGSEGV, Segmentation fault.
0x405526a0 in agfit3 (maxiter=0x0, nusedx=0x1, nvarx=0x1, start=0x8f267b8,
stop=0x8f26770, event=0x8f25f98, covar2=0x8f26728, offset=0x8f266e0,
weights=0x8f26698, nstrat=0x8ec4ab0, strata=0x8ec4a90,
sort1=0x8f25f60,
sort2=0x8f25f28, means=0x8e7ed30, Rf_beta=0x8e7ed08, u=0x8e7ece0,
imat2=0x8e7ecb8, loglik=0x8f25ef0, flag=0x8ec4a70, work=0x8f249a0,
eps=0x8e7ec90, tol_chol=0x0, sctest=0x8e7ec68) at agfit3.c:263
263 agfit3.c: No such file or directory.
in agfit3.c
-------------------------------------------------------------------------
It's of course my fault; 'control' should be a call to 'coxph.control'.
But should R crash...?
G?ran
[...]
---
G?ran Brostr?m tel: +46 90 786 5223
Department of Statistics fax: +46 90 786 6614
Ume? University http://www.stat.umu.se/egna/gb/
SE-90187 Ume?, Sweden e-mail: gb at stat.umu.se