*SAS and R code and minitab data for homework problems;

Copy and paste code one problem part at a time.


The easiest way to get mpack.txt and mrobdata.txt is to copy and paste the following 
two commands into R.

source("http://parker.ad.siu.edu/Olive/regpack.txt")
source("http://parker.ad.siu.edu/Olive/regdata.txt")


If mpack.txt and mrobdata.txt are on J drive, get lregpack and lregdata with
the following commands:

source("J:/regpack.txt") 
source("J:/regdata.txt")

#####for Linux use left most port for flash drive

source("/media/190A-4001/regpack.txt") 
source("/media/190A-4001/regdata.txt")


*problem 2.17;
options ls = 70;
data wcdata;
input x y;
cards;
30 73
20 50
60 128
80 170
40 87
50 108
60 135
30 60
70 148
60 132
;
proc print;
proc corr;
proc plot; plot y*x;
proc reg;
  model y=x;
  output out =a p = pred r = resid;
proc plot data = a;
plot resid*(pred x);
plot y*pred;
run;

*Problem 2.18;
options ls = 70;
data brand;
input  y x1 x2;
cards;
   64.0    4.0    2.0
   73.0    4.0    4.0
   61.0    4.0    2.0
   76.0    4.0    4.0
   72.0    6.0    2.0
   80.0    6.0    4.0
   71.0    6.0    2.0
   83.0    6.0    4.0
   83.0    8.0    2.0
   89.0    8.0    4.0
   86.0    8.0    2.0
   93.0    8.0    4.0
   88.0   10.0    2.0
   95.0   10.0    4.0
   94.0   10.0    2.0
  100.0   10.0    4.0
;
proc print;
proc corr;
proc plot; plot y*(x1 x2);
proc reg;
  model y=x1 x2;
  output out =a p = pred r = resid;
proc plot data = a;
plot resid*(pred x1 x2);
plot y*pred;
run;

#Problem 2.27: some R commands

a)

zbux <- cbind(buxx,buxy)
zbux <- as.data.frame(zbux)
zfull <- lm(buxy~len+nasal+bigonal+cephalic,data=zbux)
zred <- lm(buxy~len+nasal,data=zbux)
anova(zred,zfull)

d)

plot(zred$fit,buxy)
abline(0,1)

e)

plot(zred$fit,zred$resid)

f)

zbux <- zbux[-c(60,61,62,63,64,65),]
zfull <- lm(buxy~len+nasal+bigonal+cephalic,data=zbux)
zred <- lm(buxy~len+nasal,data=zbux)
anova(zred,zfull)

h)

plot(zred$fit,zbux[,5])
abline(0,1)

i)

plot(zred$fit,zred$resid)


#Problem 2.28 a) for R

#skewed data

x <- matrix(rnorm(3000),nrow=1000,ncol=3)
y <- x%*%c(1,1,1) + 2*rexp(1000)
out<-lsfit(x,y)
FIT <- y - out$resid
plot(FIT,y)
abline(0,1)
lines(lowess(FIT,y))

#Problem 2.28 b) and c) for R

plot(FIT,out$resid)
abline(0,0)
lines(lowess(FIT,out$resid))
out$coef

#Problem 3.5 commands for R

b) 

nx <- matrix(rnorm(300),nrow=100,ncol=3)
y <- exp( 4 + nx%*%c(1,1,1) + 0.5*rnorm(100) )

c) 

tplot(nx,y)

d)

out <- lsfit(nx,log(y))
ls.print(out)


#Problem 3.6 HW7 A commands for R

a)

zx <- cbrainx[,c(1,3,5,6,7,8,9,10)]
zbrain <- as.data.frame(cbind(cbrainy,zx))
zfull <- lm(cbrainy~.,data=zbrain)
summary(zfull)
back <- step(zfull)

c)

zint <- lm(cbrainy~1,data=zbrain)
forw <- step(zint,scope=list(lower=~1,
upper=~age+breadth+cephalic+circum+headht+height+len+sex),
direction="forward")

*Problem 3.30;
*From SAS User's Guide: Statistics (1985, p. 695-704,717-718);

data fitness;
      input Age Weight Oxygen RunTime RestPulse RunPulse MaxPulse @@;
      datalines;
   44 89.47 44.609 11.37 62 178 182   40 75.07 45.313 10.07 62 185 185
   44 85.84 54.297  8.65 45 156 168   42 68.15 59.571  8.17 40 166 172
   38 89.02 49.874  9.22 55 178 180   47 77.45 44.811 11.63 58 176 176
   40 75.98 45.681 11.95 70 176 180   43 81.19 49.091 10.85 64 162 170
   44 81.42 39.442 13.08 63 174 176   38 81.87 60.055  8.63 48 170 186
   44 73.03 50.541 10.13 45 168 168   45 87.66 37.388 14.03 56 186 192
   45 66.45 44.754 11.12 51 176 176   47 79.15 47.273 10.60 47 162 164
   54 83.12 51.855 10.33 50 166 170   49 81.42 49.156  8.95 44 180 185
   51 69.63 40.836 10.95 57 168 172   51 77.91 46.672 10.00 48 162 168
   48 91.63 46.774 10.25 48 162 164   49 73.37 50.388 10.08 67 168 168
   57 73.37 39.407 12.63 58 174 176   54 79.38 46.080 11.17 62 156 165
   52 76.32 45.441  9.63 48 164 166   50 70.87 54.625  8.92 48 146 155
   51 67.25 45.118 11.08 48 172 172   54 91.63 39.203 12.88 44 168 172
   51 73.71 45.790 10.47 59 186 188   57 59.08 50.545  9.93 49 148 155
   49 76.32 48.673  9.40 56 186 188   48 61.24 47.920 11.50 52 170 176
   52 82.78 47.467 10.50 53 170 172
   ;
   proc reg data=fitness;
    model Oxygen=Age Weight RunTime RunPulse RestPulse MaxPulse;
    output out =a p = pred r = resid;
      model Oxygen=Age Weight RunTime RunPulse RestPulse MaxPulse
            / selection=forward;
      model Oxygen=Age Weight RunTime RunPulse RestPulse MaxPulse
            / selection=backward;
      model Oxygen=Age Weight RunTime RunPulse RestPulse MaxPulse
            / selection=cp best = 10;
   run;

proc rsquare cp data = fitness;
model Oxygen=Age Weight RunTime RunPulse RestPulse MaxPulse;


proc plot data = a;
    plot resid*(pred);
    plot Oxygen*pred;

proc reg data=fitness;
    model Oxygen=Age  RunTime RunPulse MaxPulse;
    output out =sub p = pred r = resid;
 
proc plot data = sub;
    plot resid*(pred);
    plot Oxygen*pred;
run;

#Problem 4.1 commands for R

