f.plot=function(x,y,nr=20,nc=20, scale="raw",...) {
zx = c(1:nr,rep(1,nc),1+trunc( nr*(x- min(x))/(max(x)-min(x)) ))
zx[zx>nr] = nr
zy = c(rep(1,nr),1:nc,1+trunc( nc*(y- min(y))/(max(y)-min(y)) ))
zy[zy>nc] = nc
z = table(zx,zy); z[,1]=z[,1]-1; z[1,]=z[1,]-1;
if (scale=="l") z= log(1+z)
#image(z=t(z),x=seq(length=nr+1,from=min(x),to=max(x)),
#y= seq(length=nc+1,from=min(y),to=max(y)),
#xlab="",ylab="", col=rainbow(1000),...)
contour(t(z),x=seq(length=nr,from=min(x),to=max(x)),
y= seq(length=nc,from=min(y),to=max(y)),...)
}
x = rnorm(10000) ; y = rnorm(10000)
ux = rnorm(5000)/3
uy = ux^2 -0.5
plot(c(y,uy)+20,c(x,ux),xlab="x",ylab="y",xlim=c(18,22),ylim=c(-2,2),pch=16)
f.plot(c(y,uy)+20,c(x,ux),50,50,'l',xlim=c(18,23),ylim=c(-3.2,2.6))
f.plot(c(y,uy)+20,c(x,ux),50,50,'l',xlim=c(18,22),ylim=c(-3.2,2.6),nlevels=4)
plot(c(y,uy)+20,c(x,ux),xlab="x",ylab="y",xlim=c(18,22),ylim=c(-2,2),pch=16)
save.image("C:\\Documents and Settings\\EmirB\\Desktop\\LDesktop\\Columbia\\Computational Statistics Course\\Lecture\\.RData")
f.plot(c(y,uy)+20,c(x,ux),50,50,'l',xlim=c(18,22),ylim=c(-3.2,2.6),nlevels=4)
library(Rcmdr)
ls()
names(eff2)
plot(eff2$PAIN0,eff2$PNCHG)
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base')
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots")
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test")
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10))
abline(h=0,lty=2,col=3,lwd=2)
abline(h=0,lty=2,col="lightgrey",lwd=2)
abline(h=0,lty=2,col="lightblue",lwd=2)
?abline
abline(a=1,b=2,lty=1,col="red")
abline(a=1,b=1,lty=1,col="red")
lines(c(4,10),c(-10,4))
pdf("mygraph.pdf")
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10))
lines(c(4,10),c(-10,4))
abline(a=1,b=1,lty=1,col="red")
dev.off()
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10))
dev.off()
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10))
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10),pch=9))
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10),pch=9)
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10),pch=8)
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',main="Scatterplots \n My test",xlim=c(4,10),pch=".")
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",ylab='Change from base',col.lab=2,main="Scatterplots \n My test",xlim=c(4,10),pch=".")
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",
ylab='Change from base',col.lab=2,main="Scatterplots \n My test",
xlim=c(4,10),pch=".",font=4,fon.lab=2)
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",
ylab='Change from base',col.lab=2,main="Scatterplots \n My test",
xlim=c(4,10),pch=".",font=4,font.lab=2)
# specify the data x <- c(1:10); y <- x; z <- 10/x# create extra margin room on the right for an axis par(mar=c(5, 4, 4, 8) + 0.1)# plot x vs. y plot(x, y,type="b", pch=21, col="red",    yaxt="n", lty=3, xlab="", ylab="")# add x vs. 1/x lines(x, z, type="b", pch=22, col="blue", lty=2)# draw an axis on the left axis(2, at=x,labels=x, col.axis="red", las=2)# draw an axis on the right, with smaller text and ticks axis(4, at=z,labels=round(z,digits=2),  col.axis="blue", las=2, cex.axis=0.7, tck=-.01)# add a title for the right axis mtext("y=1/x", side=4, line=3, cex.lab=1,las=2, col="blue")# add a main title and bottom and left axis labels title("An Example of Creative Axes", xlab="X values",   ylab="Y=X") 
y
x <- c(1:10); y <- x; z <- 10/x
x
y
z
par(mar=c(5, 4, 4, 8) + 0.1)
plot(x, y,type="b", pch=21, col="red", 
   yaxt="n", lty=3, xlab="", ylab="")
