Check results for 'Design'

using R version 2.9.1 Patched (2009-07-02 r48890)
using session charset: UTF-8
checking for file 'Design/DESCRIPTION' ... OK
this is package 'Design' version '2.2-0'
checking package dependencies ... OK
checking if this is a source package ... OK
checking whether package 'Design' can be installed ... OK
checking package directory ... OK
checking for portable file names ... OK
checking for sufficient/correct file permissions ... OK
checking DESCRIPTION meta-information ... OK
checking top-level files ... OK
checking index information ... OK
checking package subdirectories ... OK
checking R files for non-ASCII characters ... OK
checking R files for syntax errors ... OK
checking whether the package can be loaded ... OK
checking whether the package can be loaded with stated dependencies ... OK
checking for unstated dependencies in R code ... OK
checking S3 generic/method consistency ... OK
checking replacement functions ... OK
checking foreign function calls ... WARNING
Foreign function calls without 'PACKAGE' argument:
.Fortran("avia", ...)
.Fortran("avia", ...)
.Fortran("lrmfit", ...)
.Fortran("lrmfit", ...)
.Fortran("lrmfit", ...)
.Fortran("lrmfit", ...)
.Fortran("matinv", ...)
.Fortran("robcovf", ...)
.C("agsurv1", ...)
.C("agsurv2", ...)
.C("S_agsurv2", ...)
.C("agsurv2", ...)
See the chapter 'System and foreign language interfaces' of the 'Writing R
Extensions' manual.
checking R code for possible problems ... NOTE
Design: possible error in eval(as.name(XDATADIST), local = FALSE):
unused argument(s) (local = FALSE)
Getlim: possible error in eval(as.name(XDATADIST), local = FALSE):
unused argument(s) (local = FALSE)
Mean.cph: multiple local function definitions for ‘f’ with different
formal arguments
bj: no visible binding for global variable ‘glm.links’
bootcov: no visible global function definition for ‘origGlmFamily’
calibrate.cph: no visible binding for global variable ‘pred.obs’
calibrate.default: no visible binding for global variable ‘.orig.cal’
calibrate.psm: no visible binding for global variable
‘survReg.distributions’
calibrate.psm: no visible binding for global variable ‘glm.links’
calibrate.psm: no visible binding for global variable ‘pred.obs’
cph: possible error in eval(yy, data = data): unused argument(s) (data
= data)
datadist: possible error in eval(as.name(nm), local = FALSE): unused
argument(s) (local = FALSE)
datadist: no visible global function definition for ‘database.object’
gendata.Design: no visible global function definition for ‘using.X’
gendata.Design: no visible global function definition for ‘ed’
gendata.Design: possible error in assign(".df.sub.", df.sub, where =
1): unused argument(s) (where = 1)
gendata.Design: no visible global function definition for ‘Edit.data’
gendata.Design: no visible binding for global variable ‘.df.sub.’
gendata.Design: possible error in get(".df.sub.", where = 1): unused
argument(s) (where = 1)
gendata.Design: no visible global function definition for ‘data.ed’
gendata.default: no visible binding for global variable ‘obj’
glsD: no visible binding for global variable ‘aNlm’
glsD: no visible global function definition for ‘glsApVar’
latex.pphsm: no visible binding for global variable ‘f’
lrm: no visible binding for global variable ‘response’
lrm: possible error in lrm.fit.strat(X, Y, Strata, offset = offs,
penalty.matrix = penalty.matrix, strata.penalty = strata.penalty, tol
= tol, weights = weights, normwt = normwt, ...): unused argument(s)
(weights = weights, normwt = normwt)
lrm: possible error in lrm.fit.strat(X, Y, Strata, penalty.matrix =
penalty.matrix, strata.penalty = strata.penalty, tol = tol, weights =
weights, normwt = normwt, ...): unused argument(s) (weights =
weights, normwt = normwt)
oldDesignFit2R: no visible binding for global variable
‘linear.predictors’
oldDesignFit2R: no visible binding for global variable ‘stats’
oldDesignFit2R: no visible binding for global variable ‘freq’
ols: no visible binding for global variable ‘response’
ols.influence: no visible global function definition for ‘left.solve’
pentrace: no visible global function definition for ‘nlminb.control’
print.lrm: no visible binding for global variable ‘est.exp’
print.psm: no visible binding for global variable
‘survReg.distributions’
print.psm: no visible global function definition for ‘summary.survReg’
print.psm: no visible global function definition for
‘print.summary.survreg’
psm: no visible global function definition for ‘survReg.control’
psm: no visible binding for global variable ‘survReg.distributions’
psm: no visible global function definition for ‘survpenal.fit’
psm: no visible global function definition for ‘survReg.fit’
residuals.psm: no visible global function definition for
‘residuals.survReg’
residuals.psm: no visible global function definition for
‘residuals.survreg’
survest.psm: no visible binding for global variable
‘survReg.distributions’
survest.psm: no visible binding for global variable ‘glm.links’
survplot.Design: no visible global function definition for ‘oldCut’
survreg.fit2: no visible binding for global variable
‘survReg.distributions’
survreg.fit2: no visible global function definition for ‘survReg.fit’
survreg.fit2: no visible global function definition for
‘survReg.control’
val.prob: possible error in legend(lp, leg, lty = lt, marks = marks,
mkh = mkh, cex = cex, bty = "n"): unused argument(s) (marks = marks,
mkh = mkh)
val.prob: no visible global function definition for ‘oldCut’
val.surv: no visible binding for global variable
‘survReg.distributions’
val.surv: no visible binding for global variable ‘glm.links’
validate.psm: no visible binding for global variable ‘inverse’
validate.tree: no visible binding for global variable ‘prune.tree’
validate.tree: no visible binding for global variable ‘response’
validate.tree: no visible global function definition for ‘tree’
checking Rd files ... OK
checking Rd files against version 2 parser ... OK
checking Rd cross-references ... OK
checking for missing documentation entries ... OK
checking for code/documentation mismatches ... OK
checking Rd \usage sections ... OK
checking line endings in C/C++/Fortran sources/headers ... OK
checking line endings in Makefiles ... OK
checking for portable use of $BLAS_LIBS ... OK
checking examples ... ERROR
Running examples in 'Design-Ex.R' failed.
The error most likely occurred in:

> ### * zzzDesignOverview
>
> flush(stderr()); flush(stdout())
>
> ### Name: DesignOverview
> ### Title: Overview of Design Library
> ### Aliases: DesignOverview Design.Overview
> ### Keywords: models
>
> ### ** Examples
>
> ######################
> # Detailed Example 1 #
> ######################
> # May want to first invoke the Hmisc store function
> # so that new variables will go into a temporary directory
> set.seed(17) # So can repeat random number sequence
> n <- 500
>
> sex <- factor(sample(c('female','male'), n, rep=TRUE))
> age <- rnorm(n, 50, 10)
> sys.bp <- rnorm(n, 120, 7)
>
> # Use two population models, one with a systolic
> # blood pressure effect and one without
>
> L <- ifelse(sex=='female', .1*(pmin(age,50)-50), .005*(age-50)^2)
> L.bp <- L + .4*(pmax(sys.bp,120)-120)
>
> dz <- ifelse(runif(n) <= plogis(L), 1, 0)
> dz.bp <- ifelse(runif(n) <= plogis(L.bp), 1, 0)
>
> # Use summary.formula in the Hmisc library to summarize the
> # data one predictor at a time
>
> s <- summary(dz.bp ~ age + sex + sys.bp)
> options(digits=3)
> print(s)
dz.bp N=500

+-------+-----------+---+-----+
| | |N |dz.bp|
+-------+-----------+---+-----+
|age |[15.1,43.6)|125|0.648|
| |[43.6,50.3)|125|0.688|
| |[50.3,56.7)|125|0.616|
| |[56.7,89.4]|125|0.760|
+-------+-----------+---+-----+
|sex |female |240|0.596|
| |male |260|0.754|
+-------+-----------+---+-----+
|sys.bp |[102,115) |125|0.512|
| |[115,120) |125|0.480|
| |[120,125) |125|0.744|
| |[125,141] |125|0.976|
+-------+-----------+---+-----+
|Overall| |500|0.678|
+-------+-----------+---+-----+
> plot(s)
>
> plsmo(age, dz, group=sex, fun=qlogis, ylim=c(-3,3))
> plsmo(age, L, group=sex, method='raw', add=TRUE, prefix='True', trim=0)
> title('Lowess-smoothed Estimates with True Regression Functions')
>
> dd <- datadist(age, sex, sys.bp)
> options(datadist='dd')
> # can also do: dd <- datadist(dd, newvar)
>
> f <- lrm(dz ~ rcs(age,5)*sex, x=TRUE, y=TRUE)
> f

Logistic Regression Model

lrm(formula = dz ~ rcs(age, 5) * sex, x = TRUE, y = TRUE)

Frequencies of Responses
0 1
235 265

Obs Max Deriv Model L.R. d.f. P C Dxy
500 1e-07 50.8 9 0 0.663 0.327
Gamma Tau-a R2 Brier
0.329 0.163 0.129 0.226

Coef S.E. Wald Z P
Intercept -7.9925 3.9390 -2.03 0.0425
age 0.1810 0.1029 1.76 0.0786
age' -0.3458 0.4323 -0.80 0.4238
age'' 1.5460 2.2145 0.70 0.4851
age''' -2.6579 3.3658 -0.79 0.4297
sex=male 14.0546 4.6913 3.00 0.0027
age * sex=male -0.3335 0.1228 -2.72 0.0066
age' * sex=male 0.8253 0.5303 1.56 0.1197
age'' * sex=male -2.8522 2.8360 -1.01 0.3146
age''' * sex=male 3.3676 4.5642 0.74 0.4606

> # x=TRUE, y=TRUE for pentrace
>
> fpred <- Function(f)
Warning in any(prm[1, -1]) :
coercing argument of type 'double' to logical
Warning in any(prm[2, -1]) :
coercing argument of type 'double' to logical
Warning in any(prm[1, -1]) :
coercing argument of type 'double' to logical
> fpred
function (age = 50.310462, sex = "male")
{
-7.9925102 + 0.18097811 * age - 0.00032804104 * pmax(age -
33.874461, 0)^3 + 0.0014667311 * pmax(age - 44.433245,
0)^3 - 0.0025216207 * pmax(age - 50.310462, 0)^3 + 0.0018029552 *
pmax(age - 55.835572, 0)^3 - 0.00042002458 * pmax(age -
66.340264, 0)^3 + 14.054577 * (sex == "male") + (sex ==
"male") * (-0.33352947 * age + 0.00078298646 * pmax(age -
33.874461, 0)^3 - 0.0027059801 * pmax(age - 44.433245,
0)^3 + 0.0031949484 * pmax(age - 50.310462, 0)^3 - 0.0016520917 *
pmax(age - 55.835572, 0)^3 + 0.00038013702 * pmax(age -
66.340264, 0)^3)
}
<environment: 0xa5a8ddc>
> fpred(age=30, sex=levels(sex))
[1] -2.56 1.49
>
> anova(f)
Wald Statistics Response: dz