e <- rnorm(9)
x <- matrix(rnorm(27),nrow=9,ncol=3)
sqrtv <- sqrt(diag(1/1:9))
Y <- 4 + x%*%c(1,2,3) + sqrtv%*%e
wtt <- 1:9
lsfit(x,Y,wtt)$coef
kinv <- sqrt(diag(1:9))
Z <- kinv%*%Y
X <- 1 + 0*1:9
X <- cbind(X,x)
U <- kinv%*%X
lsfit(U,Z,int=F)$coef

#Problem 5.14 commands for R

#a)

help(warpbreaks)
out <- aov(breaks ~ tension, data = warpbreaks)
out
summary(out)
plot(out$fit,out$residuals)
title("Residual Plot")

#c)

warpbreaks[1,]
plot(out$fit,warpbreaks[,1])
abline(0,1)
title("Response Plot")


*5.16, from SAS User's Guide: Statistics (1985, p. 126-129);

*Uses proc glm to get residual and response plots;

options ls = 70;
data clover;
input strain $ nitrogen @@;
cards;
3dok1  19.4  3dok1  32.6  3dok1  27.0  3dok1  32.1  3dok1  33.0
3dok5  17.7  3dok5  24.8  3dok5  27.9  3dok5  25.2  3dok5  24.3
3dok4  17.0  3dok4  19.4  3dok4   9.1  3dok4  11.9  3dok4  15.8
3dok7  20.7  3dok7  21.0  3dok7  20.5  3dok7  18.8  3dok7  18.6
3dok13 14.3  3dok13 14.4  3dok13 11.8  3dok13 11.6  3dok13 14.2
compos 17.3  compos 19.4  compos 19.1  compos 16.9  compos 20.8
;
proc glm;
  class strain;
  model nitrogen=strain;
  output out =a p = pred r = resid;
proc plot data = a;
  plot resid*pred;
  plot nitrogen*pred;
run;


*Problem 6.3: The program below is for the two way Anova;
*Data from Montgomery (1984, p. 198); 

options ls = 70;
data voltage;
input material $ temp $ mvoltage @@;
cards;
1  50  130 1  50  155  1  50   74 1  50  180 
1  65   34 1  65   40  1  65   80 1  65   75  
1  80   20 1  80   70  1  80   82 1  80   58   
2  50  150 2  50  188  2  50  159 2  50  126 
2  65  136 2  65  122  2  65  106 2  65  115  
2  80   25 2  80   70  2  80   58 2  80   45   
3  50  138 3  50  110  3  50  168 3  50  160 
3  65  174 3  65  120  3  65  150 3  65  139  
3  80   96 3  80  104  3  80   82 3  80   60   
;
proc glm;
  class material temp;
  model mvoltage =material|temp;
  output out =a p = pred r = resid;
proc plot data = a;
  plot resid*pred;
  plot mvoltage*pred;
proc means data = voltage nway;
class material temp;
var mvoltage;
output out=means mean=ymn;
symbol1 v=square color=black i=join;
symbol2 v=circle color=black i=join;
symbol3 v=triangle color=black i=join;
proc gplot data=means;
title"Interaction Plot";
plot ymn * material = temp;
proc gplot data=means;
title"Interaction Plot";
plot ymn * temp = material;
run;


#R commands for problem 6.5

a)

out1<-aov(stime~ptype*treat,poison)
summary(out1)
out2<-aov(stime~ptype + treat + ptype*treat,poison)
summary(out2)
out3<-aov(stime~.^2,poison)
summary(out3)

b)

plot(fitted(out1),resid(out1))
title("Residual Plot")

c)

FIT <- poison$stime - out1$resid
plot(FIT,poison$stime)
abline(0,1)
title("Response Plot")

e)

attach(poison)
out4 <- aov((1/stime)~ptype*treat,poison)
summary(out4)

f)

plot(fitted(out4),resid(out4))
title("Residual Plot")

g)

FIT <- 1/poison$stime - out4$resid
plot(FIT,(1/poison$stime))
abline(0,1)
title("Response Plot")

h) 

interaction.plot(treat,ptype,(1/stime))
detach(poison)

*Problem 7.4 program  for a one way block design;
*Data from Box, Hunter and Hunter (2005, p. 146); 

options ls = 70;
data penicillan;
input block $ treat $ yield;
cards;
1  1  89
1  2  88
1  3  97
1  4  94
2  1  84
2  2  77
2  3  92
2  4  79
3  1  81
3  2  87
3  3  87
3  4  85
4  1  87
4  2  92
4  3  89
4  4  84
5  1  79
5  2  81
5  3  80
5  4  88
;
proc glm;
  class block treat;
  model yield =block treat;
  output out =a p = pred r = resid;
proc plot data = a;
  plot resid*pred;
  plot yield*pred;
run;


*Problem 7.5 program below is for a latin squares design;
*Data from Box, Hunter and Hunter (2005, p. 157-160); 

options ls = 70;
data auto;
input rblocks $ cblocks $ additives $ emmissions;
cards;
1  1  1  19
1  2  2  24
1  3  4  23
1  4  3  26
2  1  4  23
2  2  3  24
2  3  1  19
2  4  2  30
3  1  2  15
3  2  4  14
3  3  3  15
3  4  1  16
4  1  3  19
4  2  1  18
4  3  2  19
4  4  4  16
;
proc glm;
  class rblocks cblocks additives;
  model emmissions = rblocks cblocks additives;
  output out =a p = pred r = resid;
**the following commands can be used;
**to make a response and residual plot;
proc plot data = a;
  plot resid*pred;
  plot emmissions*pred;
run;

**Problem 8.9;
*Box, Hunter and Hunter (2005, p. 183) pilot plant data;
*SAS will treat a -1 like a 0;
*T = temp, C = concentration, K = catalyst;

options ls = 70;
data pilotplant;
input T $ C $ K $ yield;
cards;
-1  -1  -1  59
 1  -1  -1  74
-1   1  -1  50
 1   1  -1  69
-1  -1   1  50
 1  -1   1  81
-1   1   1  46
 1   1   1  79
-1  -1  -1  61
 1  -1  -1  70
-1   1  -1  58
 1   1  -1  67
-1  -1   1  54
 1  -1   1  85
-1   1   1  44
 1   1   1  81
;
proc glm;
  class T C K;
  model yield = T|C|K;
  output out =a p = pred r = resid;
proc plot data = a;
plot resid*pred;
plot yield*pred;
run;

*Problem 8.10;
*this data is actually for Minitab;
*A   B   C AB  AC  BC ABC yield;

-1  -1  -1  1   1   1  -1 59
 1  -1  -1 -1  -1   1   1 74
-1   1  -1 -1   1  -1   1 50
 1   1  -1  1  -1  -1  -1 69
-1  -1   1  1  -1  -1   1 50
 1  -1   1 -1   1  -1  -1 81
-1   1   1 -1  -1   1  -1 46
 1   1   1  1   1   1   1 79
