library(epicalc) data(package="epicalc",Outbreak) attach(Outbreak) plot(vomiting~age,main="not linear!",col="purple") model.1<-lm(vomiting~age) summary(model.1) par(mfrow=c(2,2)) plot(model.1,which=1:4) sum.1<-summary(glm(vomiting~age,family=binomial(link=logit))) 1-pchisq(1452.3-1433.9,1093-1092) #simulated x and y x<-c(seq(0,100,0.5),sample(seq(0,100,0.5),10,replace=T)) y<-c() p<-c() for(i in 1:length(x)){ p[i]<-x[i]/60+rnorm(1,0,0.1) y[i]<-rbinom(1, 1, ifelse(p[i]>1,0.99,ifelse(p[i]<0,0.01,p[i])) ) } #above code generates random 0/1 outcomes dependent (a little) on the x inputs. Let's look at what we generated: plot(x,y) #now run a logistic model on these data: random.1<-glm(y~x,family=binomial(link=logit)) summary(random.1) beta0.plus.beta1.x <- random.1\$coefficients[1]+ random.1\$coefficients[2]*x #above is the model-fitted value for the log(odds) p <- exp(beta0.plus.beta1.x)/(1+exp(beta0.plus.beta1.x)) #p = estimated probability of event hapenning (getting sick, dying, etc.), given the predictor(s) X, is found through some algebra plot(x,beta0.plus.beta1.x,ylim=c(-.2,1.3),main="black is log(odds)=b0+b1X,\n red is p = b0+b1X, \n blue/purple are observed (0/1)",ylab="") points(x, y,col=ifelse(y==0,"blue","purple")) points(x,p,col="red") rm(list = ls()) library(survival) data(leukemia) detach(Outbreak) data(package="survival",leukemia) attach(leukemia) hist(time[x=="Maintained"]) hist(time[x=="Nonmaintained"]) hist(time, xlab="Length of Survival Time", main="Histogram of Survial Time in AML Paitents") survfit(Surv(time, status)~1) leukemia[order(time),] summary(survfit(Surv(time, status)~1)) surv.aml <- survfit(Surv(time,status)~1) summary(surv.aml) plot(surv.aml, main = "Plot of Survival Curve for AML Patients", xlab = "Length of Survival", ylab = "Proportion of Individuals who have Survived") print(survfit(Surv(time,status)~x)) plot(survfit(Surv(time,status)~x), main = "Plot of Survival Curves by Chemotherapy Group", xlab = "Length of Survival (t)",ylab="Proportion of Individuals who have Survived",col=c("blue","red")) legend("topright", legend=c("Maintained", "Nonmaintained"),fill=c("blue","red"),bty="n") survdiff(Surv(time,status)~x) detach(leukemia) attach(cancer) print(survfit(Surv(time,status)~sex)) plot(survfit(Surv(time,status)~sex), main = "Plot of Survival Curves by Sex", xlab = "Length of Survival (t)",ylab="Proportion of Individuals who have Survived",col=c("blue","red")) legend("topright", legend=c("Male", "Female"),fill=c("blue","red"),bty="n") coxph(Surv(time,status)~sex) coxph(formula = Surv(time, status) ~ sex + age)