hdr
chplot
library(lattice)#
pt<-data.frame(x=runif(15000,-1,1),y=runif(15000,-.5,.5))#
trend<-data.frame(x=seq(-1,1,len=1000),y=-seq(-1,1,len=1000)^2+.5 +#
runif(1000,-0.05,0.05))#
xyplot(y~x,rbind(pt,trend),pch=16,col='black')#
xyplot(y~x,rbind(pt,trend),pch=16,col='black',alpha=.2)#
xyplot(y~x,rbind(pt,trend),pch=16,col='black',cex=.2)#
#
Convex Hull#
dat1<-data.frame(x=rnorm(100,sd=.1),y=rnorm(100,sd=.1))#
dat2<-data.frame(x=rnorm(100),y=rnorm(100))#
chullpol<-function(dat){#
 ch<-chull(dat)#
 ch<-c(ch,ch[1])#
 dat[ch,]#
}#
plot(dat2,type='n')#
segments(dat1$x,dat1$y,dat2$x,dat2$y,col=grey(.85))#
polygon(chullpol(dat1),border="blue")#
polygon(chullpol(dat2),border="orange")
library(lattice)#
pt<-data.frame(x=runif(15000,-1,1),y=runif(15000,-.5,.5))#
trend<-data.frame(x=seq(-1,1,len=1000),y=-seq(-1,1,len=1000)^2+.5 +#
runif(1000,-0.05,0.05))#
xyplot(y~x,rbind(pt,trend),pch=16,col='black')
plot(dat2,type='n')
segments(dat1$x,dat1$y,dat2$x,dat2$y,col=grey(.85))
pt<-data.frame(x=runif(15000,-1,1),y=runif(15000,-.5,.5))#
trend<-data.frame(x=seq(-1,1,len=1000),y=-seq(-1,1,len=1000)^2+.5 +#
runif(1000,-0.05,0.05))#
xyplot(y~x,rbind(pt,trend),pch=16,col='black')
xyplot(y~x,rbind(pt,trend),pch=16,col='black',cex=.2)
xyplot(y~x,rbind(pt,trend),pch=16,col='black',alpha=.2)
dat1<-data.frame(x=rnorm(100,sd=.1),y=rnorm(100,sd=.1))#
dat2<-data.frame(x=rnorm(100),y=rnorm(100))#
chullpol<-function(dat){#
 ch<-chull(dat)#
 ch<-c(ch,ch[1])#
 dat[ch,]#
}
plot(dat2,type='n')#
segments(dat1$x,dat1$y,dat2$x,dat2$y,col=grey(.85))
polygon(chullpol(dat1),border="blue")
polygon(chullpol(dat2),border="orange")
install.packages("RLRsim")
library(RLRsim)
library(help=RLRsim)
?LRTsim
?extract.lmeDesign
library(lme4)
?RLRsim-package
?RLTSim
?LRTSim
g <- rep(1:10, e = 10)#
x <- rnorm(100)#
y <- 0.1 * x + rnorm(100)#
m <- lmer(y ~ x + (1|g), method="ML")#
m0 <- lm(y ~ 1)
g <- rep(1:10, e = 10)#
x <- rnorm(100)#
y <- 0.1 * x + rnorm(100)#
m <- lmer(y ~ x + (1|g), REML=FALSE)#
m0 <- lm(y ~ 1)
y
summary(y)
mean(y)
std(y)
sd(y)
g
table(g)
(obs.LRT <- 2*(logLik(m)-logLik(m0)))#
X <- m@X#
Z <- t(as.matrix(m@Zt))#
sim.LRT <- LRTSim(X, Z, 1, diag(10))#
(pval <- mean(sim.LRT > obs.LRT))
summary(m0)
print(m1)
print(m0)
print(m)
fm1 <- lmer(Reaction ~ Days + (1 | Days), sleepstudy)
exactRLT(fm1)
fm2 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
print(fm2)
exactLRT(fm2)
fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
exactRLRT(fm1)
anova(fm1, fm2)
fm3 <- lmer(Reaction ~ poly(Days,2) + (poly(Days,2) | Subject), sleepstudy)
exactRLRT(fm3)
mA <- lmer(Reaction ~ I(Days-4.5) + (1|Subject) + (0 + I(Days-4.5)|Subject), sleepstudy)#
m0 <- update(mA, . ~ . - (0 + I(Days-4.5)|Subject))#
m.slope  <- update(mA, . ~ . - (1|Subject))
m.slope
mA <- lmer(Reaction ~ I(Days-4.5) + (1|Subject) + (0 + I(Days-4.5)|Subject), sleepstudy)
m0 <- update(mA, . ~ . - (0 + I(Days-4.5)|Subject))
m1
ma
mA
m0
exactRLRT(m.slope, mA, m0)
anova(m0, mA)
fm2
anova(fm1,fm2)
fm1
exactRLTsim(fm2)
exactRLRTsim(fm2)
exactRLRT(fm2)
library(ggplot2)
cxc <- ggplot(mtcars, aes(x = factor(cyl))) +#
 geom_bar(width = 1, colour = "black")