-1  -1  -1  1   1   1  -1 61
 1  -1  -1 -1  -1   1   1 70
-1   1  -1 -1   1  -1   1 58
 1   1  -1  1  -1  -1  -1 67
-1  -1   1  1  -1  -1   1 54
 1  -1   1 -1   1  -1  -1 85
-1   1   1 -1  -1   1  -1 44
 1   1   1  1   1   1   1 81

#Problem 8.11 for R.
#Box, Hunter and Hunter (2005, p. 183 ) pilot plant data
mn <- 1 + 0*1:16
x1 <- rep(c(-1,-1,1,1),4)
x2 <- rep(c(-1,-1,-1,-1,1,1,1,1),2)
x3 <- mn
x3[1:8] <- -1
x12 <- x1*x2
x13 <- x1*x3
x23 <- x2*x3
x123 <- x12*x3
x<-cbind(x1,x2,x3,x12,x13,x23,x123)
y<-c(59,61,74,70,50,58,69,67,50,54,81,85,46,44,79,81)
out<-lsfit(x,y)
ls.print(out)
out2 <- aov(y~x1*x2*x3)
summary(out2)
summary.lm(out2)
rm(mn,x1,x2,x3,x12,x13,x23,x123,x,y,out,out2)


*Problem 9.7;
*Box, Hunter and Hunter (2005, p. 336) steel split plot data;
*Need to divide MS heat by MS wplots to get the correct
*F statistic to test heat;
options ls = 70;
data steel;
input coating $ heat $ repl $ resistance wplots $;
cards;
       1    1    1         67      1
       1    2    1         65      2
       1    3    1        155      3
       2    1    1         73      1
       2    2    1         91      2
       2    3    1        127      3
       3    1    1         83      1
       3    2    1         87      2
       3    3    1        147      3
       4    1    1         89      1
       4    2    1         86      2
       4    3    1        212      3
       1    1    2         33      4
       1    2    2        140      5
       1    3    2        108      6
       2    1    2          8      4
       2    2    2        142      5
       2    3    2        100      6
       3    1    2         46      4
       3    2    2        121      5
       3    3    2         90      6
       4    1    2         54      4
       4    2    2        150      5
       4    3    2        153      6
;
proc glm;
  class coating heat repl wplots;
  model resistance =  heat wplots coating heat*coating;
  test H = heat E = wplots;
run;


**Problem 9.8;
*SAS User's Guide: Statistics (1985, p. 131-132);
*Block 1 A 2 B 3 and 142 makes block = 1, A = 4, B = 2;
data split;
input Block 1 A 2 B 3 Response;
cards;
142 40.0
141 39.5
112 37.9
111 35.4 
121 36.7
122 38.2
132 36.4
131 34.8
221 42.7
222 41.6
212 40.3
211 41.6
241 44.5
242 47.6
231 43.6
232 42.8
;
**Use proc print to see what the data looks like;
*proc print;
proc anova;
class Block A B;
model Response = Block A Block*A B A*B;
test H=A E=Block*A;
run;

#problem 10.15 R code. Copy Gladstone data into R.
#a better model uses cau as a factor and add log(age) and log(size)
y <- cbrainx[,10]
age <- cbrainx[,1]
brd <- cbrainx[,3]
cau <- cbrainx[,4]
ceph <- cbrainx[,5]
circ <- cbrainx[,6]
hdht <- cbrainx[,7]
ht <- cbrainx[,8]
len <- cbrainx[,9]
siz <- cbrainx[,11]
outf <- glm(y~age+brd+cau+ceph+circ+hdht+ht+len+siz,family=binomial)
summary(outf)
drop1(outf,test="Chi")

#b) Test whether LR model is better than the null model
#could get G^2_o and G^2_full from anova(outf)
outn<-glm(y~1,family=binomial)
anova(outn,outf,test="Chi")

#c) reduced model with breadth, cause, height, size
outr<-glm(y~brd+cau+ht+siz,family=binomial)
anova(outr,outf,test="Chi")

#d)
ESP <- predict(outf)
ESPR <- predict(outr)
plot(ESPR,ESP)
abline(0,1)

#HW10.16  R code Copy Gladstone data into R.
y <- cbrainx[,10]
age <- cbrainx[,1]
lnag <- log(age)
brd <- cbrainx[,3]
fcau <- as.factor(cbrainx[,4])
ceph <- cbrainx[,5]
circ <- cbrainx[,6]
hdht <- cbrainx[,7]
ht <- cbrainx[,8]
len <- cbrainx[,9]
siz <- cbrainx[,11]
lnsz <- log(siz)
outf <- glm(y~age+lnag+brd+fcau+ceph+circ+hdht+ht+len+siz+lnsz,family=binomial)
summary(outf)
back <- step(outf)
outn <- glm(y~1,family=binomial)
AIC(outn)
AIC(outf)
#step(outf,direction="both")
#did backwards
#stepAIC and step would not do forwards
#forw <- step(outn,scope=list(lower~1,
#upper~age+lnag+brd+fcau+ceph+circ+hdht+ht+len+siz+lnsz),direction="forward")


*HW10.17;
*the 
*infile "J:cbrain2.txt"; 
*command needs to be changed if the file name is different
*or if the file is stored on a different drive;

options ls = 70;
data cbrain;
infile "J:cbrain2.txt";
input obs x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 sex x12;

proc rsquare;
  model sex = x2 x3 x4 x5 x6 x7 x9 x10 / cp;
  
proc reg;
  model sex = x2 x3 x4 x5 x6 x7 x9 x10;
  output out = a p = pred r = resid;
  run;


*hw 10.18;
options ls = 70;
data crabs;
* Agresti, p. 272;
input width cases satell;
cards;
22.69  14  5
23.84  14  4
24.77  28  17
25.84  39  21
26.79  22  15
27.74  24  20
28.67  18  15
30.41  14  14
;
proc logistic; model satell/cases = width;
     output out = predict p = pi_hat;
     proc print data = predict
run;


*HW 16.36, from Allison (1995, p. 49-50);

options ls = 70;
data heart;
input time status number;
cards;
25 1 16
25 0 3
75 1 11
75 0 0
150 1 4 
150 0 2
300 1 5
300 0 4
550 1 2
550 0 6
850 1 4
850 0 3
1150 1 1
1150 0 2
1450 1 1
1450 0 3
1700 1 0
1700 0 1
;
proc lifetest method = life intervals = 50 100 200 400 700
              1000 1300 1600 plots=(s,h) graphics outsurv=a;
     time time*status(0);
     freq number;
proc print data = a;
run;


*HW 16.37, from Allison (1995, p. 31-34);

options ls = 70;
data myel;
input time status;
cards;
8     1
180   1
632   1
852   0
52    1
2240  0
220   1
63    1
195   1
76    1
70    1
8     1
13    1
1990  0
1976  0
18    1
700   1
1296  0
1460  0
210   1
63    1
1328  0
1296  1
365   0
23    1
;
proc lifetest method = km  plots=(s,h) graphics outsurv=a;
     time time*status(0);
