setwd("C:/cwru/classes/epbi473_05/lec1(overView.introR)")
#1.4
search()
ls()
data(thuesen)
library(ISwR)
search()
data(thuesen)
ls()
thuesen


blood.glucose
thuesen$blood.glucose
search()
ls()
attach(thuesen)
ls()
blood.glucose
detach(thuesen)


#1.7
m<-matrix(rnorm(200),10)
apply(m,2,mean)
apply(m,1,mean)
mean(m)
sin(m)
s=apply(m,1,sum);s
qqnorm(s)


#2.1
dnorm(-5:5)
pnorm(-5:5)
plot(-5:5,pnorm(-5:5),type="l")
1-pnorm(42,mean=35,sd=6)
#2.2
qnorm(.95)
qnorm(.975)

x<-seq(-4,4,0.1)
F<-seq(0,1,0.01)
par(mfrow=c(2,2))
plot(x,dnorm(x),type='l')
plot(x,pnorm(x),type='l')
plot(F,qnorm(F),type='l')
plot(rnorm(10))
par(mfrow=c(1,1))


# chapter 4
# 4.1
data(react)
hist(react)
t.test(react,mu=0)
qqnorm(react)
wilcox.test(react,mu=0)
summary(lm1<-lm(react~1))

model.matrix(lm1)


# 5.1 
data(rmr)
attach(rmr)
lm1<-lm(metabolic.rate~body.weight);lm1
class(lm1)
str(lm1)
plot(body.weight,metabolic.rate)
abline(lm1)
legend(locator(n=1),legend=c('Line','Data'),lty=c(1,NA),pch=c(NA,1))
coef(lm1)
lm1$coef
smet<-summary(lm1)
smet$coef
tsd=lm1$coef[2,2]*qt(0.975,length(metabolic.rate)-2)
ci=c(lm1$coef[2,1]-tsd,lm1$coef[2,1]+tsd)
ci

detach(rmr)
plot(lm1)
model.matrix(lm1)
detach(rmr)
?rmr



plot(1:30,pch=1:30)
?plot
#

# 5.3
?malaria
data(malaria)
attach(malaria)
malaria
y=log(ab)
summary(a<-lm(y~age))
plot(y~age)
abline(a)
plot(fitted(a),resid(a))
qqnorm(resid(a))



# 6.2
data(lung)
lung
attach(lung)
anova(lm0<-lm(volume~method+subject))
model.matrix(lm0)
interaction.plot(method,subject,volume)
plot(volume~method)
anova(lm1<-lm(volume~method))
model.matrix(lm1)
detach(lung)


# 9.1
data(cystfibr)
cystfibr
attach(cystfibr)
pairs(cystfibr,gap=0)
# 9.3
summary(lm1<-lm(pemax~age+sex+height+weight+bmp+fev1+rv+frc+tlc))
step(lm1)
summary(lm1<-lm(pemax~weight+bmp+fev1+rv))
summary(lm1<-lm(pemax~weight+bmp+fev1))
# 10.1
summary(lm1<-lm(pemax~height+I(height^2)) )
pred=data.frame(height=seq(110,180,2))
pp=predict(lm1,interval="pred",newdata=pred);pp
pc=predict(lm1,interval="conf",newdata=pred);pc
plot(height,pemax,ylim=c(0,200)) # plots the data
matlines(pred$height,pp,lty=c(1,2,2),col="black") # CI for new data point
matlines(pred$height,pc,lty=c(1,3,3),col="black") # CI for mean value
detach(cystfibr)


#
# 10.3
#
data(red.cell.folate)
attach(red.cell.folate)
red.cell.folate
#?red.cell.folate
getOption("contrasts")
model.matrix(lm1<-lm(folate~ventilation))
summary(lm1)

options(contrasts=c("contr.sum", "contr.poly"))
getOption("contrasts")
model.matrix(lm1<-lm(folate~ventilation))
summary(lm1)

options(contrasts=c("contr.treatment", "contr.poly"))  #set back to default
detach(red.cell.folate)
#


#
# 10.4
#
data(fake.trypsin)
attach(fake.trypsin)
fake.trypsin
summary(fake.trypsin)
class(fake.trypsin)
lapply(fake.trypsin,class)
summary(lm1f<-lm(trypsin~grpf))
summary(lm1<-lm(trypsin~grp))
anova(lm1,lm1f)
m=tapply(trypsin,grpf,mean);m
stripchart(trypsin~grp,"jitter",jitter=.1,vertical=T,pch=20)
lines(1:6,m,type="b",pch=4,lty=2,cex=2)
abline(lm(trypsin~grp))
detach(fake.trypsin)


#
# 11
#
library(ISwR)
data(graft.vs.host)
graft.vs.host
?graft.vs.host
plot(graft.vs.host,gap=0)
attach(graft.vs.host)
t.test(index~gvhd) # output treated as input, should be other way around
summary(glm1<-glm(gvhd~index,family=binomial)) # logistic regression
t.test(time~gvhd) 
kruskal.test(time~gvhd) 
summary(glm1<-glm(gvhd~index+donage+rcpage+type+preg,family=binomial))
step(glm1)
summary(glm1<-glm(gvhd~index+donage+type+preg,family=binomial))
summary(glm1<-glm(gvhd~index+donage+preg,family=binomial))
summary(glm1<-glm(gvhd~index+donage,family=binomial))
#
# 12
#
library(survival)
sv=Surv(time,dead);sv
str(sv)
class(sv)
plot(sv)
summary(sf<-survfit(sv~gvhd)) 
sf
str(sf)
class(sf)
plot(sf)
sdiff<-survdiff(sv~gvhd);sdiff
class(sdiff)
scph<-coxph(sv~gvhd);scph
ftype=factor(type,labels=c("AML","ALL","CML"));ftype
scph<-coxph(sv~gvhd+index+donage+rcpage+ftype+preg);scph
step(scph)
class(scph)


library(ISwR)
library(survival)
data(melanom)
attach(melanom)
melanom
summary(melanom)
?melanom
sv=Surv(days,status==1)

par(mfrow=c(2,3))
summary(sf<-survfit(sv~sex)) 
plot(sf,main="sex")
summary(sf<-survfit(sv~ulc)) 
plot(sf,main="ulceration")
m=median(thick)
fthick=cut(thick,breaks=c(0,m,10000));fthick
summary(sf<-survfit(sv~fthick)) 
plot(sf,main="thickness")
summary(sf<-survfit(sv~sex+ulc)) 
plot(sf,main="sex+ulc")
summary(sf<-survfit(sv~sex+fthick)) 
plot(sf,main="sex+fthick")
summary(sf<-survfit(sv~sex+fthick+ulc)) 
plot(sf,main="sex+fthick+ulc")
par(mfcol=c(1,1))

coxph(sv~sex)
coxph(sv~strata(ulc)+sex)
coxph(sv~strata(thick)+sex)
coxph(sv~strata(ulc)+strata(thick)+sex)

detach(melanom)


setwd("C:/cwru/classes/epbi473_05/lec1(overView.introR)")
dat=read.table("melanom.txt",header=T);dat # get this from the ISwR installation Data directory
#                                             e.g. C:\Program Files\R\rw2010\library\ISwR\data
dat[1:5,]=0
dat[1:10,]  # dump top 10 rows
write.table(dat,"newMelanom.txt",sep='\t',row.names=F)
ndat=read.table("newMelanom.txt",header=T)
ndat[1:10,]