cxc
cxc + geom_hline(yintercept = 5) + coord_polar()
cxc + geom_hline(yintercept = 5)
cxc + geom_segment(x = 0.5, xend = 3.5, y = 5, yend = 5) + coord_polar()
library(ggplot2)#
library(quantreg)#
library(Hmisc)#
#
rm(list=ls())#
#
# http://had.co.nz/ggplot2/stat_summary.html#
stat_sum_single <- function(fun, geom="point", ...) { #
   stat_summary(fun=fun, colour="blue", geom=geom, size = 2, ...) #
}#
#
stat_sum_df <- function(fun, geom="smooth", ...) {#
  stat_summary(fun=fun, colour="red", geom=geom, method=lm, width=0.2, size=1,...)#
}
load("PSC275_fpf.sfc_R.rda")#
#
d$ldur = log(d$dur)#
#
tmp <- by(d, d$id, function(x) rq(ldur ~ f + f1 + f2 + p + p1 + p2 + l + l1 + l2#
                                       + lsr + o + I(o^2) + ao + I(f^2) + I(f^3) + I(p^2) + I(p^3),#
                                       tau=seq(10,90,10)/100,data=x))#
# Frequency and predictability#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.05, to=+0.02, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.05, +0.02) )#
fp.plot <- p + geom_hline()
fp.plot
qr.1 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao#
                                        + I(f:f1) + I(f:l)   + I(f2:l) + I(p2:l), #
                                       tau=seq(10,90,10)/100,data=x))#
#
# 14-parameter main effects only model (Kliegl et al., 2006; using (1/launch site) and 1/ao)#
qr.2 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao, #
                                       tau=seq(10,90,10)/100,data=x))
qr.1 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao#
                                        + I(f*f1)  + I(f*l)   + I(f2*l) + I(p2*l), #
                                       tau=seq(10,90,10)/100,data=x))
tmp <- qr.1#
#
# Frequency and predictability#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.05, to=+0.02, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.05, +0.03) )#
rq1.fp.plot <- p + geom_hline()
rq1.fp.plot
tmp <- qr.2#
#
# Frequency and predictability#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.05, to=+0.02, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.05, +0.03) )#
rq2.fp.plot <- p + geom_hline()
qr.3 <- by(d, d$id, function(x) rq(ldur ~ f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao #
                                        + I(f^2)   + I(f^3)   + I(p^2)  + I(p^3)  + I(l^2) #
                                        + I(f1^2)             + I(p1^2)           + I(l1^2) #
                                        + I(f2^2)  + I(f2^3)                      + I(l2^2) #
                                        + I(lsr^2) + I(lsr^3) #
                                                   + I(f:l), #
									   tau=seq(10,90,10)/100,data=x))
qr.3 <- by(d, d$id, function(x) rq(ldur ~ f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao #
                                        + I(f^2)   + I(f^3)   + I(p^2)  + I(p^3)  + I(l^2) #
                                        + I(f1^2)             + I(p1^2)           + I(l1^2) #
                                        + I(f2^2)  + I(f2^3)                      + I(l2^2) #
                                        + I(lsr^2) + I(lsr^3) #
                                                   + I(f*l), #
									   tau=seq(10,90,10)/100,data=x))