proc print data = a;
run;


*HW 16.38 data from Miller (1981, p. 49-50) and Smith (2002, p. 174); 

data Leukemia;
input Weeks Status Group @@;
      datalines;
   1 1 0  3 1 0  3 1 0  6 1 0  7 1 0 7 1 0
  10 1 0 12 1 0 14 1 0 15 1 0 18 1 0
  19 1 0 22 1 0 26 1 0 28 0 0 29 1 0
  34 1 0 40 1 0 48 0 0 49 0 0
   1 1 1  1 1 1  2 1 1  2 1 1  3 1 1
   4 1 1  5 1 1  8 1 1  8 1 1  9 1 1
  11 1 1 12 1 1 14 1 1 16 1 1 18 1 1
  21 1 1 27 0 1 31 1 1 38 0 1 44 1 1
 ;
 proc phreg data=Leukemia;
      model Weeks*Status(0)=Group;
   run;


*HW 16.39  data from SAS/STAT User's Guide (1999)
*www.helsinki.fi/atk/sas/sashtml/stat/chap49/index.htm;

data Myeloma;
      input Time VStatus LogBUN HGB Platelet Age LogWBC Frac
            LogPBM Protein SCalc;
      label Time='Survival Time'
            VStatus='0=Alive 1=Dead';
      datalines;
    1.25  1  2.2175   9.4  1  67  3.6628  1  1.9542  12  10
    1.25  1  1.9395  12.0  1  38  3.9868  1  1.9542  20  18
    2.00  1  1.5185   9.8  1  81  3.8751  1  2.0000   2  15
    2.00  1  1.7482  11.3  0  75  3.8062  1  1.2553   0  12
    2.00  1  1.3010   5.1  0  57  3.7243  1  2.0000   3   9
    3.00  1  1.5441   6.7  1  46  4.4757  0  1.9345  12  10
    5.00  1  2.2355  10.1  1  50  4.9542  1  1.6628   4   9
    5.00  1  1.6812   6.5  1  74  3.7324  0  1.7324   5   9
    6.00  1  1.3617   9.0  1  77  3.5441  0  1.4624   1   8
    6.00  1  2.1139  10.2  0  70  3.5441  1  1.3617   1   8
    6.00  1  1.1139   9.7  1  60  3.5185  1  1.3979   0  10
    6.00  1  1.4150  10.4  1  67  3.9294  1  1.6902   0   8
    7.00  1  1.9777   9.5  1  48  3.3617  1  1.5682   5  10
    7.00  1  1.0414   5.1  0  61  3.7324  1  2.0000   1  10
    7.00  1  1.1761  11.4  1  53  3.7243  1  1.5185   1  13
    9.00  1  1.7243   8.2  1  55  3.7993  1  1.7404   0  12
   11.00  1  1.1139  14.0  1  61  3.8808  1  1.2788   0  10
   11.00  1  1.2304  12.0  1  43  3.7709  1  1.1761   1   9
   11.00  1  1.3010  13.2  1  65  3.7993  1  1.8195   1  10
   11.00  1  1.5682   7.5  1  70  3.8865  0  1.6721   0  12
   11.00  1  1.0792   9.6  1  51  3.5051  1  1.9031   0   9
   13.00  1  0.7782   5.5  0  60  3.5798  1  1.3979   2  10
   14.00  1  1.3979  14.6  1  66  3.7243  1  1.2553   2  10
   15.00  1  1.6021  10.6  1  70  3.6902  1  1.4314   0  11
   16.00  1  1.3424   9.0  1  48  3.9345  1  2.0000   0  10
   16.00  1  1.3222   8.8  1  62  3.6990  1  0.6990  17  10
   17.00  1  1.2304  10.0  1  53  3.8808  1  1.4472   4   9
   17.00  1  1.5911  11.2  1  68  3.4314  0  1.6128   1  10
   18.00  1  1.4472   7.5  1  65  3.5682  0  0.9031   7   8
   19.00  1  1.0792  14.4  1  51  3.9191  1  2.0000   6  15
   19.00  1  1.2553   7.5  0  60  3.7924  1  1.9294   5   9
   24.00  1  1.3010  14.6  1  56  4.0899  1  0.4771   0   9
   25.00  1  1.0000  12.4  1  67  3.8195  1  1.6435   0  10
   26.00  1  1.2304  11.2  1  49  3.6021  1  2.0000  27  11
   32.00  1  1.3222  10.6  1  46  3.6990  1  1.6335   1   9
   35.00  1  1.1139   7.0  0  48  3.6532  1  1.1761   4  10
   37.00  1  1.6021  11.0  1  63  3.9542  0  1.2041   7   9
   41.00  1  1.0000  10.2  1  69  3.4771  1  1.4771   6  10
   41.00  1  1.1461   5.0  1  70  3.5185  1  1.3424   0   9
   51.00  1  1.5682   7.7  0  74  3.4150  1  1.0414   4  13
   52.00  1  1.0000  10.1  1  60  3.8573  1  1.6532   4  10
   54.00  1  1.2553   9.0  1  49  3.7243  1  1.6990   2  10
   58.00  1  1.2041  12.1  1  42  3.6990  1  1.5798  22  10
   66.00  1  1.4472   6.6  1  59  3.7853  1  1.8195   0   9
   67.00  1  1.3222  12.8  1  52  3.6435  1  1.0414   1  10
   88.00  1  1.1761  10.6  1  47  3.5563  0  1.7559  21   9
   89.00  1  1.3222  14.0  1  63  3.6532  1  1.6232   1   9
   92.00  1  1.4314  11.0  1  58  4.0755  1  1.4150   4  11
    4.00  0  1.9542  10.2  1  59  4.0453  0  0.7782  12  10
    4.00  0  1.9243  10.0  1  49  3.9590  0  1.6232   0  13
    7.00  0  1.1139  12.4  1  48  3.7993  1  1.8573   0  10
    7.00  0  1.5315  10.2  1  81  3.5911  0  1.8808   0  11
    8.00  0  1.0792   9.9  1  57  3.8325  1  1.6532   0   8
   12.00  0  1.1461  11.6  1  46  3.6435  0  1.1461   0   7
   11.00  0  1.6128  14.0  1  60  3.7324  1  1.8451   3   9
   12.00  0  1.3979   8.8  1  66  3.8388  1  1.3617   0   9
   13.00  0  1.6628   4.9  0  71  3.6435  0  1.7924   0   9
   16.00  0  1.1461  13.0  1  55  3.8573  0  0.9031   0   9
   19.00  0  1.3222  13.0  1  59  3.7709  1  2.0000   1  10
   19.00  0  1.3222  10.8  1  69  3.8808  1  1.5185   0  10
   28.00  0  1.2304   7.3  1  82  3.7482  1  1.6721   0   9
   41.00  0  1.7559  12.8  1  72  3.7243  1  1.4472   1   9
   53.00  0  1.1139  12.0  1  66  3.6128  1  2.0000   1  11
   57.00  0  1.2553  12.5  1  66  3.9685  0  1.9542   0  11
   77.00  0  1.0792  14.0  1  60  3.6812  0  0.9542   0  12
   ;