Factor Chi-Square d.f. P
age (Factor+Higher Order Factors) 30.2 8 0.0002
All Interactions 16.6 4 0.0024
Nonlinear (Factor+Higher Order Factors) 23.3 6 0.0007
sex (Factor+Higher Order Factors) 25.6 5 0.0001
All Interactions 16.6 4 0.0024
age * sex (Factor+Higher Order Factors) 16.6 4 0.0024
Nonlinear 16.0 3 0.0011
Nonlinear Interaction : f(A,B) vs. AB 16.0 3 0.0011
TOTAL NONLINEAR 23.3 6 0.0007
TOTAL NONLINEAR + INTERACTION 23.5 7 0.0014
TOTAL 37.2 9 <.0001
>
> p <- plot(f, age=NA, sex=NA, conf.int=FALSE, ylim=c(-3,3))
> datadensity(p, age, sex)
> scat1d(age)
>
> plsmo(age, L, group=sex, method='raw', add=TRUE, prefix='True', trim=0)
> title('Spline Fits with True Regression Functions')
>
> f.bp <- lrm(dz.bp ~ rcs(age,5)*sex + rcs(sys.bp,5))
>
> for(method in c('persp','image'))
+ p <- plot(f.bp, age=NA, sys.bp=NA, method=method)
> # Legend(p) # NOTE: Needs subplot - not in R
>
> cat('Doing 25 bootstrap repetitions to validate model\n')
Doing 25 bootstrap repetitions to validate model
> validate(f, B=25) # in practice try to use 150
Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5
Iteration 6
Iteration 7
Iteration 8
Iteration 9
Iteration 10
Iteration 11
Iteration 12
Iteration 13
Iteration 14
Iteration 15
Iteration 16
Iteration 17
Iteration 18
Iteration 19
Iteration 20
Iteration 21
Iteration 22
Iteration 23
Iteration 24
Iteration 25

index.orig training test optimism index.corrected n
Dxy 0.3272 0.336 0.30141 0.03445 0.29272 25
R2 0.1290 0.142 0.11279 0.02953 0.09944 25
Intercept 0.0000 0.000 0.00886 -0.00886 0.00886 25
Slope 1.0000 1.000 0.86917 0.13083 0.86917 25
Emax 0.0000 0.000 0.03274 0.03274 0.03274 25
D 0.0996 0.111 0.08629 0.02470 0.07490 25
U -0.0040 -0.004 0.00186 -0.00586 0.00186 25
Q 0.1036 0.115 0.08444 0.03056 0.07305 25
B 0.2258 0.223 0.22953 -0.00637 0.23221 25
>
> cat('Doing 25 bootstrap reps to check model calibration\n')
Doing 25 bootstrap reps to check model calibration
> cal <- calibrate(f, B=25) # use 150 in practice
Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5
Iteration 6
Iteration 7
Iteration 8
Iteration 9
Iteration 10
Iteration 11
Iteration 12
Iteration 13
Iteration 14
Iteration 15
Iteration 16
Iteration 17
Iteration 18
Iteration 19
Iteration 20
Iteration 21
Iteration 22
Iteration 23
Iteration 24
Iteration 25

> plot(cal)

n=500 Mean absolute error=0.0185 Mean squared error=0.000536
0.9 Quantile of absolute error=0.0319

> title('Calibration of Unpenalized Model')
>
> p <- if(.R.) pentrace(f, penalty=c(.009,.009903,.02,.2,.5,1)) else
+ pentrace(f, penalty=1, method='optimize')
>
> f <- update(f, penalty=p$penalty)
> f

Logistic Regression Model

lrm(formula = dz ~ rcs(age, 5) * sex, x = TRUE, y = TRUE, penalty = p$penalty)

Frequencies of Responses
0 1
235 265

Penalty factors:

simple nonlinear interaction nonlinear.interaction
0.02 0.02 0.02 0.02

Final penalty on -2 log L: 1.71

Obs Max Deriv Model L.R. d.f. P C Dxy
500 2e-05 49.8 7.95 0 0.664 0.327
Gamma Tau-a R2 Brier
0.329 0.163 0.122 0.226

Coef S.E. Wald Z P Penalty Scale
Intercept -4.940899 2.47309 -2.00 0.0457 0.0000
age 0.099822 0.06297 1.59 0.1129 1.4410
age' -0.008624 0.24123 -0.04 0.9715 1.4759
age'' -0.049210 1.25333 -0.04 0.9687 0.4255
age''' -0.520005 2.06692 -0.25 0.8014 0.1320
sex=male 9.646487 2.96111 3.26 0.0011 0.1000
age * sex=male -0.216102 0.07559 -2.86 0.0042 3.6160
age' * sex=male 0.328440 0.30085 1.09 0.2750 1.0759
age'' * sex=male -0.454827 1.68334 -0.27 0.7870 0.2820
age''' * sex=male 0.077775 3.01079 0.03 0.9794 0.0841

> specs(f,long=TRUE)
Warning in any(i1) : coercing argument of type 'double' to logical
Warning in any(i2) : coercing argument of type 'double' to logical
Warning in any(i1) : coercing argument of type 'double' to logical
Warning in any(i2) : coercing argument of type 'double' to logical
Warning in any(i1) : coercing argument of type 'double' to logical
Warning in any(i2) : coercing argument of type 'double' to logical
lrm(formula = dz ~ rcs(age, 5) * sex, x = TRUE, y = TRUE, penalty = p$penalty)

Assumption Parameters d.f.
age rcspline 33.874 44.433 50.31 55.836 66.34 4
sex category female male 1
age * sex interaction nonlinear x linear - f(A)B 4