tmp <- qr.3#
#
# Frequency and predictability#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.05, to=+0.02, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.05, +0.03) )#
rq3.fp.plot <- p + geom_hline()
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.05, to=+0.02, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
rq3.fp.plot <- p + geom_hline()
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
rq3.fp.plot <- p + geom_hline()
tmp <- qr.2#
#
# Frequency and predictability#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
rq2.fp.plot <- p + geom_hline()
tmp <- qr.3
set <- c(2, 15, 16, 5, 17, 18)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f", "I(f^2", "I(f^3)", "p", "I(p^2", "I(p^3)"))#
levels(df$par) <- c("N Freq", "N Freq^2", "N Freq^3", "N Pred", "N Pred^2", "N Pred^3")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
rq3.fp.poly.plot <- p + geom_hline()
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f", "I(f^2)", "I(f^3)", "p", "I(p^2)", "I(p^3)"))#
levels(df$par) <- c("N Freq", "N Freq^2", "N Freq^3", "N Pred", "N Pred^2", "N Pred^3")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
rq3.fp.poly.plot <- p + geom_hline()
rq3.fp.plot
rq2.fp.plot
rq3.fp.poly.plot
tmp <- qr.1  # use qr.1, qr.2, qr.3 #
#
set <- c(8:14)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)", "ao"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Peak", "IOVP Curvature", "Saccade length")
set
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par, scales="free_x")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient") #, breaks=seq(from=-0.05, to=+0.02, by=0.01) )#
#p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
qr1.loa.plot <- p + geom_hline()
qr1.loa.plot
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(par ~ ., scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient") #, breaks=seq(from=-0.05, to=+0.02, by=0.01) )#
#p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(par ~ ., scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.05, to=+0.05, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.5, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
tmp <- qr.2  # use qr.1, qr.2, qr.3 #
#
set <- c(8:14)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)", "ao"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Peak", "IOVP Curvature", "Saccade length")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.5, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr2.loa.plot <- p + geom_hline())
tmp <- qr.1  # use qr.1, qr.2, qr.3 #
#
set <- c(8:14)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)", "ao"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Peak", "IOVP Curvature", "Saccade length")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
(qr1.loa.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.5, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
tmp <- qr.3  # use qr.1, qr.2, qr.3 #
#
set <- c(8:14)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)", "ao"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Peak", "IOVP Curvature", "Saccade length")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.5, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr3.loa.plot <- p + geom_hline())
tmp <- qr.3  # use qr.1, qr.2, qr.3 #
#
set <- c(8:14)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)", "ao"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Peak", "IOVP Curvature", "Saccade length")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.5, +0.5) )#
(qr3.loa.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr3.loa.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr3.loa.plot <- p + geom_hline())
tmp <- qr.1  # use qr.1, qr.2, qr.3 #
#
set <- c(8:14)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)", "ao"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Peak", "IOVP Curvature", "Saccade length")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(2,1))) #
#
print(qr1.loa.plot, vp=vplayout(1,1)) #
print(qr3.loa.plot, vp=vplayout(2,1))
# Length, launch site, relative fixation position (curvature)#
tmp <- qr.1  # use qr.1, qr.2, qr.3 #
#
set <- c(8:11, 13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "I(o^2)"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
tmp <- qr.3  # use qr.1, qr.2, qr.3 #
#
set <- c(8:11, 13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "I(o^2)"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_grid(. ~ par , scales="free_x", space="free")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr3.loa.plot <- p + geom_hline())
# Combine the plots#
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(2,1))) #
#
print(qr1.loa.plot, vp=vplayout(1,1)) #
print(qr3.loa.plot, vp=vplayout(2,1))
# Reference 40-parameter model (full 3rd degree polynomial + significant interactions); won't converge obviously#
qr.0 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao #
                                        + I(f^2)   + I(f^3)   + I(p^2)  + I(p^3)  + I(l^2) + I(l^3)#
                                        + I(f1^2)  + I(f1^3)  + I(p1^2) + I(p1^3) + I(l1^2) + I(l1^3)#
                                        + I(f2^2)  + I(f2^3)  + I(p2^2) + I(p2^3) + I(l2^2) + I(l2^3)#
                                        + I(lsr^2) + I(lsr^3) + I(ao^2) + I(ao^3)#
                                        + I(f*f1)  + I(f*l)   + I(f2*l) + I(p2*l), #
                                       tau=seq(10,90,10)/100,data=x))
tmp <- qr.0  # use qr.1, qr.2 for simpler models#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
qr0.fp.plot <- p + geom_hline()
qr0.fp.plot
tmp <- qr.1  # use qr.1, qr.2 for simpler models#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr1.fp.plot <- p + geom_hline())
tmp <- qr.2  # use qr.1, qr.2 for simpler models#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr2.fp.plot <- p + geom_hline())
tmp <- qr.3  # use qr.1, qr.2 for simpler models#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr3.fp.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) ) + title("28-Parameter Model")#
qr3.fp.poly3.plot <- p + geom_hline()
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) ) + geom_title("28-Parameter Model")#
qr3.fp.poly3.plot <- p + geom_hline()
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) ) + main("28-Parameter Model")#
qr3.fp.poly3.plot <- p + geom_hline()
p <- ggplot(df, aes(x=Quantile, y=value, main("28-Parameter Model") ) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) ) #
qr3.fp.poly3.plot <- p + geom_hline()
qr3.fp.poly3.plot
p <- qplot(x=Quantile, y=value, df, main="28-Parameter Model", facet_wrap( ~ par, ncol=3))
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- qplot(x=Quantile, y=value, data=df, main="40-Parameter Model") + facet_wrap( ~ par, ncol=3)
p <- ggplot(df, aes(x=Quantile, y=value, main=title) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value), main=title ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title=title, aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)
p
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7, 12:14)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2", "o", "I(o^2)", "ao"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred", "IOVP Peak", "IOVP Curvature", "Saccade")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
set <- c("f1", "f", "f2", "p1", "p", "p2")#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]
# 14-parameter main effects only model (Kliegl et al., 2006; using (1/launch site) and 1/ao)#
qr.1 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao, #
									    tau=seq(10,90,10)/100,data=x))