proc phreg data=Myeloma;
      model Time*VStatus(0)=LogBUN HGB Platelet Age LogWBC
                            Frac LogPBM Protein SCalc
                            / selection=backward 
                              slstay=0.15 details;
   run;
proc phreg data=Myeloma;
      model Time*VStatus(0)=LogBUN HGB Platelet Age LogWBC
                            Frac LogPBM Protein SCalc
                            / selection=forward slentry=0.25
                              details;
   run;
   proc phreg data=Myeloma;
      model Time*VStatus(0)=LogBUN HGB Platelet Age LogWBC
                            Frac LogPBM Protein SCalc
                            / selection=stepwise slentry=0.25
                              slstay=0.15 details;
   run;
proc phreg data=Myeloma;
      model Time*VStatus(0)=LogBUN HGB Platelet Age LogWBC
                            Frac LogPBM Protein SCalc
                            / selection=score best=3;
   run;



*HW 16.40 Data described in Chapter 3 of Allison (1999);
*need to save recid.txt on e: drive;

data recid;
infile "e:recid.txt";
input week arrest fin age race wexp mar paro prio educ emp1-emp52;
proc phreg data = recid;
model week*arrest(0)=fin age race wexp mar paro prio;
run;

proc phreg data = recid;
model week*arrest(0)=fin age race wexp mar paro prio
 / selection=backward 
                              slstay=0.15 details;
run;

proc phreg data = recid;
model week*arrest(0)=fin age race wexp mar paro prio
                            / selection=forward slentry=0.25
                              details;
run;

proc phreg data = recid;
model week*arrest(0)=fin age race wexp mar paro prio  
                            / selection=stepwise slentry=0.25
                              slstay=0.15 details;
run;

proc phreg data = recid;
model week*arrest(0)=fin age race wexp mar paro prio
                            / selection=score best=3;
run;


* HW 16.41  Collett (2003, p. 344-346);
* use model survtime*status(0)=;
* to get log L(none)

data ovarian;
input patient survtime status treat age;
datalines;
1  156  1  1  66
2  1040 0  1  38
3   59  1  1  72
4  421  0  2  53
5  329  1  1  43
6  769  0  2  59
7  365  1  2  64
8  770  0  2  57
9  1227 0  2  59
10 268  1  1  74
11 475  1  2  59
12 1129 0  2  53
13 464  1  2  56
14 1206 0  2  44
15 638  1  1  56
16 563  1  2  55
17 1106 0  1  44
18 431  1  1  50
19 855  0  1  43
20 803  0  1  39
21 115  1  1  74
22 744  0  2  50
23 477  0  1  64
24 448  0  1  56
25 353  1  2  63
26 377  0  2  58
;
proc lifereg data = ovarian;
class treat;
model survtime*status(0)= age treat;
run;


# HW 16.42 Allison(1995, p. 158);

data myel;
input id dur status treat renal;
cards;
1  8 1 1 1
2 180 1 2 0
3 632 1 2 0
4 852 0 1 0
5 52  1 1 1
6 2240 0 2 0
7 220 1 1  0
8 63  1 1  1
9 195 1 2 0
10 76 1 2 0
11 70 1 2 0
12 8 1 1 0
13 13 1 2 1
14 1990 0 2 0
15 1976 0 1 0
16 18 1 2 1
17 700 1 2 0
18 1296 0 1 0
19 1460 0 1 0
20  210 1 2 0
21  63 1 1 1
22  1328 0 1 0
23  1296 1 2 0
24  365 0 1 0
25  23 1 2 1
;
proc phreg data=myel;
   model dur*status(0)=treat / ties=exact;
   strata renal;
run;
proc phreg data=myel;
   model dur*status(0)=treat / ties=exact;
run;
proc phreg data=myel;
   model dur*status(0)=treat renal / ties=exact;
run;


*16.60 (= HW11.1)   See Allison (1995, p. 156);

data recid;
infile "e:recid.txt";
input week arrest fin age race wexp mar paro prio educ emp1-emp52;
proc phreg data = recid;
model week*arrest(0)=fin age race wexp mar paro prio finlt agelt racelt
wexplt marlt parolt priolt / ties =efron;
finlt = fin*log(week);
agelt = age*log(week);
racelt = race*log(week);
wexplt = wexp*log(week);
marlt = mar*log(week);
parolt = paro*log(week);
priolt = prio*log(week);
proc phreg data = recid;
model week*arrest(0)=fin age race wexp mar paro prio / ties =efron;
run;



*16.61 (=HW11.2) Allison (1995, ch 5);

options yearcutoff=1900;
data stan;
input dob mmddyy9. doa mmddyy9. dot mmddyy9. dls mmddyy9.
   id age dead dur surg trans wtime m1 m2 m3 reject;