age sex
Low:effect 43.6 <NA>
Adjust to 50.3 male
High:effect 56.6 <NA>
Low:prediction 30.4 female
High:prediction 72.4 male
Low 15.1 female
High 89.4 male
> edf <- effective.df(f)

Original and Effective Degrees of Freedom

Original Penalized
All 9 7.95
Simple Terms 2 1.90
Interaction or Nonlinear 7 6.05
Nonlinear 6 5.12
Interaction 4 3.53
Nonlinear Interaction 3 2.61
>
> p <- plot(f, age=NA, sex=NA, conf.int=FALSE, ylim=c(-3,3))
> datadensity(p, age, sex)
> scat1d(age)
>
> plsmo(age, L, group=sex, method='raw', add=TRUE, prefix='True', trim=0)
> title('Penalized Spline Fits with True Regression Functions')
>
> options(digits=3)
> s <- summary(f)
> s
Effects Response : dz

Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
age 43.6 56.6 13.0 0.80 0.35 0.12 1.49
Odds Ratio 43.6 56.6 13.0 2.23 NA 1.12 4.42
sex - female:male 2.0 1.0 NA -0.07 0.31 -0.67 0.53
Odds Ratio 2.0 1.0 NA 0.93 NA 0.51 1.70

Adjusted to: age=50.3 sex=male

> plot(s)
>
> s <- summary(f, sex='male')
> plot(s)
>
> fpred <- Function(f)
Warning in any(prm[1, -1]) :
coercing argument of type 'double' to logical
Warning in any(prm[2, -1]) :
coercing argument of type 'double' to logical
Warning in any(prm[1, -1]) :
coercing argument of type 'double' to logical
> fpred
function (age = 50.310462, sex = "male")
{
-4.9408986 + 0.099821811 * age - 8.1818002e-06 * pmax(age -
33.874461, 0)^3 - 4.6687755e-05 * pmax(age - 44.433245,
0)^3 - 0.00049334968 * pmax(age - 50.310462, 0)^3 + 0.0008754865 *
pmax(age - 55.835572, 0)^3 - 0.00032726726 * pmax(age -
66.340264, 0)^3 + 9.6464867 * (sex == "male") + (sex ==
"male") * (-0.21610185 * age + 0.00031160482 * pmax(age -
33.874461, 0)^3 - 0.00043151337 * pmax(age - 44.433245,
0)^3 + 7.3788465e-05 * pmax(age - 50.310462, 0)^3 - 0.00017574468 *
pmax(age - 55.835572, 0)^3 + 0.00022186476 * pmax(age -
66.340264, 0)^3)
}
<environment: 0xb36d040>
> fpred(age=30, sex=levels(sex))
[1] -1.95 1.22
> sascode(fpred)
-4.9408986 + 0.099821811 * age - 8.1818002e-06 * max(age -
33.874461, 0)**3 - 4.6687755e-05 * max(age - 44.433245,
0)**3 - 0.00049334968 * max(age - 50.310462, 0)**3 + 0.0008754865 *
max(age - 55.835572, 0)**3 - 0.00032726726 * max(age -
66.340264, 0)**3 + 9.6464867 * (sex = "male") + (sex =
"male") * (-0.21610185 * age + 0.00031160482 * max(age -
33.874461, 0)**3 - 0.00043151337 * max(age - 44.433245,
0)**3 + 7.3788465e-05 * max(age - 50.310462, 0)**3 - 0.00017574468 *
max(age - 55.835572, 0)**3 + 0.00022186476 * max(age -
66.340264, 0)**3);
>
> cat('Doing 40 bootstrap reps to validate penalized model\n')
Doing 40 bootstrap reps to validate penalized model
> validate(f, B=40)
Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5
Iteration 6
Iteration 7
Iteration 8
Iteration 9
Iteration 10
Iteration 11
Iteration 12
Iteration 13
Iteration 14
Iteration 15
Iteration 16
Iteration 17
Iteration 18
Iteration 19
Iteration 20
Iteration 21
Iteration 22
Iteration 23
Iteration 24
Iteration 25
Iteration 26
Iteration 27
Iteration 28
Iteration 29
Iteration 30
Iteration 31
Iteration 32
Iteration 33
Iteration 34
Iteration 35
Iteration 36
Iteration 37
Iteration 38
Iteration 39
Iteration 40

index.orig training test optimism index.corrected n
Dxy 0.3273 0.339 0.29923 0.03963 0.28764 40
R2 0.1224 0.139 0.11053 0.02883 0.09356 40
Intercept 0.0000 0.000 0.00409 -0.00409 0.00409 40
Slope 1.0000 1.000 0.89841 0.10159 0.89841 40
Emax 0.0000 0.000 0.02459 0.02459 0.02459 40
D 0.0976 0.108 0.08447 0.02401 0.07356 40
U -0.0040 -0.004 0.00200 -0.00600 0.00200 40
Q 0.1016 0.112 0.08247 0.03001 0.07156 40
B 0.2263 0.222 0.23027 -0.00798 0.23425 40
>
> cat('Doing 40 bootstrap reps to check penalized model calibration\n')
Doing 40 bootstrap reps to check penalized model calibration
> cal <- calibrate(f, B=40)
Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5
Iteration 6
Iteration 7
Iteration 8
Iteration 9
Iteration 10
Iteration 11
Iteration 12
Iteration 13
Iteration 14
Iteration 15
Iteration 16
Iteration 17
Iteration 18
Iteration 19
Iteration 20
Iteration 21
Iteration 22
Iteration 23
Iteration 24
Iteration 25
Iteration 26
Iteration 27
Iteration 28
Iteration 29
Iteration 30
Iteration 31
Iteration 32
Iteration 33
Iteration 34
Iteration 35
Iteration 36
Iteration 37
Iteration 38
Iteration 39
Iteration 40