# 19-parameter main effects and interactions model (Kliegl et al., 2006; excl. f:f2; using (1/launch site) and 1/ao)#
qr.2 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao#
                                        + I(f*f1)  + I(f*l)   + I(f2*l) + I(p2*l), #
                                       tau=seq(10,90,10)/100,data=x))
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr1.fp.plot <- p + geom_hline())
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14 vs. 40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr1.fp.plot <- p + geom_hline())#
(qr0.fp.plot <- p %+% qr.1
)
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14 vs. 40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
tmp <- qr.1  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr1.fp.plot <- p + geom_hline())
set <- c(8:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.loa.plot <- p + geom_hline())
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(8:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.loa.plot <- p + geom_hline())
tmp <- qr.1  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(8:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr1.loa.plot <- p + geom_hline())
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(2,2))) #
#
print(qr0.fp.plot, vp=vplayout(1,1)) #
print(qr1.fp.plot, vp=vplayout(1,2)) #
print(qr0.loa.plot, vp=vplayout(2,1)) #
print(qr1.loa.plot, vp=vplayout(2,2))
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$value10[1:nrow(df)/2] <- df$value[1:nrow(df)/2]*10  # Multiply f and p coefficients by 10 to get on common scale#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.fp.plot <- p + geom_hline())#
 #
# Combine plots#
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr0.fp.plot, vp=vplayout(1,1)) #
print(qr1.fp.plot, vp=vplayout(1,2))
tmp <- qr.1  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$value10[1:nrow(df)/2] <- df$value[1:nrow(df)/2]*10  # Multiply f and p coefficients by 10 to get on common scale#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr1.fp.plot <- p + geom_hline())#
 #
# Combine plots#
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr0.fp.plot, vp=vplayout(1,1))
?any
?set
?union
1:10 %in% c(1,3,5,9)
ix <- c("f", "f1", "f2", "p", "p1", "p2") %in% df$par
ix
df$par %in% c("f", "f1", "f2", "p", "p1", "p2")
ix <- df$par %in% c("f", "f1", "f2", "p", "p1", "p2")
length(ix)
nrow(df)
head(df)
ix <- df$par %in% c("N Freq", "N-1 Frea", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
df$value10[ix] <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.plot <- p + geom_hline())
df$value <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.plot <- p + geom_hline())
tmp <- qr.1  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
df$value <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr1.plot <- p + geom_hline())#
 #
# Combine plots#
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr0.fp.plot, vp=vplayout(1,1)) #
print(qr1.fp.plot, vp=vplayout(1,2))
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr0.plot, vp=vplayout(1,1)) #
print(qr1.plot, vp=vplayout(1,2))
tmp <- qr.3  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
df$value <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="28-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr3.plot <- p + geom_hline())
tmp <- qr.2  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
df$value <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="28-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr2.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="19-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr2.plot <- p + geom_hline()
tmp <- qr.3  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
df$value <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="28-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr3.plot <- p + geom_hline()
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr3.plot, vp=vplayout(1,1)) #
print(qr2.plot, vp=vplayout(1,2))
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
df$value.10 <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value.10) ) + opts(title="28-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr0.plot <- p + geom_hline()
rm(list=ls())#
#
# http://had.co.nz/ggplot2/stat_summary.html#
stat_sum_single <- function(fun, geom="point", ...) { #
   stat_summary(fun=fun, colour="blue", geom=geom, size = 2, ...) #
}#
#
stat_sum_df <- function(fun, geom="smooth", ...) {#
  stat_summary(fun=fun, colour="red", geom=geom, method=lm, width=0.2, size=1,...)#
}#
#
#
load("PSC275_fpf.sfc_R.rda")#
#
d$ldur = log(d$dur)#
#
# Reference 40-parameter model (full 3rd degree polynomial + significant interactions)-- surprisingly, it converges!#
qr.0 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao #
                                        + I(f^2)   + I(f^3)   + I(p^2)  + I(p^3)  + I(l^2) + I(l^3)#
                                        + I(f1^2)  + I(f1^3)  + I(p1^2) + I(p1^3) + I(l1^2) + I(l1^3)#
                                        + I(f2^2)  + I(f2^3)  + I(p2^2) + I(p2^3) + I(l2^2) + I(l2^3)#
                                        + I(lsr^2) + I(lsr^3) + I(ao^2) + I(ao^3)#
                                        + I(f*f1)  + I(f*l)   + I(f2*l) + I(p2*l), #
                                       tau=seq(10,90,10)/100,data=x))#
#
# 14-parameter main effects only model (Kliegl et al., 2006; using (1/launch site) and 1/ao)#
qr.1 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao, #
									    tau=seq(10,90,10)/100,data=x))#
									#