surv1=dls-doa;
surv2=dls-dot;
wait=dot-doa;
agetrans=(dot-dob)/365.25;
ageaccpt=(doa-dob)/365.25;
if wait = . then wait = 10000;
agels=(dls-dob)/365.25;
cards;
01/10/37 11/15/67     .    01/03/68 1 30 1 50 0 0 . . . . .
03/02/16 01/02/68     .    01/07/68 2 51 1 6 0 0 . . . . .
09/19/13 01/06/68 01/06/68 01/21/68 3 54 1 16 0 1 1 2 0 1.110 0
12/23/27 03/28/68 05/02/68 05/05/68 4 40 1 39 0 1 36 3 0 1.660 0
07/28/47 05/10/68     .    05/27/68 5 20 1 18 0 0 . . . . .
11/08/13 06/13/68     .    06/15/68 6 54 1 3 0 0 . . . . .
08/29/17 07/12/68 08/31/68 05/17/70 7 50 1 675 0 1 51 4 0 1.320 1
03/27/23 08/01/68     .    09/09/68 8 45 1 40 0 0 . . . . .
06/11/21 08/09/68     .    11/01/68 9 47 1 85 0 0 . . . . .
02/09/26 08/11/68 08/22/68 10/07/68 10 42 1 58 0 1 12 2 0 0.610 1
08/22/20 08/15/68 09/09/68 01/14/69 11 47 1 153 0 1 26 1 0 0.360 0
07/09/15 09/17/68     .    09/24/68 12 53 1 8 0 0 . . . . .
02/22/14 09/19/68 10/05/68 12/08/68 13 54 1 81 0 1 17 3 0 1.890 1
09/16/14 09/20/68 10/26/68 07/07/72 14 53 1 1386 0 1 37 1 0 0.870 1
12/04/14 09/27/68     .    09/27/68 15 53 1 1 1 0 . . . . .
05/16/19 10/26/68 11/22/68 08/29/69 16 49 1 308 0 1 28 2 0 1.120 1
06/29/48 10/28/68     .    12/02/68 17 20 1 36 0 0 . . . . .
12/27/11 11/01/68 11/20/68 12/13/68 18 56 1 43 0 1 20 3 0 2.050 0
10/04/09 11/18/68     .    12/24/68 19 59 1 37 0 0 . . . . .
10/19/13 01/29/69 02/15/69 02/25/69 20 55 1 28 0 1 18 3 1 2.760 1
09/29/25 02/01/69 02/08/69 11/29/71 21 43 1 1032 0 1 8 2 0 1.130 1
06/05/26 03/18/69 03/29/69 05/07/69 22 42 1 51 0 1 12 3 0 1.380 1
12/02/10 04/11/69 04/13/69 04/13/71 23 58 1 733 0 1 3 3 0 0.960 1
07/07/17 04/25/69 07/16/69 11/29/69 24 51 1 219 0 1 83 3 1 1.620 1
02/06/36 04/28/69 05/22/69 04/01/74 25 33 0 1799 0 1 25 2 0 1.060 0
10/18/38 05/01/69     .    03/01/73 26 30 0 1400 0 0 . . . . .
07/21/60 05/04/69     .    01/21/70 27 8 1 263 0 0 . . . . .
05/30/15 06/07/69 08/16/69 08/17/69 28 53 1 72 0 1 71 2 0 0.470 0
02/06/19 07/14/69     .    08/17/69 29 50 1 35 0 0 . . . . .
09/20/24 08/19/69 09/03/69 12/18/71 30 44 1 852 0 1 16 4 0 1.580 1
10/04/14 08/23/69     .    09/07/69 31 54 1 16 0 0 . . . . .
04/02/05 08/29/69 09/14/69 11/13/69 32 64 1 77 0 1 17 4 0 0.690 1
01/01/21 11/27/69 01/16/70 04/01/74 33 48 0 1586 0 1 51 3 0 0.910 0
05/24/29 12/12/69 01/03/70 04/01/74 34 40 0 1571 0 1 23 2 0 0.380 0
08/04/26 01/21/70     .    02/01/70 35 43 1 12 0 0 . . . . .
05/01/21 04/04/70 05/19/70 07/12/70 36 48 1 100 0 1 46 2 0 2.090 1
10/24/08 04/25/70 05/13/70 06/29/70 37 61 1 66 0 1 19 3 1 0.870 1
11/14/28 05/05/70 05/08/70 05/09/70 38 41 1 5 0 1 5 3 0 0.870 0
11/12/19 05/20/70 05/21/70 07/11/70 39 50 1 53 0 1 2 . . . .
11/30/21 05/25/70 07/04/70 04/01/74 40 48 0 1407 1 1 41 4 0 0.750 0
04/30/25 08/19/70 10/15/70 04/01/74 41 45 0 1321 1 1 58 2 0 0.980 0
03/13/34 08/21/70     .    08/23/70 42 36 1 3 0 0 . . . . .
06/01/27 10/22/70     .    10/23/70 43 43 1 2 1 0 . . . . .
05/02/28 11/30/70     .    01/08/71 44 42 1 40 1 0 . . . . .
10/30/34 01/05/71 01/05/71 02/18/71 45 36 1 45 0 1 1 1 0 0.000 0
06/01/22 01/10/71 01/11/71 10/01/73 46 48 1 995 1 1 2 2 0 0.810 1
12/28/23 02/02/71 02/22/71 04/14/71 47 47 1 72 0 1 21 3 0 1.380 1
01/23/15 02/05/71     .    02/13/71 48 56 1 9 0 0 . . . . .
06/21/34 02/15/71 03/22/71 04/01/74 49 36 0 1141 1 1 36 4 0 1.350 0
03/28/25 02/15/71 05/08/71 10/21/73 50 45 1 979 1 1 83 . . . .
06/29/22 03/24/71 04/24/71 01/02/72 51 48 1 285 0 1 32 4 1 1.080 1
01/24/30 04/25/71     .    08/04/71 52 41 1 102 0 0 . . . . .
02/27/24 07/02/71 08/11/71 01/05/72 53 47 1 188 0 1 41 . . . .
09/16/23 07/02/71     .    07/04/71 54 47 1 3 0 0 . . . . .
02/24/19 08/09/71 08/18/71 10/08/71 55 52 1 61 0 1 10 2 0 1.510 1
12/05/32 09/03/71 11/08/71 04/01/74 56 38 0 941 0 1 67 4 0 0.980 0
06/08/30 09/13/71     .    02/08/72 57 41 1 149 0 0 . . . . .
09/17/23 09/23/71 10/13/71 08/30/72 58 47 1 342 1 1 21 2 1 1.820 1
05/12/30 09/29/71 12/15/71 04/01/74 59 41 0 915 1 1 78 2 0 0.190 0
10/29/22 11/18/71 11/20/71 01/24/72 60 49 1 68 0 1 3 3 0 0.660 1
05/12/19 12/04/71     .    12/05/71 61 52 1 2 0 0 . . . . .
08/01/32 12/09/71     .    02/15/72 62 39 1 69 0 0 . . . . .
04/15/39 12/12/71 01/07/72 04/01/74 63 32 0 841 0 1 27 3 1 1.930 0
04/09/23 02/01/72 03/04/72 09/06/73 64 48 1 583 1 1 32 1 0 0.120 0
11/19/20 03/06/72 03/17/72 05/22/72 65 51 1 78 0 1 12 2 0 1.120 1
01/02/19 03/20/72     .    04/20/72 66 53 1 32 0 0 . . . . .
09/03/52 03/23/72 05/18/72 01/01/73 67 19 1 285 0 1 57 3 0 1.020 0
01/10/27 04/07/72 04/09/72 06/13/72 68 45 1 68 0 1 3 3 1 1.680 1
06/05/24 06/01/72 06/10/72 04/01/74 69 47 0 670 0 1 10 2 0 1.200 0
06/17/19 06/17/72 06/21/72 07/16/72 70 52 1 30 0 1 5 3 1 1.680 1
02/22/25 07/21/72 08/20/72 04/01/74 71 47 0 620 0 1 31 3 0 0.970 0
11/22/45 08/14/72 08/17/72 04/01/74 72 26 0 596 0 1 4 3 1 1.460 0
05/13/16 09/11/72 10/07/72 12/09/72 73 56 1 90 0 1 27 3 1 2.160 1
07/20/43 09/18/72 09/22/72 10/04/72 74 29 1 17 0 1 5 1 0 0.610 0
07/25/20 09/29/72     .    09/30/72 75 52 1 2 0 0 . . . . .
09/03/20 10/04/72 11/18/72 04/01/74 76 52 0 545 1 1 46 3 1 1.700 0
08/27/31 10/06/72     .    10/26/72 77 41 1 21 0 0 . . . . .
02/20/24 11/03/72 05/31/73 04/01/74 78 48 0 515 0 1 210 3 0 0.810 0
02/18/19 11/30/72 02/04/73 03/05/73 79 53 1 96 0 1 67 2 0 1.080 1
06/27/26 12/06/72 12/31/72 04/01/74 80 46 0 482 1 1 26 3 0 1.410 0
02/21/20 01/12/73 01/17/73 04/01/74 81 52 0 445 0 1 6 4 1 1.940 0
08/19/42 11/01/71     .    01/01/73 82 29 0 427 0 0 . . . . .
10/04/19 01/24/73 02/24/73 04/13/73 83 53 1 80 0 1 32 4 0 3.050 0
05/13/30 01/30/73 03/07/73 12/29/73 84 42 1 334 0 1 37 4 0 0.600 1
02/13/25 02/06/73     .    02/10/73 85 47 1 5 0 0 . . . . .
03/30/24 03/01/73 03/08/73 04/01/74 86 48 0 397 0 1 8 3 1 1.440 0
12/19/26 03/21/73 05/19/73 07/08/73 87 46 1 110 0 1 60 2 0 2.250 1
11/16/18 03/28/73 04/27/73 04/01/74 88 54 0 370 0 1 31 3 0 0.680 0
03/19/22 04/05/73 08/21/73 10/28/73 89 51 1 207 0 1 139 4 1 1.330 1
03/25/21 04/06/73 09/12/73 10/08/73 90 52 1 186 1 1 160 3 1 0.820 0
09/08/25 04/13/73     .    03/18/74 91 47 1 340 0 0 . . . . .
05/03/28 04/27/73 03/02/74 04/01/74 92 44 0 340 0 1 310 1 0 0.160 0
10/10/25 07/11/73 08/07/73 04/01/74 93 47 0 265 0 1 28 2 0 0.330 0
11/11/29 09/14/73 09/17/73 02/25/74 94 43 1 165 1 1 4 3 0 1.200 1
06/11/33 09/22/73 09/23/73 10/07/73 95 40 1 16 0 1 2 . . . .
02/09/47 10/04/73 10/16/73 04/01/74 96 26 0 180 0 1 13 2 0 0.460 0
04/11/50 11/22/73 12/12/73 04/01/74 97 23 0 131 0 1 21 3 1 1.780 0
04/28/45 12/14/73 03/19/74 04/01/74 98 28 0 109 0 1 96 4 1 0.770 0
02/24/24 12/25/73     .    01/14/74 99 49 1 21 0 0 . . . . .
01/31/39 02/22/74 03/31/74 04/01/74 100 35 0 39 1 1 38 3 0 0.670 0
08/25/24 03/02/74     .    04/01/74 101 49 0 31 0 0 . . . . .
10/30/33 03/22/74     .    04/01/74 102 40 0 11 0 0 . . . . .
05/20/28 09/13/67     .    09/18/67 103 39 1 6 0 0 . . . . .
;