lines(x, z, type="b", pch=22, col="blue", lty=2)
axis(2, at=x,labels=x, col.axis="red", las=2)
?AXIS
?axis
attach(mtcars)par(mfrow=c(2,2))plot(wt,mpg, main="Scatterplot of wt vs. mpg")plot(wt,disp, main="Scatterplot of wt vs disp")hist(wt, main="Histogram of wt")boxplot(wt, main="Boxplot of wt")
par(mfrow=c(2,2))
plot(wt,mpg, main="Scatterplot of wt vs. mpg")
plot(wt,disp, main="Scatterplot of wt vs disp")
hist(wt, main="Histogram of wt")
boxplot(wt, main="Boxplot of wt")
attach(mtcars)
detach()
data(mtcars)
ls()
summary(mtcars)
mtcasr
mtcars
attach(mtcars)
par(mfrow=c(2,2))
plot(wt,mpg, main="Scatterplot of wt vs. mpg")
plot(wt,disp, main="Scatterplot of wt vs disp")
hist(wt, main="Histogram of wt")
boxplot(wt, main="Boxplot of wt")
par(mfrow=c(3,1)) hist(wt)hist(mpg)hist(disp) 
hist(wt)
hist(mpg)
hist(disp) 
par(mfrow=c(3,1))
hist(wt)
hist(mpg)
hist(disp) 
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE), 
   widths=c(3,1), heights=c(1,2))
hist(wt)
hist(mpg)
hist(disp)
library(lattice)
splom(mtcars[c(1,3,5,6)], groups=cyl, data=mtcars,
   panel=panel.superpose, 
   key=list(title="Three Cylinder Options",
   columns=3,
   points=list(pch=super.sym$pch[1:3],
   col=super.sym$col[1:3]),
   text=list(c("4 Cylinder","6 Cylinder","8 Cylinder")))) 
data(mtcars)
library(lattice)
splom(mtcars[c(1,3,5,6)], groups=cyl, data=mtcars,
   panel=panel.superpose, 
   key=list(title="Three Cylinder Options",
   columns=3,
   points=list(pch=super.sym$pch[1:3],
   col=super.sym$col[1:3]),
   text=list(c("4 Cylinder","6 Cylinder","8 Cylinder")))) 
super.sym$col[1:3]
splom(mtcars[c(1,3,5,6)], groups=cyl, data=mtcars,
   panel=panel.superpose, 
   key=list(title="Three Cylinder Options",
   columns=3,
   points=list(pch=mtcars$pch[1:3],
   col= mtcars$col[1:3]),
   text=list(c("4 Cylinder","6 Cylinder","8 Cylinder")))) 
splom(mtcars[c(1,3,5,6)], groups=cyl, data=mtcars,
   panel=panel.superpose, 
   key=list(title="Three Cylinder Options",
   columns=3,
   text=list(c("4 Cylinder","6 Cylinder","8 Cylinder")))) 
splom(mtcars[c(1,3,5,6)], groups=cyl, data=mtcars,
   panel=panel.superpose, 
   key=list(title="Three Cylinder Options",
   columns=3,
   text=list(c("4 Cylinder","6 Cylinder","8 Cylinder")))) 
splom(mtcars[c(1,3,5,6)], groups=mtcars$cyl, data=mtcars,
   panel=panel.superpose, 
   key=list(title="Three Cylinder Options",
   columns=3,
   text=list(c("4 Cylinder","6 Cylinder","8 Cylinder")))) 
library(car)
scatterplot.matrix(~mpg+disp+drat+wt|cyl, data=mtcars,
   main="Three Cylinder Options")