# 19-parameter main effects and interactions model (Kliegl et al., 2006; excl. f:f2; using (1/launch site) and 1/ao)#
qr.2 <- by(d, d$id, function(x) rq(ldur ~ 1 + f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao#
                                        + I(f*f1)  + I(f*l)   + I(f2*l) + I(p2*l), #
                                       tau=seq(10,90,10)/100,data=x))#
#
# 14+14-parameter + allowing for obvious polynomial trends (Kliegl, 2007; Kliegl, 2009) plus f:l interaction#
qr.3 <- by(d, d$id, function(x) rq(ldur ~ f + f1 + f2 + p + p1 + p2 + l + l1 + l2 + lsr + o + I(o^2) + ao #
                                        + I(f^2)   + I(f^3)   + I(p^2)  + I(p^3)  + I(l^2) #
                                        + I(f1^2)             + I(p1^2)           + I(l1^2) #
                                        + I(f2^2)  + I(f2^3)                      + I(l2^2) #
                                        + I(lsr^2) + I(lsr^3) #
                                                   + I(f*l), #
									   tau=seq(10,90,10)/100,data=x))#
#
# Frequency and predictability. #
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
df$value.10 <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale#
#
p <- ggplot(df, aes(x=Quantile, y=value.10) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr0.plot <- p + geom_hline()
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(8:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.loa.plot <- p + geom_hline())
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr0.fp.plot, vp=vplayout(1,1)) #
print(qr0.loa.plot, vp=vplayout(1,2))
tmp <- qr.1  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr1.fp.plot <- p + geom_hline())#
#
# Length, launch site, relative fixation position (curvature)#
tmp <- qr.1  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(8:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr1.lo.plot <- p + geom_hline())#
#
# Combine plots#
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr1.fp.plot, vp=vplayout(1,1)) #
print(qr1.lo.plot, vp=vplayout(1,2))
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")
tapply(df$value, df$par, mean)
cast( . ~ par, mean, df)
cast( . ~ par, mean, data=df)
options(digits=3)
options(digits=2)
options(digits=1)
cast(par ~ ., mean, data=df)
table(ix)
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)
value.10 <- df$value[ix]*10         # Multiply f and p coefficients by 10 to get on common scale
table(ix,par)
table(ix,df$par)
table(df$par,ix)
length(value.10)
value.10 <- ifelse(ix, df$value[ix]*10, df$value)
tapply(value.10, df$par, mean)
tapply(value.10, list(df$par,ix), mean)
tapply(value.10, list(df$par,ix), sd)
df$value.10 <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)
str(dv)
str(df)
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)
p <- ggplot(df, aes(x=Quantile, y=value.10) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr0.plot <- p + geom_hline()
qr0.plot
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ Quantile, function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ Quantile, function(x) c(M=mean(x)), data=df)
tmp <- qr.1 # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
p <- ggplot(df, aes(x=Quantile, y=value.10) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr1.plot <- p + geom_hline()
qr1.plot
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr0.plot, vp=vplayout(1,1)) #
print(qr1.plot, vp=vplayout(1,2))
tmp <- qr.0 # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr0.plot <- p + geom_hline()
tmp <- qr.0 # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr0.plot <- p + geom_hline()
tmp <- qr.1 # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ ., function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
qr1.plot <- p + geom_hline()
vplayout <- function(x, y) #
viewport(layout.pos.row = x, layout.pos.col = y) #
#
grid.newpage() #
pushViewport(viewport(layout = grid.layout(1,2))) #
#
print(qr0.plot, vp=vplayout(1,1)) #
print(qr1.plot, vp=vplayout(1,2))
tmp <- qr.0 # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ decile, function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ decile, function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ quantile, function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ Quantile, function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ Quantile, function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ Quantile, function(x) c(M=mean(x)), data=df)
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
cast(par ~ Quantile, function(x) c(M=mean(x)), data=df)
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
cast(par ~ Quantile, function(x) c(M=mean(x)), data=df)
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
cast(par ~ Quantile, function(x) c(M=mean(x)), data=df)#
#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]
tmp <- qr.0#
#
set <- c(2, 15, 16, 5, 17, 18, 3, 21, 22, 6, 23, 24, 6, 28, 29, 7, 30, 31)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f",  "I(f^2)",  "I(f^3)",  "p",  "I(p^2)",  "I(p^3)", #
                                  "f1", "I(f1^2)", "I(f1^3)", "p1", "I(p1^2)", "I(p1^3)",#
                                  "f2", "I(f2^2)", "I(f2^3)", "p1", "I(p2^2)", "I(p2^3)") )#
#
levels(df$par) <- c("N Freq",   "N Freq^2",   "N Freq^3",   "N Pred",   "N Pred^2", "N Pred^3", #
                    "N-1 Freq", "N-1 Freq^2", "N-1 Freq^3", "N-1 Pred", "N-1 Pred^2", "N-1 Pred^3",#
                    "N+1 Freq", "N+1 Freq^2", "N+1 Freq^3", "N+1 Pred", "N+1 Pred^2", "N+1 Pred^3")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