> plot(cal)

n=500 Mean absolute error=0.0170 Mean squared error=0.000442
0.9 Quantile of absolute error=0.0287

> title('Calibration of Penalized Model')
>
> nomogram(f.bp, fun=plogis,
+ funlabel='Prob(dz)',
+ fun.at=c(.15,.2,.3,.4,.5,.6,.7,.8,.9,.95,.975),
+ fun.side=c(1,3,1,3,1,3,1,3,1,3,1))
> options(datadist=NULL)
>
> #####################
> #Detailed Example 2 #
> #####################
> # Simulate the data.
> n <- 1000 # define sample size
> set.seed(17) # so can reproduce the results
> treat <- factor(sample(c('a','b','c'), n, TRUE))
> num.diseases <- sample(0:4, n, TRUE)
> age <- rnorm(n, 50, 10)
> cholesterol <- rnorm(n, 200, 25)
> weight <- rnorm(n, 150, 20)
> sex <- factor(sample(c('female','male'), n, TRUE))
> label(age) <- 'Age' # label is in Hmisc
> label(num.diseases) <- 'Number of Comorbid Diseases'
> label(cholesterol) <- 'Total Cholesterol'
> label(weight) <- 'Weight, lbs.'
> label(sex) <- 'Sex'
> units(cholesterol) <- 'mg/dl' # uses units.default in Hmisc
>
> # Specify population model for log odds that Y=1
> L <- .1*(num.diseases-2) + .045*(age-50) +
+ (log(cholesterol - 10)-5.2)*(-2*(treat=='a') +
+ 3.5*(treat=='b')+2*(treat=='c'))
> # Simulate binary y to have Prob(y=1) = 1/[1+exp(-L)]
> y <- ifelse(runif(n) < plogis(L), 1, 0)
> cholesterol[1:3] <- NA # 3 missings, at random
>
> ddist <- datadist(cholesterol, treat, num.diseases,
+ age, weight, sex)
> # Could have used ddist <- datadist(data.frame.name)
> options(datadist="ddist") # defines data dist. to Design
> cholesterol <- impute(cholesterol) # see impute in Hmisc library
> # impute, describe, and several other basic functions are
> # distributed as part of the Hmisc library
>
> fit <- lrm(y ~ treat*log(cholesterol - 10) +
+ scored(num.diseases) + rcs(age))
>
> describe(y ~ treat + scored(num.diseases) + rcs(age))
y ~ treat + scored(num.diseases) + rcs(age)

4 Variables 1000 Observations
--------------------------------------------------------------------------------
y
n missing unique Sum Mean
1000 0 2 522 0.522
--------------------------------------------------------------------------------
treat
n missing unique
1000 0 3

a (352, 35%), b (323, 32%), c (325, 32%)
--------------------------------------------------------------------------------
scored(num.diseases) : Number of Comorbid Diseases
n missing unique
1000 0 5

0 1 2 3 4
Frequency 187 197 221 231 164
% 19 20 22 23 16
--------------------------------------------------------------------------------
age : Age
n missing unique Mean .05 .10 .25 .50 .75 .90
1000 0 1000 49.93 33.04 37.05 43.06 49.56 57.05 63.12
.95
65.92

lowest : 21.43 22.64 23.60 23.89 25.13, highest: 74.64 74.90 77.43 77.74 80.32
--------------------------------------------------------------------------------
age' : Age
n missing unique Mean .05 .10 .25 .50
1000 0 951 8.578 4.214e-12 5.953e-02 9.286e-01 4.169e+00
.75 .90 .95
1.280e+01 2.415e+01 3.005e+01

lowest : 0.000e+00 4.436e-12 1.637e-06 8.878e-06 1.872e-05
highest: 4.854e+01 4.909e+01 5.444e+01 5.510e+01 6.057e+01
--------------------------------------------------------------------------------
age'' : Age
n missing unique Mean .05 .10 .25 .50 .75 .90
1000 0 726 1.702 0.0000 0.0000 0.0000 0.1735 2.1325 5.9507
.95
8.0651