library(hexbin)x <- rnorm(1000)y <- rnorm(1000)bin<-hexbin(x, y, xbins=50) plot(bin, main="Hexagonal Binning") 
install.package("hexbin")
help.search("install")
utils:::menuInstallPkgs()
utils:::menuInstallPkgs()
utils:::menuInstallPkgs("hexbin")
utils:::menuInstallPkgs(hexbin)
utils:::menuInstallPkgs("USA(IA)")
?utils:::menuInstallPkgs
??utils:::menuInstallPkgs
utils:::menuInstallPkgs()
utils:::menuInstallPkgs('http://streaming.stat.iastate.edu/CRAN/bin/windows/contrib/2.9/hexbin.zip')
utils:::menuInstallPkgs('http://streaming.stat.iastate.edu/CRAN/bin/windows/contrib/2.9/')
utils:::menuInstallPkgs('http://streaming.stat.iastate.edu/')
install = function() {
install.packages(c("moments", “hexbin“, “Rcmdr“)
repos="http://lib.stat.cmu.edu/R/CRAN")
}
install()
install.packages(c("moments", "hexbin", "Rcmdr")
install = function() {
install.packages(c("moments", “hexbin“, “Rcmdr“)
install = function() {
install.packages(c("moments", "hexbin", "Rcmdr")
repos="http://lib.stat.cmu.edu/R/CRAN")
install = function() {
install.packages(c("moments","randomForest","nnet","e1071","gbm","lars","lasso2","boot","cluster","graphics","mclust","ARF"),
repos="http://lib.stat.cmu.edu/R/CRAN")
}
install()
install = function() {
install.packages(c("moments","graphics","Rcmdr","hexbin"),
repos="http://lib.stat.cmu.edu/R/CRAN")
}
install()
library(hexbin)x <- rnorm(1000)y <- rnorm(1000)bin<-hexbin(x, y, xbins=50) plot(bin, main="Hexagonal Binning") 
x <- rnorm(1000)
y <- rnorm(1000)
bin<-hexbin(x, y, xbins=50) 
plot(bin, main="Hexagonal Binning") 
library(hexbin)
bin<-hexbin(x, y, xbins=50) 
plot(bin, main="Hexagonal Binning") 
with(mydata, t.test(y1,y2))
with(mtcars, t.test())
names(mydata)
names(mtcars)
with(mtcars, t.test(mpg,hp))
install = function() {
install.packages(c("pwr"),
repos="http://lib.stat.cmu.edu/R/CRAN")
}
install()
library(pwr)# range of correlationsr <- seq(.1,.5,.01)nr <- length(r)# power valuesp <- seq(.4,.9,.1)np <- length(p)# obtain sample sizessamsize <- array(numeric(nr*np), dim=c(nr,np))for (i in 1:np){  for (j in 1:nr){    result <- pwr.r.test(n = NULL, r = r[j],    sig.level = .05, power = p[i],    alternative = "two.sided")    samsize[j,i] <- ceiling(result$n)  }}
library(pwr)
# range of correlations
r <- seq(.1,.5,.01)
nr <- length(r)
# power values
p <- seq(.4,.9,.1)
np <- length(p)
# obtain sample sizes
samsize <- array(numeric(nr*np), dim=c(nr,np))
for (i in 1:np){
  for (j in 1:nr){
    result <- pwr.r.test(n = NULL, r = r[j],
    sig.level = .05, power = p[i],
    alternative = "two.sided")
    samsize[j,i] <- ceiling(result$n)
  }
}
xrange <- range(r)
yrange <- round(range(samsize))
colors <- rainbow(length(p))
plot(xrange, yrange, type="n",
  xlab="Correlation Coefficient (r)",
  ylab="Sample Size (n)" )
for (i in 1:np){
  lines(r, samsize[,i], type="l", lwd=2, col=colors[i])
}
abline(v=0, h=seq(0,yrange[2],50), lty=2, col="grey89")
abline(h=0, v=seq(xrange[1],xrange[2],.02), lty=2,
   col="grey89")
title("Sample Size Estimation for Correlation Studies\n
  Sig=0.05 (Two-tailed)")
legend("topright", title="Power", as.character(p),
   fill=colors)
ls()
names(eff2)
lr = glm( PNRSP30 ~ TRT , data= eff2, family=binomial())
lr
summary(lr)
save.image("C:\\Documents and Settings\\EmirB\\Desktop\\LDesktop\\Columbia\\Computational Statistics Course\\Lecture\\.RData")
growth=read.csv("./GrowthCurve.csv",T)
growth=read.csv("./growthCurve.csv",T)
growth=read.csv("./Datasets/growthCurve.csv",T)
growth
plot(growth)
plot(growth,type="b"))
plot(growth,type="b"))
plot(growth,type="b")
fit=nlr(y~a-b*(g^x),start = c(ed50=25,e0=5,emax=20), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=edd)
summary(nls.fit)
fit=nlr(y~a-b*(g^x),start = c(a=2,b=1,g=.8), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=growth)
?nlr
help.search("nlr")
fit=nls(y~a-b*(g^x),start = c(a=2,b=1,g=.8), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=growth)
summary(nls)
summary(fit)
plot(fit)
history(20)
plot(growth)
ls()
edd=read.table("./DataSets/edd.txt",T)
ffit = function(doses =c(0,100), ed50=75.76,e0=1.885,emax=15.85,lambda=0.75,pl=T,...) {
dose = seq(doses[1],doses[2],length=100)
yy=e0 +(emax * dose^lambda)/(dose^lambda+ed50^lambda)
if (pl) plot(dose,yy,type="l",...) else lines(dose,yy,...)
}
i = sort.list(edd$dose)
doses = edd$dose[i]
y = edd$ch8[i]
u =  split(y,doses)
sd1 = sapply(u,sd)
mu1 = sapply(u,mean)
n1 = sapply(u,length)
cls = mu1+ cbind(-sd1,sd1)*qt(0.95,n1-1)/sqrt(n1) 
# EMAX 3 and 4 param  plot
ffit(doses=c(0,15),ed50=5,e0=2.0984807,emax=10,lambda=.5,lty=2,pl=T,ylim=c(0,15),ylab="Change in EDD",xlab="Dose (mg)",col=2)
fit=nls(y~a-b*(g^x),start = c(a=2,b=1,g=.8), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=growth)
summary(fit)
fit=nls(y~a-b*(g^x),start = c(a=1,b=1,g=.5), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=growth)
?nls
fit=nls(y~a-b*(g^x),start = c(a=1,b=1,g=.5), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=growth,lower=(g >0 & g < 1))
fit=nls(y~a-b*(g^x),start = c(a=1,b=1,g=.5), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=growth,alg="port",lower=c( g >0 & g < 1))
maxiter = 100),trace=T,na.action=na.omit,data=growth,alg="port",lower=c(g >0),upper=c(g < 1))
maxiter = 100),trace=T,na.action=na.omit,data=growth,alg="port",lower=c(g >0),upper=c(g < 1))
ffit()
ffit(doses=c(0,75),ed50=5,e0=1,emax=10,lambda=1,lty=1,pl=T,ylim=c(0,10),main="Lambda <= 1",ylab="Response",xlab="Dose (mg)",col=1,lwd=2)
ffit(doses=c(0,75),ed50=5,e0=1,emax=10,lambda=.75,lty=2,pl=F,ylim=c(0,10),main="",ylab="",xlab="Dose (mg)",col=2,lwd=2)
ffit(doses=c(0,75),ed50=5,e0=1,emax=10,lambda=1,lty=1,pl=T,ylim=c(0,10),main="Lambda <= 1",ylab="Response",xlab="Dose (mg)",col=1,lwd=2)
ffit(doses=c(0,75),ed50=5,e0=1,emax=10,lambda=.75,lty=2,pl=F,ylim=c(0,10),main="",ylab="",xlab="Dose (mg)",col=2,lwd=2)
ffit(doses=c(0,75),ed50=5,e0=1,emax=10,lambda=.5,lty=3,pl=F,ylim=c(0,10),main="",ylab="",xlab="Dose (mg)",col=3,lwd=2)
history(1000)
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",
ylab='Change from base',col.lab=2,main="Scatterplots \n My test",
xlim=c(4,10),pch=".",font=4,fon.lab=2)
plot(eff2$PAIN0,eff2$PNCHG,xlab="Baseline Pain",
ylab='Change from base',col.lab=2,main="Scatterplots \n My test",
xlim=c(4,10),pch=".",font=4,font.lab=2)
library(Rcmdr)
summary(eff2)
eff2$TRT=as.factor(eff2$TRT)
eff2$CSEX=as.factor(eff2$CSEX)
summary(eff2)
table(eff$TRT)
table(eff2$TRT)
AnovaModel.1 <- aov(PNCHG ~ TRT, data=eff2)
aov(PNCHG ~ TRT, data=eff2)
summary(AnovaModel.1)
numSummary(eff2$PNCHG , groups=eff2$TRT, statistics=c("mean", "sd"))
?aov
?lm
edd
names(edd)
edd[1:10,]
plot(edd$dose,edd$ch8)
table(edd$dose)
plotMeans(edd$ch8, edd$dose, error.bars="conf.int", level=0.95)
plotMeans(edd$ch8, edd$dose)
plotMeans(edd$ch8, as.factor(edd$dose))
?nls
nls.fit<-nls(ch8 ~ e0 + (emax * dose)/(dose + ed50), 
start = c(ed50=25,e0=5,emax=20), control = nls.control(
maxiter = 100),trace=T,na.action=na.omit,data=edd)
summary(nls.fit)
ffit = function(doses =c(0,100), ed50=75.76,e0=1.885,emax=15.85,lambda=0.75,pl=T,...) {
dose = seq(doses[1],doses[2],length=100)
yy=e0 +(emax * dose^lambda)/(dose^lambda+ed50^lambda)
if (pl) plot(dose,yy,type="l",...) else lines(dose,yy,...)
}
ffit()
plot(edd$dose,edd$ch8)
par(new=T)
ffit()
ffit(doses=c(0,75),ed50=5,e0=1,emax=10,lambda=1,lty=1,pl=T,ylim=c(0,10),main="",ylab="Response",xlab="Dose (mg)",col=1,lwd=2)
example(plot)
save.image("C:\\Documents and Settings\\EmirB\\Desktop\\LDesktop\\Columbia\\Computational Statistics Course\\Lecture\\Jan14.RData")