qr0.fp.poly3.plot <- p + geom_hline()
tmp <- qr.0#
#
set <- c(2, 15, 16, 5, 17, 18, 3, 21, 22, 6, 23, 24, 4, 27, 28, 7, 29, 30)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f",  "I(f^2)",  "I(f^3)",  "p",  "I(p^2)",  "I(p^3)", #
                                  "f1", "I(f1^2)", "I(f1^3)", "p1", "I(p1^2)", "I(p1^3)",#
                                  "f2", "I(f2^2)", "I(f2^3)", "p1", "I(p2^2)", "I(p2^3)") )#
#
levels(df$par) <- c("N Freq",   "N Freq^2",   "N Freq^3",   "N Pred",   "N Pred^2", "N Pred^3", #
                    "N-1 Freq", "N-1 Freq^2", "N-1 Freq^3", "N-1 Pred", "N-1 Pred^2", "N-1 Pred^3",#
                    "N+1 Freq", "N+1 Freq^2", "N+1 Freq^3", "N+1 Pred", "N+1 Pred^2", "N+1 Pred^3")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
qr0.fp.poly3.plot <- p + geom_hline()
tmp <- qr.0#
#
set <- c(2, 15, 16, 5, 17, 18, 3, 21, 22, 6, 23, 24, 4, 27, 28, 7, 29, 30)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f",  "I(f^2)",  "I(f^3)",  "p",  "I(p^2)",  "I(p^3)", #
                                  "f1", "I(f1^2)", "I(f1^3)", "p1", "I(p1^2)", "I(p1^3)",#
                                  "f2", "I(f2^2)", "I(f2^3)", "p2", "I(p2^2)", "I(p2^3)") )#
#
levels(df$par) <- c("N Freq",   "N Freq^2",   "N Freq^3",   "N Pred",   "N Pred^2", "N Pred^3", #
                    "N-1 Freq", "N-1 Freq^2", "N-1 Freq^3", "N-1 Pred", "N-1 Pred^2", "N-1 Pred^3",#
                    "N+1 Freq", "N+1 Freq^2", "N+1 Freq^3", "N+1 Pred", "N+1 Pred^2", "N+1 Pred^3")#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
qr0.fp.poly3.plot <- p + geom_hline()
names(qr.0)
levels(qr.0$par)
str(tmp)
names(tmp)
head(tmp)
tmp <- qr.0#
#
set <- c(2:7, 15, 21, 27, 17, 23, 29, 16, 22, 28,  18, 24, 30)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f",       "f1",      "f2",      "p",     "p1",      "p2",#
                                  "I(f^2)", "I(f1^2)", "I(f2^2)", "I(p^2)", "I(p1^2)", "I(p2^2)",#
                                  "I(f^3)", "I(f1^3)", "I(f2^3)", "I(p^3)", "I(p1^3)", "I(p2^3)")#
#
levels(df$par) <- c("N Freq",   "N-1 Freq",   "N+1 Freq",   "N Pred",   "N-1 Pred",    "N+1 Pred", #
                    "N Freq^2", "N-1 Freq^2", "N+1 Freq^2", "N Pred^2", "N-1 Pred^2",  "N+1 Pred^2",#
                    "N Freq^3", "N-1 Freq^3", "N+1 Freq^3", "N Pred^3", "N-1 Pred^3",  "N+1 Pred^3" )#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
qr0.fp.poly3.plot <- p + geom_hline()
tmp <- qr.0#
#
set <- c(2:7, 15, 21, 27, 17, 23, 29, 16, 22, 28,  18, 24, 30)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f",       "f1",      "f2",      "p",     "p1",      "p2",#
                                  "I(f^2)", "I(f1^2)", "I(f2^2)", "I(p^2)", "I(p1^2)", "I(p2^2)",#
                                  "I(f^3)", "I(f1^3)", "I(f2^3)", "I(p^3)", "I(p1^3)", "I(p2^3)") )#
#
levels(df$par) <- c("N Freq",   "N-1 Freq",   "N+1 Freq",   "N Pred",   "N-1 Pred",    "N+1 Pred", #
                    "N Freq^2", "N-1 Freq^2", "N+1 Freq^2", "N Pred^2", "N-1 Pred^2",  "N+1 Pred^2",#
                    "N Freq^3", "N-1 Freq^3", "N+1 Freq^3", "N Pred^3", "N-1 Pred^3",  "N+1 Pred^3" )
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
qr0.fp.poly3.plot <- p + geom_hline()
qr0.fp.poly3.plot
set <- c(8:11, 14, 12, 19, 25, 31, 33, 35, 13, 20, 26, 32, 34, 36)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l",      "l1",      "l2",      "lrs",      "ao",      "o",  #
                                  "I(l^2)", "I(l1^2)", "I(l2^2)", "I(lrs^2)", "I(ao^2)", "I(o^2)",#
                                  "I(l^3)", "I(l1^3)", "I(l2^3)", "I(lrs^3)", "I(ao^3)") )#