lowest : 0.000e+00 1.347e-11 2.236e-09 1.647e-08 6.044e-07
highest: 1.471e+01 1.491e+01 1.684e+01 1.707e+01 1.904e+01
--------------------------------------------------------------------------------
age''' : Age
n missing unique Mean .05 .10 .25 .50
1000 0 501 0.4724 0.000e+00 0.000e+00 0.000e+00 3.646e-10
.75 .90 .95
3.877e-01 1.802e+00 2.644e+00

lowest : 0.000e+00 7.291e-10 9.830e-07 1.677e-06 2.905e-06
highest: 5.304e+00 5.384e+00 6.154e+00 6.248e+00 7.035e+00
--------------------------------------------------------------------------------
> # or use describe(formula(fit)) for all variables used in fit
> # describe function (in Hmisc) gets simple statistics on variables
> #fit <- robcov(fit) # Would make all statistics which follow
> # use a robust covariance matrix
> # would need x=TRUE, y=TRUE in lrm
> specs(fit) # Describe the design characteristics
Warning in any(i1) : coercing argument of type 'double' to logical
Warning in any(i2) : coercing argument of type 'double' to logical
Warning in any(i1) : coercing argument of type 'double' to logical
Warning in any(i2) : coercing argument of type 'double' to logical
Warning in any(i1) : coercing argument of type 'double' to logical
Warning in any(i2) : coercing argument of type 'double' to logical
lrm(formula = y ~ treat * log(cholesterol - 10) + scored(num.diseases) +
rcs(age))

Units Label Transformation
treat treat
cholesterol mg/dl Total Cholesterol log(cholesterol - 10)
num.diseases Number of Comorbid Diseases
age Age
treat * cholesterol treat * cholesterol
Assumption Parameters d.f.
treat category a b c 2
cholesterol asis 1
num.diseases scored 0 1 2 3 4 4
age rcspline 33.043 43.838 49.562 56.28 65.924 4
treat * cholesterol interaction linear x linear - AB 2
> a <- anova(fit)
> print(a, which='subscripts') # print which parameters being tested
Wald Statistics Response: y

Factor Chi-Square d.f. P
treat (Factor+Higher Order Factors) 15.86 4 0.0032
All Interactions 10.77 2 0.0046
cholesterol (Factor+Higher Order Factors) 12.76 3 0.0052
All Interactions 10.77 2 0.0046
num.diseases 14.92 4 0.0049
Nonlinear 0.73 3 0.8662
age 67.91 4 <.0001
Nonlinear 1.11 3 0.7736
treat * cholesterol (Factor+Higher Order Factors) 10.77 2 0.0046
TOTAL NONLINEAR 2.03 6 0.9167
TOTAL NONLINEAR + INTERACTION 13.19 8 0.1056
TOTAL 90.62 13 <.0001
Tested
1-2,12-13
12-13
3,12-13
12-13
4-7
5-7
8-11
9-11
12-13
5-7,9-11
5-7,9-13
1-13

Subscripts correspond to:
[1] treat=b treat=c cholesterol
[4] num.diseases num.diseases=2 num.diseases=3
[7] num.diseases=4 age age'
[10] age'' age''' treat=b * cholesterol
[13] treat=c * cholesterol
> plot(anova(fit)) # Depict Wald statistics graphically
> anova(fit, treat, cholesterol) # Test these 2 by themselves
Wald Statistics Response: y

Factor Chi-Square d.f. P
treat (Factor+Higher Order Factors) 15.9 4 0.0032
All Interactions 10.8 2 0.0046
cholesterol (Factor+Higher Order Factors) 12.8 3 0.0052
All Interactions 10.8 2 0.0046
TOTAL 17.7 5 0.0034
> summary(fit) # Estimate effects using default ranges
Effects Response : y

Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
cholesterol 183.8 217 32.9 -0.24 0.15 -0.54 0.06
Odds Ratio 183.8 217 32.9 0.79 NA 0.58 1.06
age 43.1 57 14.0 0.84 0.18 0.48 1.20
Odds Ratio 43.1 57 14.0 2.31 NA 1.61 3.30
treat - b:a 1.0 2 NA 0.40 0.16 0.08 0.73
Odds Ratio 1.0 2 NA 1.50 NA 1.09 2.07
treat - c:a 1.0 3 NA 0.20 0.16 -0.11 0.52
Odds Ratio 1.0 3 NA 1.23 NA 0.89 1.69
num.diseases - 0:2 3.0 1 NA -0.47 0.21 -0.89 -0.06
Odds Ratio 3.0 1 NA 0.62 NA 0.41 0.94
num.diseases - 1:2 3.0 2 NA -0.35 0.21 -0.76 0.05
Odds Ratio 3.0 2 NA 0.70 NA 0.47 1.06
num.diseases - 3:2 3.0 4 NA 0.08 0.20 -0.31 0.47
Odds Ratio 3.0 4 NA 1.08 NA 0.73 1.60
num.diseases - 4:2 3.0 5 NA 0.24 0.22 -0.20 0.67
Odds Ratio 3.0 5 NA 1.27 NA 0.82 1.96

Adjusted to: treat=a cholesterol=200

> plot(summary(fit)) # Graphical display of effects with C.L.
> summary(fit, treat="b", age=60)
Effects Response : y

Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
cholesterol 183.8 217 32.9 0.50 0.17 0.17 0.83
Odds Ratio 183.8 217 32.9 1.65 NA 1.18 2.30
age 43.1 57 14.0 0.84 0.18 0.48 1.20
Odds Ratio 43.1 57 14.0 2.31 NA 1.61 3.30
treat - a:b 2.0 1 NA -0.40 0.16 -0.73 -0.08
Odds Ratio 2.0 1 NA 0.67 NA 0.48 0.92
treat - c:b 2.0 3 NA -0.20 0.17 -0.53 0.13
Odds Ratio 2.0 3 NA 0.82 NA 0.59 1.14
num.diseases - 0:2 3.0 1 NA -0.47 0.21 -0.89 -0.06
Odds Ratio 3.0 1 NA 0.62 NA 0.41 0.94
num.diseases - 1:2 3.0 2 NA -0.35 0.21 -0.76 0.05
Odds Ratio 3.0 2 NA 0.70 NA 0.47 1.06
num.diseases - 3:2 3.0 4 NA 0.08 0.20 -0.31 0.47
Odds Ratio 3.0 4 NA 1.08 NA 0.73 1.60
num.diseases - 4:2 3.0 5 NA 0.24 0.22 -0.20 0.67
Odds Ratio 3.0 5 NA 1.27 NA 0.82 1.96

Adjusted to: treat=b cholesterol=200