data stan;
   input dob mmddyy9. doa mmddyy9. dot mmddyy9. dls mmddyy9.
         dead surg m1 m2 m3;
   surv1=dls-doa;
   surv2=dls-dot;
   ageaccpt=(doa-dob)/365.25;
   agetrans=(dot-dob)/365.25;
   wait=dot-doa;
   if dot=. then trans=0; else trans=1;
run;

proc phreg data=stan;
   model surv1*dead(0)=plant surg ageaccpt / ties=exact;
   if wait>surv1 or wait=. then plant=0; else plant=1;
run;




####R code for ch 16 problems

#HW 16.43 Miller p. 49 time until relapse under maintenance chemotherapy
#AML leukemia data
library(survival)
time <- c(9,13,13,18,23,28,31,34,45,48,161)
status <- c(1,1,0,1,1,0,1,1,0,1,0)
z <- survfit(Surv(time,status)~1,conf.type="plain")
summary(z)
#plot KM estimator hat(S)_K(t) and the CIs for S(t)
plot(z)


#HW 16.44a type library(survival) if necessary
zy <- rexp(20)
zz <- rexp(20,rate=.2)
zt <- pmin(zy,zz)
zdelta <- as.numeric(zz>zy)
out <- survfit(Surv(zt, zdelta)~1, conf.type = "plain")
st <- exp( - out$time)
plot(out)
points(out$time,st)

#HW 16.44b type library(survival) if necessary
#repeat with n = 200
zy <- rexp(200)
zz <- rexp(200,rate=.2)
zt <- pmin(zy,zz)
zdelta <- as.numeric(zz>zy)
out <- survfit(Surv(zt, zdelta)~1, conf.type = "plain")
st <- exp( - out$time)
plot(out)
points(out$time,st)

#HW16.45
source("http://lagrange.math.siu.edu/Olive/regpack.txt")
library(survival)
kmsim2(10)

#HW 16.46a

library(survival)
z <- coxph(Surv(futime,fustat)~age + ecog.ps, ovarian)
z
plot(survfit(z))

# 16.46b
new<-data.frame(futime=365,fustat=1,age=66,resid.ds=2,rx=1,ecog.ps=4)
plot(survfit(z,newdata=new))


#HW 16.48

library(survival)

zfull <- coxph(Surv(time,status)~strata(sex) + age +ph.ecog + ph.karno+pat.karno+wt.loss, 
data = lung,na.action=na.omit)

zred1 <- coxph(Surv(time,status)~strata(sex) +  ph.ecog + ph.karno+pat.karno+wt.loss, 
data = lung,na.action=na.omit)

zred2 <- coxph(formula = Surv(time, status) ~ strata(sex) + ph.ecog + 
    +pat.karno + wt.loss, data = lung, na.action = na.omit)




#HW 16.49 need to download bphgfit function, Smith (2002, p. 174) Lawless (1982)
#leukemia data

lweeks <- c(1, 3, 3, 6, 7, 7, 10, 12, 14, 15, 18, 19, 22, 26, 28, 29, 
  34,  40, 48,  49,  1, 1, 2, 2, 3, 4, 5, 8, 8, 9, 11, 12, 14, 16, 18,
  21,  27, 31,  38,  44 )
 
lstatus <- 1 + 0*1:40

lstatus[c(15,19,20,37,39)] <- 0

lg <- 0*1:20
lg2 <- lg + 1
lgroup <- c(lg,lg2)

leuk <-  cbind(lweeks,lstatus,lgroup)

rm(lweeks,lstatus,lg,lg2,lgroup)

zout <- coxph(Surv(leuk[,1], leuk[,2]) ~ leuk[,3])