#
levels(df$par) <- c("N Length",   "N-1 Length",   "N+1 Length",   "Launch site",   "Saccade",  "IOVP Peak", #
                    "N Length^2", "N-1 Length^2", "N+1 Length^2", "Launch site^2", "Saccade^2",  "IOVP Curve",#
                    "N Length^3", "N-1 Length^3", "N+1 Length^3", "Launch site^3", "Saccade^3" )#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.3, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.3) )#
qr0.lo.poly.plot <- p + geom_hline()
par
levels(df$par)
set <- c(8:11, 14, 12)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l",      "l1",      "l2",      "lrs",      "ao",      "o") )#
#
levels(df$par) <- c("N Length",   "N-1 Length",   "N+1 Length",   "Launch site",   "Saccade",  "IOVP Peak" )#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.3, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.3) )#
p + geom_hline()
tmp <- qr.0#
#
set <- c(8:11, 14, 12)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l",      "l1",      "l2",      "lrs",      "ao",      "o") )#
#
levels(df$par) <- c("N Length",   "N-1 Length",   "N+1 Length",   "Launch site",   "Saccade",  "IOVP Peak" )#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.3, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.3) )#
p + geom_hline()
tmp <- qr.0#
#
set <- c(8:11, 14, 12)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l",      "l1",      "l2",      "lsr",      "ao",      "o") )#
#
levels(df$par) <- c("N Length",   "N-1 Length",   "N+1 Length",   "Launch site",   "Saccade",  "IOVP Peak" )#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.3, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.3) )#
p + geom_hline()
tmp <- qr.0#
#
set <- c(8:11, 14, 12, 19, 25, 31, 33, 35, 13, 20, 26, 32, 34, 36)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("l",      "l1",      "l2",      "lsr",      "ao",      "o",  #
                                  "I(l^2)", "I(l1^2)", "I(l2^2)", "I(lsr^2)", "I(ao^2)", "I(o^2)",#
                                  "I(l^3)", "I(l1^3)", "I(l2^3)", "I(lsr^3)", "I(ao^3)") )#
#
levels(df$par) <- c("N Length",   "N-1 Length",   "N+1 Length",   "Launch site",   "Saccade",  "IOVP Peak", #
                    "N Length^2", "N-1 Length^2", "N+1 Length^2", "Launch site^2", "Saccade^2",  "IOVP Curve",#
                    "N Length^3", "N-1 Length^3", "N+1 Length^3", "Launch site^3", "Saccade^3" )#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + facet_wrap( ~ par, ncol=6)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.3, by=0.1) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.3) )#
qr0.lo.poly.plot <- p + geom_hline()
qr0.lo.poly.plot
tmp <- qr.0  # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:7)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
df$Quantile <- as.numeric(df$Quantile)#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred")#
#
cast(par ~ Quantile, function(x) c(M=mean(x)), data=df)#
#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="40-Parameter Model", aspect=1) + facet_wrap( ~ par, ncol=3)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.07, to=+0.03, by=0.01) )#
p <- p + coord_cartesian(ylim=c(-0.07, +0.03) )#
(qr0.fp.plot <- p + geom_hline())
tmp <- qr.0 # use qr.1, qr.2, qr.3 for comparisons with simpler models #
# Note: Modify plot title in opts() to match qr model!#
#
set <- c(2:13)#
id  <- names(tmp)#
par <- names(tmp[[1]]$coef[ ,1])[set]#
#
df1 <- data.frame(id = rep(id, each=length(set)), par = rep(par, length(id)) )#
df2 <- data.frame()#
for (i in id) df2 <- rbind(df2, tmp[[i]]$coef[set, ])#
names(df2) <- 1:9#
 #
df <- melt(cbind(df1, df2), id=c("id", "par"), variable_name="Quantile")#
#
df$Quantile <- as.numeric(df$Quantile)#
#
df$par <- factor(df$par, levels=c("f1", "f", "f2", "p1", "p", "p2",#
                                  "l1", "l", "l2", "lsr", "o", "I(o^2)"))#
levels(df$par) <- c("N-1 Freq", "N Freq", "N+1 Freq", "N-1 Pred", "N Pred", "N+1 Pred",#
                    "N-1 Length", "N Length", "N+1 Length", "Launch site", "IOVP-Peak", "IOVP Curvature")#