> # Specify reference cell and adjustment val
>
> summary(fit, age=c(50,70)) # Estimate effect of increasing age from
Effects Response : y

Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
cholesterol 184 217 32.9 -0.24 0.15 -0.54 0.06
Odds Ratio 184 217 32.9 0.79 NA 0.58 1.06
age 50 70 20.0 1.18 0.32 0.54 1.81
Odds Ratio 50 70 20.0 3.25 NA 1.72 6.13
treat - b:a 1 2 NA 0.40 0.16 0.08 0.73
Odds Ratio 1 2 NA 1.50 NA 1.09 2.07
treat - c:a 1 3 NA 0.20 0.16 -0.11 0.52
Odds Ratio 1 3 NA 1.23 NA 0.89 1.69
num.diseases - 0:2 3 1 NA -0.47 0.21 -0.89 -0.06
Odds Ratio 3 1 NA 0.62 NA 0.41 0.94
num.diseases - 1:2 3 2 NA -0.35 0.21 -0.76 0.05
Odds Ratio 3 2 NA 0.70 NA 0.47 1.06
num.diseases - 3:2 3 4 NA 0.08 0.20 -0.31 0.47
Odds Ratio 3 4 NA 1.08 NA 0.73 1.60
num.diseases - 4:2 3 5 NA 0.24 0.22 -0.20 0.67
Odds Ratio 3 5 NA 1.27 NA 0.82 1.96

Adjusted to: treat=a cholesterol=200

> # 50 to 70
> summary(fit, age=c(50,60,70)) # Increase age from 50 to 70,
Effects Response : y

Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
cholesterol 184 217 32.9 -0.24 0.15 -0.54 0.06
Odds Ratio 184 217 32.9 0.79 NA 0.58 1.06
age 50 70 20.0 1.18 0.32 0.54 1.81
Odds Ratio 50 70 20.0 3.25 NA 1.72 6.13
treat - b:a 1 2 NA 0.40 0.16 0.08 0.73
Odds Ratio 1 2 NA 1.50 NA 1.09 2.07
treat - c:a 1 3 NA 0.20 0.16 -0.11 0.52
Odds Ratio 1 3 NA 1.23 NA 0.89 1.69
num.diseases - 0:2 3 1 NA -0.47 0.21 -0.89 -0.06
Odds Ratio 3 1 NA 0.62 NA 0.41 0.94
num.diseases - 1:2 3 2 NA -0.35 0.21 -0.76 0.05
Odds Ratio 3 2 NA 0.70 NA 0.47 1.06
num.diseases - 3:2 3 4 NA 0.08 0.20 -0.31 0.47
Odds Ratio 3 4 NA 1.08 NA 0.73 1.60
num.diseases - 4:2 3 5 NA 0.24 0.22 -0.20 0.67
Odds Ratio 3 5 NA 1.27 NA 0.82 1.96

Adjusted to: treat=a cholesterol=200

> # adjust to 60 when estimating
> # effects of other factors
> # If had not defined datadist, would have to define
> # ranges for all var.
>
> # Estimate and test treatment (b-a) effect averaged
> # over 3 cholesterols
> contrast(fit, list(treat='b',cholesterol=c(150,200,250)),
+ list(treat='a',cholesterol=c(150,200,250)),
+ type='average')
Contrast S.E. Lower Upper Z Pr(>|z|)
1 0.291 0.166 -0.0343 0.616 1.75 0.0796
> # Remove type='average' to get 3 separate contrasts for b-a
>
> # Plot effects. plot(fit) plots effects of all predictors,
> # showing values used for interacting factors as subtitles
> # The ref.zero parameter is helpful for showing effects of
> # predictors on a common scale for comparison of strength
> plot(fit, ref.zero=TRUE, ylim=c(-2,2))
>
> plot(fit, age=seq(20,80,length=100), treat=NA, conf.int=FALSE)
> # Plots relationship between age and log
> # odds, separate curve for each treat, no C.I.
> plot(fit, age=NA, cholesterol=NA)
> # 3-dimensional perspective plot for age, cholesterol, and
> # log odds using default ranges for both variables
> plot(fit, num.diseases=NA, fun=function(x) 1/(1+exp(-x)), #or fun=plogis
+ ylab="Prob", conf.int=.9)
> # Plot estimated probabilities instead of log odds
> # Again, if no datadist were defined, would have to
> # tell plot all limits
> logit <- predict(fit, expand.grid(treat="b",num.diseases=1:3,
+ age=c(20,40,60),
+ cholesterol=seq(100,300,length=10)))
> #logit <- predict(fit, gendata(fit, nobs=12))
> # Interactively specify 12 predictor combinations using UNIX
> # For UNIX or Windows, generate 9 combinations with other variables
> # set to defaults, get predicted values
> logit <- predict(fit, gendata(fit, age=c(20,40,60),
+ treat=c('a','b','c')))
>
> # Since age doesn't interact with anything, we can quickly and
> # interactively try various transformations of age,
> # taking the spline function of age as the gold standard. We are
> # seeking a linearizing transformation. Here age is linear in the
> # population so this is not very productive. Also, if we simplify the
> # model the total degrees of freedom will be too small and
> # confidence limits too narrow
>
> ag <- 10:80
> logit <- predict(fit, expand.grid(treat="a",
+ num.diseases=0, age=ag,
+ cholesterol=median(cholesterol)),
+ type="terms")[,"age"]
> # Note: if age interacted with anything, this would be the age
> # "main effect" ignoring interaction terms
> # Could also use
> # logit <- plot(f, age=ag, ...)$x.xbeta[,2]
> # which allows evaluation of the shape for any level
> # of interacting factors. When age does not interact with
> # anything, the result from
> # predict(f, ..., type="terms") would equal the result from
> # plot if all other terms were ignored
> # Could also use
> # logit <- predict(fit, gendata(fit, age=ag, cholesterol=median...))
>
> plot(ag^.5, logit) # try square root vs. spline transform.
> plot(ag^1.5, logit) # try 1.5 power
>
> # w <- latex(fit) # invokes latex.lrm, creates fit.tex
> # print(w) # display or print model on screen
>
> # Draw a nomogram for the model fit
> nomogram(fit, fun=plogis, funlabel="Prob[Y=1]")
>
> # Compose S function to evaluate linear predictors from fit
> g <- Function(fit)
Warning in any(prm[1, -1]) :
coercing argument of type 'double' to logical
Warning in any(prm[2, -1]) :
coercing argument of type 'double' to logical
Warning in any(prm[1, -1]) :
coercing argument of type 'double' to logical
> g(treat='b', cholesterol=260, age=50)
[1] 1.30
> # Leave num.diseases at reference value
>
> # Use the Hmisc dataRep function to summarize sample
> # sizes for subjects as cross-classified on 2 key
> # predictors
> drep <- dataRep(~ roundN(age,10) + num.diseases)
> print(drep, long=TRUE)

