R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

library(SparseM)
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
library(PolynomF)
library(modelfree)
library(psyphy)
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, round.POSIXt, trunc.POSIXt, units
pseData <- read.table("all_0727.txt", header=TRUE)
b <- pseData$h
r <- pseData$c
t <- pseData$f
m <- 168
p1 <- b/m
p2 <- r/m
p3 <- t/m

x <- c(750,833,900,983,1050,1133,1200,1283,1350)

d1 = data.frame(
  x, p1
  , sd1 = c(0.021735022,
            0.015368982,
            0.043470044,
            0.059125652,
            0.06547619,
            0.034193825,
            0.029761905,
            0.016835876,
            0.01497983
  )
)

d2 = data.frame(
  x, p2
  , sd2 = c(0.005952381,
            0.006873217,
            0.021735022,
            0.10304097,
            0.078064744,
            0.050507627,
            0.04051702,
            0.020331252,
            0.016835876
  )
)

d3 = data.frame(
  x, p3
  , sd3 = c(0.005952381,
            0.028129855,
            0.029761905,
            0.076150163,
            0.080154801,
            0.055306984,
            0.03504667,
            0.039333785,
            0.011397942
            
  )
)



#blue & red

#blue model
userlink<-weibull_link( 200 )
model <- glm(p1~x, family=binomial(link=userlink)) # fitting with lodistic 
## Warning: non-integer #successes in a binomial glm!
numxfit <- 999 # Number of points to be generated minus 1
xfit <- (max(x)-min(x)) * (0:numxfit) / numxfit + min(x)
lengthyfit <- predict(model, data.frame(x=xfit), type="response")
#red model
rmodel <- glm(p2~x, family=binomial(link=userlink)) # fitting with lodistic 
## Warning: non-integer #successes in a binomial glm!
numxfit <- 999 # Number of points to be generated minus 1
rxfit <- (max(x)-min(x)) * (0:numxfit) / numxfit + min(x)
rlengthyfit <- predict(rmodel, data.frame(x=rxfit), type="response")

#control 2.077
cmodel <- glm(p3~x, family=binomial(link=userlink)) # fitting with lodistic 
## Warning: non-integer #successes in a binomial glm!
numxfit <- 999 # Number of points to be generated minus 1
cfit <- (max(x)-min(x)) * (0:numxfit) / numxfit + min(x)
clengthyfit <- predict(cmodel, data.frame(x=cfit), type="response")

plot(p1~x, xlab="o", ylab="p",  xlim = c(750, 1350),ylim = c(0,1.0),axes=FALSE,pch=16, cex=3)#,col="gray53",)
with (data = d1, expr = errbar(x, p1, p1+sd1, p1-sd1, add=T, pch=10, cap=.1))
lines(xfit, lengthyfit, lwd=5, lty=4)#, cex.axis=1.5, col="gray53",)

par(new=TRUE,cex.axis=1.5)

plot(p2~x,xlab="o", ylab="p", ylim= c(0,1.0),axes=FALSE,pch="o", cex=3)#, col="chartreuse3")
with (data = d2, expr = errbar(x, p2, p2+sd2, p2-sd2, add=T, pch=10, cap=.1))
lines(rxfit, rlengthyfit, lwd=5)#col="chartreuse3",)

par(new=TRUE,cex.axis=1.5)

plot(p3~x,xlab="o", ylab="p", ylim= c(0,1.0),axes=FALSE,pch=17, cex=3)#, col="green"
with (data = d3, expr = errbar(x, p3, p3+sd3, p3-sd3, add=T, pch=10, cap=.1))
lines(cfit, clengthyfit,lwd=5, lty=5)#, col="green"

axis(1, at=c(750,1050,1350),las=0)
axis(2, at=c(0,0.5,0.75,1.0),las=0)

#750,833,900,983,1050,1133,1200,1283,1350

OpseResult<-threshold_slope(lengthyfit, xfit, 0.75)
pse1<-threshold_slope(lengthyfit, xfit, 0.5)
pse <- pse1$x_th
OpseResult
## $x_th
## [1] 1128.979
## 
## $slope
## [1] 0.002384068
pse
## [1] 1028.078
rOpseResult<-threshold_slope(rlengthyfit, rxfit, 0.75)
rpse<-threshold_slope(rlengthyfit, rxfit, 0.5)
rpse <- rpse$x_th
rOpseResult
## $x_th
## [1] 1150.601
## 
## $slope
## [1] 0.002533996
rpse
## [1] 1056.306
cOpseResult<-threshold_slope(clengthyfit, cfit, 0.75)
cpse<-threshold_slope(clengthyfit, cfit, 0.5)
cpse <- cpse$x_th
cOpseResult
## $x_th
## [1] 1146.997
## 
## $slope
## [1] 0.002562801
cpse
## [1] 1053.303

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.