zout

bphgfit(leuk[,1],leuk[,2],leuk[,3])


# 16.50 a) status = 1 + usual status variable

#library(survival) if necessary

flung<-lung[lung$sex==1,]
nflung<-na.omit(flung[,c("time","status","ph.ecog","ph.karno",
"pat.karno","wt.loss")])
out <- coxph(Surv(time, status) 
~ ph.ecog + ph.karno + pat.karno + wt.loss,data = nflung)
        x <- as.matrix(nflung[,c(3,4,5,6)])
        ESP <- x %*% out$coef
        Time <- as.vector(nflung[,1])
        status <- as.vector(nflung[,2])
        plot(ESP,Time,pch = ' ')
        points(ESP[status == 2], Time[status == 2], pch = 3)
	points(ESP[status == 1], Time[status == 1], pch = 1)
        title("ESSP for males")


# 16.50 b) 

flung<-lung[lung$sex==2,]
nflung<-na.omit(flung[,c("time","status","ph.ecog","ph.karno",
"pat.karno","wt.loss")])
out <- coxph(Surv(time, status) 
~ ph.ecog + ph.karno + pat.karno + wt.loss,data = nflung)
        x <- as.matrix(nflung[,c(3,4,5,6)])
        ESP <- x %*% out$coef
        time <- as.vector(nflung[,1])
        status <- as.vector(nflung[,2])
        plot(ESP,time,pch = ' ')
        points(ESP[status == 2], time[status == 2], pch = 3)
	points(ESP[status == 1], time[status == 1], pch = 1)
        title("ESSP for Females")

#16.50 c)

vlung2(1)
title("males")

#16.50 d)

vlung2(2)
title("females")


# 16.51a)  Splus p. 269-275

library(survival)

nlung<-na.omit(lung[,c("time","status","sex","ph.ecog","ph.karno",
"pat.karno","wt.loss")])
par(mfrow=c(2,2))
attach(nlung)
fit1 <- coxph(Surv(time, status) 
~ strata(sex)+ ph.karno + pat.karno + wt.loss,data = nlung)
scatter.smooth(ph.ecog,resid(fit1))
fit2 <- coxph(Surv(time, status) 
~ strata(sex)+ ph.ecog + pat.karno + wt.loss,data = nlung)
scatter.smooth(ph.karno,resid(fit2))
fit3 <- coxph(Surv(time, status) 
~ strata(sex)+ ph.ecog + ph.karno + wt.loss,data = nlung)
scatter.smooth(pat.karno,resid(fit3))
fit4 <- coxph(Surv(time, status) 
~ strata(sex)+ ph.ecog + ph.karno + pat.karno,data = nlung)
scatter.smooth(wt.loss,resid(fit4))



# 16.51c)

lung.fit2 <- coxph(Surv(time, status) 
~ strata(sex)+ ph.ecog + ph.karno + pat.karno + wt.loss,data = nlung)

plot(cox.zph(lung.fit2))

cox.zph(lung.fit2)

par(mfrow=c(1,1))

#HW16.52

source("http://lagrange.math.siu.edu/Olive/regpack.txt")
library(survival)
par(mfrow=c(1,1))

#a)

phsim(k=1)

phsim(k=0.5)

phsim(k=5)

#b)

phsim2(n=100,k=1)

phsim2(n=200,k=1)

#c)

bphsim3(n=10,k=1)

#d)

phsim5(n=50,k=1)

phsim5(n=60,k=1)

#HW7 2)

source("http://lagrange.math.siu.edu/Olive/regdata.txt")
library(MASS)
library(survival)
alung<-as.data.frame(alung)
zc <- coxph(Surv(alung[,1],alung[,2])~perf+age+ttoent+size+type+ttype+trt,
data=alung)
stepAIC(zc) #default is backward

#HW 16.53 ovarian data Collett (2003, p. 335-6)

library(survival)

time <- c(156, 1040, 59, 421, 329, 769, 365, 770, 1227, 
268, 475, 1129, 464, 1206, 638, 563, 1106, 431, 855, 
803, 115, 744, 477, 448, 353, 377)

status <- c(1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1,
1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0)

treat <- c(1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,
0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0)

age <- c(66, 38, 72, 53, 43, 59, 64, 57, 59, 74, 59, 53,
56, 44, 56, 55, 44, 50, 43, 39, 74, 50, 64, 56, 63, 58)

ovar <- cbind(time,status,treat,age)

rm(time,status,treat,age)

z <- survreg(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4],dist="weibull")
summary(z)
zfit <- 11.548 -.561*ovar[,3] -.079*ovar[,4]
plot(zfit,log(ovar[,1]),pch = ' ')
points(zfit[ovar[,2]==0],log(ovar[ovar[,2]==0,1]),pch=1)
points(zfit[ovar[,2]==1],log(ovar[ovar[,2]==1,1]),pch=3)
abline(0,1) 
abline( lsfit( zfit[ovar[,2]==1],log(ovar[ovar[,2]==1,1]) )$coef)

# HW 16.55a ovarian data Collett (2003, p. 335-6)

library(survival)

time <- c(156, 1040, 59, 421, 329, 769, 365, 770, 1227, 
268, 475, 1129, 464, 1206, 638, 563, 1106, 431, 855, 
803, 115, 744, 477, 448, 353, 377)

status <- c(1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1,
1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0)

treat <- c(1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,
0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0)

age <- c(66, 38, 72, 53, 43, 59, 64, 57, 59, 74, 59, 53,
56, 44, 56, 55, 44, 50, 43, 39, 74, 50, 64, 56, 63, 58)

ovar <- cbind(time,status,treat,age)

rm(time,status,treat,age)

z <- survreg(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4],dist="weibull")
zx <- cbind(ovar[,3],ovar[,4])
zc <- coxph(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4])
ESPC <- zx%*%zc$coef
ESPW <- zx%*%z$coef[-1]
sighat<-z$scale
plot( - ESPW/sighat, ESPC)
abline(0,1)
title("a) Weibull")
cor( - ESPW/sighat, ESPC)

#HW 16.55b

z <- survreg(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4],dist="exponential")
ESPE <- zx%*%z$coef[-1]
shat<-z$scale
plot( - ESPE, ESPC)
abline(0,1)
title("b) Exponential")


#HW 16.59  lung data, Splus p. 269-275

library(survival)

sfit <-  coxph(Surv(time, status)~ strata(sex)+ ph.ecog  
+ pat.karno + wt.loss,data = lung,na.action=na.omit)

summary(sfit)

plot(survfit(sfit), lty = 2:3)
title("F Survival is Top Curve, Ave. Covariates")


#HW 16.62  (= HW11 3) R stanford heart data, generalzed Cox regression

library(survival)
options(contrast=c("contr.treatment","contr.poly"))
z <- coxph(Surv(start,stop,event)~transplant + surgery + age,data=heart)
z