|
|
|
||||
|
|
|||||
SPLUS Textbook Examples
Applied Survival Analysis by Hosmer and Lemeshow
Chapter 1: Introduction
Sections noted by R indicate differences between R and Splus.
Splus program.
R program.
Delimited data sets to be used in R:
hmohiv.csv To input the data sets in R use the following code which must be executed before running the R program:
>hmohiv <- read.csv("d:/hmohiv.csv", header=T)
Here is the same data set in SPLUS hmohiv.
Data analysis in SPLUS can be done either through a point and click interface or through syntax. In this page we have shown have to use both syntax and point and click.
Fig. 1.1, p. 6
Note1: The argument marks correspond to the argument pch
in the plot and points statements. Also,
the first two numbers are the coordinates of the upper left corner of the
legend.
Note2: SPLUS is case sensitive.
This means that the variable name (the part after the dollar sign) has to have
exactly the same case as in the dataset else the program will fail to run.
For simplicity we have changed all variable names in the data sets to be lower
case.
R-Note that in R the marks argument in the legend function is called pch and that the numbers for the symbols are different between the two packages.
plot(hmohiv$age[hmohiv$censor==1],
hmohiv$time[hmohiv$censor==1], xlab="Age",
ylab="Time" , pch=5, ylim=c(0, 65), xlim=c(15,55))
points(hmohiv$age[hmohiv$censor==0], hmohiv$time[hmohiv$censor==0], pch=18)
legend(40, 60, c("Censor=1", "Censor=0") , marks=c(5,18) )
Creating the new variable age1 which is equal to 1000/age.
hmohiv$age1 <- 1000/hmohiv$age
Point and Click:
Creating a new variable, age1, which is equal to 1000/age.
DATA
TRANSFORM
Make sure that you have hmohiv displayed in the data set field. Type age1 in the Target Column field, this will be the name of the new variable we are creating and in the Expression field type 1000/age which is the transformation we are applying to age. Now, click on the OK button and age1 will be the last variable in your data set.
Movie for creation of age1 using point and click.
Fig. 1.2, p. 6
plot(hmohiv$age1[hmohiv$censor==1],
hmohiv$time[hmohiv$censor==1], xlab="1000/Age",
ylab="Time" , pch=5, ylim=c(0, 65), xlim=c(15,55))
points(hmohiv$age1[hmohiv$censor==0], hmohiv$time[hmohiv$censor==0], pch=18)
legend(15, 60, c("Censor=1", "Censor=0") , marks=c(5,18) )
Table 1.2, p. 14
Note: In order to get the complete output we must first create the object "test" which is the
fitted survival regression and then get the summary of the object by using the summary command.
test <- survReg( Surv(time, censor) ~ age, hmohiv, dist="exponential")
summary(test)
Value Std.
Error z
p
(Intercept) 5.8590 0.5853
10.01 1.37e-023
age -0.0939
0.0158
-5.96 2.59e-009
Scale fixed at 1
Exponential distribution
Loglik(model)= -275 Loglik(intercept only)=
-292.3
Chisq= 34.5 on 1 degrees of freedom, p= 4.3e-009
Number of Newton-Raphson Iterations: 4
n= 100
Correlation of Coefficients:
(Intercept)
age -0.982
Point and Click:
STATISTICS
SURVIVAL
PARAMETRIC SURVIVAL...
In the Distribution field use the pull down menu and select "exponential". Then click on the Create Formula button which will bring up a new menu. Select time by clicking on time in the Choose Variables menu and then click on the Time1 button on the right. Then select censor and click on the Censor Codes button on the right. Then click on Add Response and you will see Surv(time, censor, type='right')~ 1 in the formula field. Now, we need to add the independent variable, age, to the formula. Select age in the Choose Variables menu and click on Main Effect(+) button in the Add Terms square. Now you will see Surv(time, censor, type='right')~ age in the formula field which is the formula that we would like use. So, click on the OK button at the bottom which brings you back to the original menu. Click the OK button and the survival regression will be executed.
Movie for table 1.2 and fig. 1.3 using point and click.
The list of all the distribution available in SPLUS.
names(survReg.distributions)
[1] "extreme" "logistic" "gaussian" "weibull"
"exponential"
[6] "rayleigh" "loggaussian" "lognormal" "loglogistic"
"t"
Obtaining the predicted values.
hmohiv$p <- predict(test, type="response")
Sorting time, age, censor and p by age. We want to order by ascending order of the variable age so we use age as the argument for the order function. Then we create the ordered variables age.1, time.1, censor.1 and p.1.
hmohiv$age.1<-hmohiv$age[order(hmohiv$age)]
hmohiv$time.1<-hmohiv$time[order(hmohiv$age)]
hmohiv$censor.1<-hmohiv$censor[order(hmohiv$age)]
hmohiv$p.1<-hmohiv$p[order(hmohiv$age)]
Point and Click:
Sorting time, age, censor and p by age.Note: The variables time.1, age.1, censor.1 and p.1 have already been created. If you choose to do the sort you will actually create four new variables namely time.2, age.2, censor.2 and p.2 which will be identical to time.1, age.1, censor.1 and p.1. SPLUS does not over-write these variables thereby allowing you to have variables sorting in multiple ways available at the same time.
DATA
RESTRUCTURE
SORT
First, make sure that you have hmohiv in the dataset slot. Select age, time, censor and p in the Columns window (use control key and click to select more than one variable) and select age in the Sort By Columns window. Click on the OK button at the bottom.
Movie of sorting time, age, censor and p by age using point and click.
Four new variables have been created: time.1, age.1, censor.1 and p.1 which have been sorted by age. These new variables will be used in the following graph.
Fig. 1.3, p. 16.
plot(hmohiv$age.1[hmohiv$censor.1==1], hmohiv$time.1[hmohiv$censor.1==1], xlab="Age",
ylab="Time" , pch=5, ylim=c(0, 65), xlim=c(15,55))
points(hmohiv$age.1[hmohiv$censor.1==0], hmohiv$time.1[hmohiv$censor.1==0], pch=18)
lines(hmohiv$age.1, hmohiv$p.1)
legend(40, 50, c("Censor=1", "Censor=0") , marks=c(5,18) )
UCLA Researchers are invited to our Statistical Consulting Services
We recommend others to our list of Other Resources for Statistical Computing Help
These pages are Copyrighted (c) by UCLA Academic Technology Services