#
ix <- df$par %in% c("N Freq", "N-1 Freq", "N+1 Freq", "N Pred", "N-1 Pred", "N+1 Pred")#
cast(par ~ Quantile, function(x) c(M=mean(x)), data=df)#
#
#df$value <- ifelse(ix, df$value[ix]*10, df$value)   # compensate for sd * 10 difference#
cast(par ~ Quantile, function(x) c(M=mean(x), SD=sd(x), N=length(x)), data=df)#
#
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2)#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient", breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2, scales="free_y")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient") # breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.plot <- p + geom_hline())
p <- ggplot(df, aes(x=Quantile, y=value) ) + opts(title="14-Parameter Model", aspect=1) + facet_wrap( ~ par, nrow=2, scales="free_y")#
p <- p + stat_sum_df("mean_cl_normal") + stat_sum_single("mean", size=2)#
p <- p + scale_x_continuous("Quantile of log(single-fixation duration)", breaks=1:9 )#
p <- p + scale_y_continuous("QR Coefficient") # breaks=seq(from=-0.7, to=+0.5, by=0.2) )#
#p <- p + coord_cartesian(ylim=c(-0.7, +0.5) )#
(qr0.plot <- p + geom_hline())
df$FP <- ix#
levels(d$FP) <- c("Freq/Pred", "Length/Launch/IOVP")
setwd('/Users/kliegl/documents/projects/EyeMoves/KNE_R8')
df$FP <- ifelse(ix, 1, 0)#
levels(d$FP) <- c("Freq/Pred", "Length/Launch/IOVP")
rm(list=ls())#
ps <- read.csv("scopus.all.psycsci.csv", header=TRUE)#
ps <- ps[ps$Year > 1995, -c(3, 5)]
ps <- read.csv("scopus.all.psycsci.csv", header=TRUE)
dir()
library(reshape)#
library(ggplot2)#
#
Year <- 2007:1996#
all <-      c(329, 421, 393, 372, 276, 254, 262, 219, 227, 204, 208, 165)#
D.abs <-    c(25, 30, 25, 22, 25, 13, 16, 10, 7, 5, 5, 11)#
UK.abs <-   c(39, 53, 48, 48, 41, 27, 21, 26, 22, 16, 16, 12)#
NL.abs <-   c(16, 20, 17, 25, 14, 11, 7, 4, 2, 6, 3, 3)#
CA.abs <-   c(31, 42, 46, 40, 27, 27, 33, 31, 23, 19, 18, 18)#
ANZ.abs <-  c(21, 18, 26, 19, 13, 15, 16, 14, 9, 10, 6, 4)  # Australia, New Zealand#
F.abs <-    c(8, 14, 14, 14, 7, 6, 11, 3, 5, 7, 4, 2) #
AS.abs <-   c(5, 5, 2, 4, 1, 0, 1, 0, 0, 1, 0, 0)           # Austria, Switzerland#
SC.abs <-   c(6, 4, 6, 3, 7, 2, 6, 4, 1, 3, 3, 1)           # SE, FI, NO, DK, IL#
SE.abs <-   c(12, 20, 11, 5, 10, 9, 4, 6, 4, 4, 0, 4)       # Italy, Spain, P, Greece#
ASIA.abs <- c(7, 15, 7, 4, 5, 5, 6, 3, 1, 4, 2, 0)          # China, Taiwan, Japan, India#
IS.abs <-   c(5, 13, 20, 9, 14, 11, 4, 5, 5, 4, 7, 2)       # Israel#
#
d <- data.frame(Year, all, D.abs, UK.abs, NL.abs, CA.abs, ANZ.abs, F.abs, AS.abs,#
	    SC.abs, SE.abs, ASIA.abs,IS.abs,#
	    D=D.abs/all, UK=UK.abs/all, NL=NL.abs/all, CA=CA.abs/all, ANZ=ANZ.abs/all,#
	    F=F.abs/all, AS=AS.abs/all, SC=SC.abs/all, SE=SE.abs/all, ASIA=ASIA.abs/all, IS=IS.abs/all)#
#
d$US <- 1-rowSums(d[,14:24])	#
ix <- which(d$Year >= 1997 & d$Year <=2006)  # 10-year slot#
d <- d[ix, ]#
d.rs <- melt(d, id="Year", measure=c("D", "UK", "CA", "ANZ", "IS"), var="Country")#
qplot(Year, value, data=d.rs, colour=Country, geom=c("point", "smooth"), method=lm, formula=y~poly(x,2),#
       ylab="Proportion of Publications")
d.rs <- melt(d, id="Year", measure=c("D.abs", "UK.abs",  "CA.abs", "ANZ.abs", "F.abs", "NL.abs"), var="Country")#
qplot(Year, value, data=d.rs, colour=Country, geom=c("point", "smooth"), method=lm, formula=y~poly(x,2))
setwd('/Users/kliegl/documents/manuscripts/APS')
d.rs <- melt(d, id="Year", measure=c("D", "UK", "CA", "ANZ", "US", "IS"), var="Country")#
qplot(Year, value, data=d.rs, colour=Country, geom=c("point", "smooth"), method=lm, formula=y~poly(x,2),#
       ylab="Proportion of Publications")