Data Representativeness n=1000

dataRep(formula = ~roundN(age, 10) + num.diseases)

Specifications for Matching

Type Parameters
age round to nearest 10
num.diseases exact numeric

Unique Combinations of Descriptor Variables

age num.diseases Frequency
2 30 0 11
3 40 0 44
4 50 0 68
5 60 0 53
6 70 0 11
8 20 1 1
9 30 1 10
10 40 1 57
11 50 1 75
12 60 1 43
13 70 1 11
15 20 2 1
16 30 2 21
17 40 2 47
18 50 2 85
19 60 2 43
20 70 2 22
21 80 2 2
22 20 3 2
23 30 3 17
24 40 3 48
25 50 3 88
26 60 3 63
27 70 3 13
30 30 4 8
31 40 4 47
32 50 4 61
33 60 4 40
34 70 4 7
35 80 4 1
>
> # Some approaches to making a plot showing how
> # predicted values vary with a continuous predictor
> # on the x-axis, with two other predictors varying
>
> fit <- lrm(y ~ log(cholesterol - 10) +
+ num.diseases + rcs(age) + rcs(weight) + sex)
>
> combos <- gendata(fit, age=10:100,
+ cholesterol=c(170,200,230),
+ weight=c(150,200,250))
> # num.diseases, sex not specified -> set to mode
> # can also used expand.grid
>
> combos$pred <- predict(fit, combos)
> library(lattice)
> xyplot(pred ~ age | cholesterol*weight, data=combos)
> xYplot(pred ~ age | cholesterol, groups=weight,
+ data=combos, type='l') # in Hmisc
Loading required package: grid
> xYplot(pred ~ age, groups=interaction(cholesterol,weight),
+ data=combos, type='l')
>
> # Can also do this with plot.Design but a single
> # plot may be busy:
> ch <- c(170, 200, 230)
> plot(fit, age=NA, cholesterol=ch, weight=150,
+ conf.int=FALSE)
> plot(fit, age=NA, cholesterol=ch, weight=200,
+ conf.int=FALSE, add=TRUE)
> plot(fit, age=NA, cholesterol=ch, weight=250,
+ conf.int=FALSE, add=TRUE)
>
> #Here we use plot.Design to make 9 separate plots, with CLs
> d <- expand.grid(cholesterol=c(170,200,230),
+ weight=c(150,200,250))
> for(i in 1:nrow(d)) {
+ plot(fit, age=NA, cholesterol=d$cholesterol[i],
+ weight=d$weight[i])
+ title(paste('Chol=',format(d$cholesterol[i]),' ',
+ 'Wt=',format(d$weight[i]),sep=''))
+ }
> options(datadist=NULL)
>
> ######################
> # Detailed Example 3 #
> ######################
> n <- 2000
> set.seed(731)
> age <- 50 + 12*rnorm(n)
> label(age) <- "Age"
> sex <- factor(sample(c('Male','Female'), n,
+ rep=TRUE, prob=c(.6, .4)))
> cens <- 15*runif(n)
> h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
> t <- -log(runif(n))/h
> label(t) <- 'Follow-up Time'
> e <- ifelse(t<=cens,1,0)
> t <- pmin(t, cens)
> units(t) <- "Year"
> age.dec <- cut2(age, g=10, levels.mean=TRUE)
> dd <- datadist(age, sex, age.dec)
> options(datadist='dd')
> Srv <- Surv(t,e)
>
> # Fit a model that doesn't assume anything except
> # that deciles are adequate representations of age
> f <- cph(Srv ~ strat(age.dec)+strat(sex), surv=TRUE)
> # surv=TRUE speeds up computations, and confidence limits when
> # there are no covariables are still accurate.
>
> # Plot log(-log 3-year survival probability) vs. mean age
> # within age deciles and vs. sex
> plot(f, age.dec=NA, sex=NA, time=3,
+ loglog=TRUE, val.lev=TRUE, ylim=c(-5,-1))
>
> # Fit a model assuming proportional hazards for age and
> # absence of age x sex interaction
> f <- cph(Srv ~ rcs(age,4)+strat(sex), surv=TRUE)
> survplot(f, sex=NA, n.risk=TRUE)
> # Add ,age=60 after sex=NA to tell survplot use age=60
> # Validate measures of model performance using the bootstrap
> # First must add data (design matrix and Srv) to fit object
> f <- update(f, x=TRUE, y=TRUE)
> validate(f, B=10, dxy=TRUE, u=5) # use t=5 for Dxy (only)
Iteration 1
Iteration 2
Iteration 3
Iteration 4
Iteration 5
Iteration 6
Iteration 7