m$yloca.norm <- round(m$y.loca/max(data$y.loca),digits=1)
m$xloca.norm <- round(data$x.loca/max(data$x.loca),digits=1)#
m$yloca.norm <- round(m$y.loca/max(data$y.loca),digits=1)
m <- data#
m$xloca.norm <- round(data$x.loca/max(data$x.loca),digits=1)#
m$yloca.norm <- round(m$y.loca/max(data$y.loca),digits=1)
m$xy <- interaction(m$xloca.norm,m$yloca.norm,drop=TRUE,sep=":")#
#m$freq <- rep(1,length(data$z.loca))#
ids <- c("xloca.norm","yloca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
ed(m)#
m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
m <- data#
m$xloca.norm <- round(data$x.loca/max(data$x.loca),digits=1)#
m$yloca.norm <- round(m$y.loca/max(data$y.loca),digits=1)#
m$xy <- interaction(m$xloca.norm,m$yloca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("xloca.norm","yloca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
ed(m)#
m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
m <- data#
m$xloca.norm <- round(data$xloca/max(data$xloca),digits=1)#
m$yloca.norm <- round(m$yloca/max(data$yloca),digits=1)#
m$xy <- interaction(m$xloca.norm,m$yloca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("xloca.norm","yloca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
#ed(m)#
m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
data$xloca
data
data$xloca
data$yloca
data$zloca
colnames(data)
rm(data)
ls()
rm(list=ls())
ls()
data <- read.delim("/Users/ackman/Desktop/120703_01_STATS-Centroids.txt",sep=" ")
nrow(data)
data$xloca
data
colnames(data)
data$x.loca
length(data$x.loca)
length(data$y.loca)
length(data$z.loca)
data$x.loca/max(data$x.loca)
length(data$x.loca)
is.numeric(data$x.loca)
max(data$x.loca)
rm(data)
data <- read.delim(pipe("pbpaste"))
colnames(data)
nrow(data)
max(data$x.loca)
data$x.loca
length(data$x.loca)
mean(data$x.loca)
is.vector(data$x.loca)
is.numeric(data$x.loca)
quit()
require(reshape)
data <- read.delim(pipe("pbpaste"))
data
with(data,mean(x.loca))
with(data,is.numeric(x.loca))
help(is.numeric)
with(data,is.double(x.loca))
data <- read.delim(pipe("pbpaste"))
with(data,max(x.loca))
with(data,is.double(x.loca))
with(data,length(x.loca))
with(data,mean(x.loca))
with(data,is.double(x.loca))
with(data,lenght(x.loca))
with(data,length(x.loca))
data
rm(data)
data <- read.delim(pipe("pbpaste"),sep=" ")
with(data,length(x.loca))
with(data,max(x.loca))
with(data,mean(x.loca))
require(reshape)#
#
#-----------R code for normalized frequency contour plot-------------#
quartz();#
m <- data#
m$x.loca.norm <- round(data$x.loca/max(data$x.loca),digits=1)#
m$y.loca.norm <- round(m$y.loca/max(data$y.loca),digits=1)#
m$xy <- interaction(m$x.loca.norm,m$y.loca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("x.loca.norm","y.loca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
#ed(m)#
m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
m <- data#
m$x.loca.norm <- round(data$x.loca/max(data$x.loca),digits=2)#
m$y.loca.norm <- round(m$y.loca/max(data$y.loca),digits=2)#
m$xy <- interaction(m$x.loca.norm,m$y.loca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("x.loca.norm","y.loca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
#ed(m)#
m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
max(m$freq)
max(m$x.loca)
max(m$x.loca.norm)
mean(m$x.loca)
mean(m$x.loca.norm)
colnames(d4)
colnames(m)
m$x.loca.norm
m <- data#
m$x.loca.norm <- round(data$x.loca/max(data$x.loca),digits=3)#
m$y.loca.norm <- round(m$y.loca/max(data$y.loca),digits=3)#
m$xy <- interaction(m$x.loca.norm,m$y.loca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("x.loca.norm","y.loca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
#ed(m)#
#m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.001)#
y1 <- seq(0,1.0,by=0.001)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
matchidx
x1
m$x.loca.norm
x1
m$x.loca.norm
m <- data#
m$x.loca.norm <- round(data$x.loca/max(data$x.loca),digits=1)#
m$y.loca.norm <- round(m$y.loca/max(data$y.loca),digits=1)#
m$xy <- interaction(m$x.loca.norm,m$y.loca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("x.loca.norm","y.loca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
#ed(m)#
m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
m <- data#
m$x.loca.norm <- round(data$x.loca/max(data$x.loca),digits=1)#
m$y.loca.norm <- round(m$y.loca/max(data$y.loca),digits=1)#
m$xy <- interaction(m$x.loca.norm,m$y.loca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("x.loca.norm","y.loca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
#ed(m)#
#m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"),asp=1)
m <- data#
m$x.loca.norm <- round(data$x.loca/max(data$x.loca),digits=1)#
m$y.loca.norm <- round(m$y.loca/max(data$y.loca),digits=1)#
m$xy <- interaction(m$x.loca.norm,m$y.loca.norm,drop=TRUE,sep=":")#
m$freq <- rep(1,length(data$z.loca))#
ids <- c("x.loca.norm","y.loca.norm","xy")#
meas <- c("freq")#
m <- as.data.frame(m)#
d3 <- with(m,melt(m,id.var=ids,measure.var=meas))#
d4 <- cast(d3,fun.aggregate=sum)#
m <- d4#
#ed(m)#
#m$freq <- m$freq/max(m$freq)#
#
#make continous frame for a matrix for inputs to contour#
x1 <- seq(0,1.0,by=0.1)#
y1 <- seq(0,1.0,by=0.1)#
df <- expand.grid(x=x1,y=y1)#
df$xy <- interaction(df$x,df$y,drop=TRUE,sep=":")#
matchidx <- match(df$xy,m$xy)#
tmp <- m$freq[matchidx]#
tmp[which(is.na(tmp))] <- 0#
df$z <- tmp#
tmp <- data.frame(x=df$x,y=df$y,z=df$z)#
mat <- matrix(tmp$z,nrow=sqrt(nrow(tmp)))#
#filled.contour(mat,color=terrain.colors)#
filled.contour(mat,color=colorRampPalette(c("blue", "white", "red"), space = "rgb"))
max(data$x.loca)
max(data$y.loca)
help(matrix)
tmp <– matrix(data = rep(0,512*696), nrow = 512, ncol = 696)
tmp <– matrix(nrow = 512, ncol = 696)
rm(tmp)
tmp <– matrix(nrow = 512, ncol = 696)
tmp<–matrix(nrow = 512, ncol = 696)
help(matrix)
tmp<–matrix(data = 0,nrow = 512, ncol = 696)
matrix(data = 0,nrow = 512, ncol = 696)
mat <- matrix(data = 0,nrow = 512, ncol = 696)
tmp
tmp <- matrix(data = 0,nrow = 512, ncol = 696)
rm(tmp)
dim(mat)
image(tmp)
image(mat)
rm(mat)
rm(tmp)
rm(mat)
m <- data#
m$x.loca.norm <- round(data$x.loca)#
m$y.loca.norm <- round(data$y.loca)#
tmp <– matrix(data = 0, nrow = 512, ncol = 696)#
#
for(i in length(m$x.loca.norm)) {#
tmp[m$y.loca.norm(i),m$x.loca.norm(i)] <- tmp[m$y.loca.norm(i),m$x.loca.norm(i)] + 1#
}
tmp <– matrix(data = 0, nrow = 512, ncol = 696)
matrix(data = 0, nrow = 512, ncol = 696)
tmp = matrix(data = 0, nrow = 512, ncol = 696)
for(i in length(m$x.loca.norm)) {#
tmp[m$y.loca.norm(i),m$x.loca.norm(i)] <- tmp[m$y.loca.norm(i),m$x.loca.norm(i)] + 1#
}
for (i in 1:length(m$x.loca.norm)) {#
tmp[m$y.loca.norm(i),m$x.loca.norm(i)] <- tmp[m$y.loca.norm(i),m$x.loca.norm(i)] + 1#
}
1:length(m$x.loca.norm))
1:length(m$x.loca.norm)
for (i in 1:length(m$x.loca.norm)) {#
#tmp[m$y.loca.norm(i),m$x.loca.norm(i)] <- tmp[m$y.loca.norm(i),m$x.loca.norm(i)] + 1#
tmp[m$y.loca.norm(i),m$x.loca.norm(i)] <- 1#
}
for (i in 1:length(m$x.loca.norm)) {#
tmp[m$y.loca.norm[i],m$x.loca.norm[i]] <- tmp[m$y.loca.norm[i],m$x.loca.norm[i]] + 1#
}
image(tmp)
help(image)
heatmap(tmp)
quit()
NA  <- 0.5
NA <- 0.5
numapp <- 0.5
ls <- 200
lambda <- 925
n <- 1.33
((2*pi*(numapp^2))/(lambda*n*ls))*(z^2)*exp(-2*z/ls)
z <- 1760
((2*pi*(numapp^2))/(lambda*n*ls))*(z^2)*exp(-2*z/ls)
quit()
help(igraph)
require(igraph)
library(help="igraph")
help(walktrap.community)
g <- erdos.renyi.game(10, 5/10) %du% erdos.renyi.game(9, 5/9)
g <- add.edges(g, c(0, 11))
g <- subgraph(g, subcomponent(g, 0))
spinglass.community(g, spins=2)
spinglass.community(g, vertex=0)
g
(19*18)/2
40/171
g <- graph.full(5) %du% graph.full(5) %du% graph.full(5)#
g <- add.edges(g, c(0,5, 0,10, 5, 10))#
wtc <- walktrap.community(g)#
memb <- community.to.membership(g, wtc$merges, steps=12)#
modularity(g, memb$membership)
quit()
demo()
demo(image)
demo(graphics)
quartz()
x <- rand(1,100)
x <- randn(1,100)
x <- rnorm(1,100)
x
help(rnorm)
x <- rnorm(100)
x
y <- rnorm(100)
plot(x,y)
par <- opar
plot(x,y)
par <- opar
x <- rnorm(100)
y <- rnorm(100)
plot(x,y)
quit()
help("Deprecated")
help(find.package)
quit()
quit(0
)
quit()
data <- read.delim("landdata_states.txt")
head(data)
tail(data)
require(ggplot2)#
require(reshape)
require(reshape2)
df <- subset(data,STATE=="CT"|STATE=="VT"|STATE=="NH"|STATE=="KY")
p <- ggplot(df, aes(x=Date,y=Land.Value,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE)) +  scale_shape_manual(value = c(15, 16)) + theme_bw() + opts(aspect.ratio=1) + scale_colour_brewer(palette="Set1") #+ opts(axis.text.x=theme_text(angle=-90, hjust=0))
p <- ggplot(df, aes(x=Date,y=Land.Value,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE))
is.factor(df$Date)
is.date(df$Date)
data <- read.delim("landdata_states.txt")
is.factor(df$Date)
df <- subset(data,STATE=="CT"|STATE=="VT"|STATE=="NH"|STATE=="KY")
is.factor(df$Date)
p <- ggplot(df, aes(x=Date,y=Land.Value.usd,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE)) +  scale_shape_manual(value = c(15, 16)) + theme_bw()
p <- ggplot(df, aes(x=Date,y=Land.Value.usd,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE))
is.factor(df$Land.Value.usd)
head(df)
is.factor(df$Date)
is.factor(df$Home.Value.usd)
as.numeric(df$Home.Value.usd)
as.vector(df$Home.Value.usd)
help(as)
as.integer(df$Home.Value.usd)
data <- read.delim("landdata_states.txt")
df <- subset(data,STATE=="CT"|STATE=="VT"|STATE=="NH"|STATE=="KY")
is.factor(df$Date)
is.factor(df$Home.Value.usd)
as.integer(df$Home.Value.usd)
is.integer(df$Home.Value.usd)
is.numeric(df$Home.Value.usd)
p <- ggplot(df, aes(x=Date,y=Land.Value.usd,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE))
is.numeric(df$Land.Value.usd)
is.factor(df$Land.Value.usd)
tail(df)
as.integer(df$Land.Value.usd)
as.integer(data$Land.Value.usd)
tail(data)
levels(df$Land.Value.usd)
df$Land.Value.usd
levels(df$Land.Value.usd)
df$Land.Value.usd
tail(df)
df$Land.Value.usd <- as.numeric(level(df$Land.Value.usd))
df$Land.Value.usd <- as.numeric(levels(df$Land.Value.usd))
tail(df)
p <- ggplot(df, aes(x=Date,y=Land.Value.usd,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE))
is.factor(df$Land.Value.usd)
as.character(df$Land.Value.usd)
as.numeric(as.character(df$Land.Value.usd))
is.factor(as.numeric(as.character(df$Land.Value.usd)))
is.vector(as.numeric(as.character(df$Land.Value.usd)))
df <- subset(data,STATE=="CT"|STATE=="VT"|STATE=="NH"|STATE=="KY")#
df$Land.Value.usd <- as.numeric(as.character(df$Land.Value.usd))
p <- ggplot(df, aes(x=Date,y=Land.Value.usd,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE))
data <- read.delim("landdata_states.txt")   #first removed the " $ % , symbols from the dataset and changed headers
df <- subset(data,STATE=="CT"|STATE=="VT"|STATE=="NH"|STATE=="KY")
p <- ggplot(df, aes(x=Date,y=Land.Value.usd,group=STATE)) #
p + geom_point(aes(color=STATE,shape=STATE)) + geom_line(aes(color=STATE))
p <- ggplot(df, aes(x=Date,y=Land.Value.usd,group=STATE)) #
p + geom_line(aes(color=STATE))
p + geom_line(aes(color=STATE)) + theme_bw() + opts(aspect.ratio=1) + scale_colour_brewer(palette="Set1")
ggsave(file="geomline_landvalue-time-STATE.pdf")
p <- ggplot(df, aes(x=Date,y=Land.Share.frac*100,group=STATE)) #
p + geom_line(aes(color=STATE)) + theme_bw() + opts(aspect.ratio=1) + scale_colour_brewer(palette="Set1")
ggsave(file="geomline_landshare-time-STATE.pdf")
summary(cars)
plot(cars)
library(devtools)
require(googleVis)
cars
asis
install.packages('googleVis')
require(googleVis)
demo(googleVis)
fruits
help(gvisMotionChart)
M1
Fruits
Volcano
volcano
quit(_
quit()
actvFraction <- c(0.2838, 0.31114, 0.14775, 0.058811, 0.104, 0.079014, 0.015483)
motorState <- c("quiet","quiet","quiet","active","active","active","active")
df <-  data.frame(actvFraction, motorState)
df
summary(df)
help(summary)
table(df)
lm.df <- with(df, lm(actvFraction ~ motorState))
summary(lm.df)
sum(df$actvFraction)
is.factor(df$motorState)
help(sum)
help(pylr)
help.start(0
help.start(0)
help.start()
require(plyr)
dfx <- data.frame(#
  group = c(rep('A', 8), rep('B', 15), rep('C', 6)),#
  sex = sample(c("M", "F"), size = 29, replace = TRUE),#
  age = runif(n = 29, min = 18, max = 54)#
)
dfx
ddply(dfx, .(group, sex), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2))
ddply(dfx, .(group, sex), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2), len = length(age))
help(summarize)
ddply(dfx, (group, sex), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2), len = length(age))
ddply(dfx, .(group, sex), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2))
ddply(dfx, c("group", "sex"), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2))
ddply(dfx, c("group", "sex"), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2), N = length(age))
dfx
sd(1)
sd(c(1,2))
help(sd)
ddply(dfx, c("group", "sex"), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2), N = length(age), se = sd/sqrt(N))
ddply(dfx, c("group", "sex"), summarize,#
 mean = round(mean(age), 2),#
 sd = round(sd(age), 2), N = sum(!is.na(age)), se = sd/sqrt(N))
help(summarize)
ddply(df, c("motorState"), summarize,#
sum = round(sum(actvFraction)),#
mean = round(mean(actvFraction), 2),#
sd = round(sd(actvFraction), 2), #
N = sum(!is.na(age)), #
se = sd/sqrt(N))
ddply(df, c("motorState"), summarize,#
sum = round(sum(actvFraction)),#
mean = round(mean(actvFraction), 2),#
sd = round(sd(actvFraction), 2), #
N = sum(!is.na(actvFraction)), #
se = sd/sqrt(N))
ddply(df, c("motorState"), summarize,#
sum = round(sum(actvFraction),2),#
mean = round(mean(actvFraction), 2),#
sd = round(sd(actvFraction), 2), #
N = sum(!is.na(actvFraction)), #
se = sd/sqrt(N))
ddply(df, c("motorState"), summarize,#
sum = round(sum(actvFraction),2),#
mean = round(mean(actvFraction), 2),#
sd = round(sd(actvFraction), 2), #
N = length(actvFraction), #
se = sd/sqrt(N))
help(aggregate)
with(df,aggregate(actvFraction,by=motorState,mean))
with(df,aggregate(actvFraction ~ motorState,mean))
aggregate(actvFraction ~ motorState,df,mean)
aggregate(actvFraction ~ motorState,df,sum)
myFormula <- (actvFraction ~ motorState)#
#
	ddply(df, myFormula, summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction), 2),#
	sd = round(sd(actvFraction), 2), #
	N = length(actvFraction), #
	se = sd/sqrt(N))
ddply(df, ~ motorState, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(df, actvFraction ~ motorState, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
dfx
ddply(df, .(actvFraction ~ motorState), summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(df, actvFraction ~ motorState, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(df, ~ motorState, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(dfx, ~ group, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(dfx, sex ~ group, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(dfx, age ~ group, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(dfx, ~ group + age, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(dfx, ~ group + sez, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(dfx, ~ group + sex, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(dfx, ~ group + sex, summarize,#
	sum = round(sum(age),2),#
	N = length(age))
ddply(dfx, sex ~ group, summarize,#
	sum = round(sum(age),2),#
	N = length(age))
as.quoted(sex,group)
myFormula <- ~ group + sex
myFormula
ddply(dfx, myFormula, summarize,#
	sum = round(sum(age),2),#
	N = length(age))
myFormula <- ~ motorState#
	ddply(df, myFormula, summarize,#
	sum = round(sum(actvFraction),2),#
	N = length(actvFraction))
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction), 2),#
	med = round(median(actvFraction),2),#
	std = round(sd(actvFraction), 2), #
	N = sum(!is.na(actvFraction), #
	sem = std/sqrt(N),#
	CI95 = qnorm(0.975)*sem,#
	medSD = round(mad(actvFraction)),#
	seMed<-1.25*sem)
)
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction), 2),#
	med = round(median(actvFraction),2),#
	std = round(sd(actvFraction), 2), #
	N = sum(!is.na(actvFraction), #
	sem = std/sqrt(N),#
	CI95 = qnorm(0.975)*sem,#
	medSD = round(mad(actvFraction),2),#
	seMed = 1.25*sem)
))
help(mad)
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction),2),#
	med = round(median(actvFraction),2),#
	std = round(sd(actvFraction),2), #
	N = sum(!is.na(actvFraction), #
	sem = std/sqrt(N))
)
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction), 2),#
	sd = round(sd(actvFraction), 2), #
	N = length(actvFraction), #
	se = sd/sqrt(N))
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction), 2),#
	sd = round(sd(actvFraction), 2), #
	N = length(actvFraction), #
	sem = sd/sqrt(N))
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction), 2),#
	std = round(sd(actvFraction), 2), #
	N = length(actvFraction), #
	se = std/sqrt(N))
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction),2),#
	med = round(median(actvFraction),2),#
	std = round(sd(actvFraction),2), #
	N = sum(!is.na(actvFraction)), #
	sem = std/sqrt(N),#
	CI95 = qnorm(0.975)*sem,#
	medSD = round(mad(actvFraction),2),#
	seMed = 1.25*sem)
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction),2),#
	std = round(sd(actvFraction),2), #
	N = sum(!is.na(actvFraction)), #
	sem = std/sqrt(N),#
	CI95 = qnorm(0.975)*sem,#
	median = round(median(actvFraction),2),#
	medSD = round(mad(actvFraction),2),#
	seMed = 1.25*sem)
all.vars(myFormula)
ddply(df, c("motorState"), summarize,#
	sum = round(sum(actvFraction),2),#
	mean = round(mean(actvFraction),2),#
	std = round(sd(actvFraction),2), #
	N = sum(!is.na(actvFraction)), #
	sem = std/sqrt(N),#
	CI95 = qnorm(0.975)*sem,#
	median = round(median(actvFraction),2),#
	medSD = round(mad(actvFraction),2),#
	seMed = 1.25*sem)
help(t.test)
measVar = "actvFraction"#
groupVars = c("motorState")
myFormula <- as.formula(paste(measVar, paste(groupVars, collapse=" + "), sep=" ~ "))
myFormula
t.test(myFormula,df)
help(t.test)
t.test(myFormula,df,alternative="less")
t.test(myFormula,df,alternative="greater")
printSummary <- function(data=NULL, measVar, groupVars=NULL, na.rm=FALSE,#
                      conf.interval=.95, .drop=TRUE) {#
	require(plyr)#
#
	#Handle NAs#
	len <- function (x, na.rm=FALSE) {#
	    if (na.rm) sum(!is.na(x))#
	    else       length(x)#
	}#
#
	#Main summary#
	data2 <- ddply(data, groupVars, .drop=.drop,#
	      .fun = function(xx, col) {#
	        c(N    = len(xx[[col]], na.rm=na.rm),#
	          mean = mean   (xx[[col]], na.rm=na.rm),#
	          sd   = sd     (xx[[col]], na.rm=na.rm)#
	        )#
	      },#
	      measVar#
	    )#
#
	#Standard error of the mean#
	data2$se <- data2$sd / sqrt(data2$N)#
#
	#Confidence interval of 95%#
	data2$CI95 = qnorm(0.975)*se#
#
	return(data2)#
}
printSummary(df, measVar,groupVars)
printSummary <- function(data=NULL, measVar, groupVars=NULL, na.rm=FALSE,#
                      conf.interval=.95, .drop=TRUE) {#
	require(plyr)#
#
	#Handle NAs#
	len <- function (x, na.rm=FALSE) {#
	    if (na.rm) sum(!is.na(x))#
	    else       length(x)#
	}#
#
	#Main summary#
	data2 <- ddply(data, groupVars, .drop=.drop,#
	      .fun = function(xx, col) {#
	        c(N    = len(xx[[col]], na.rm=na.rm),#
	          mean = mean   (xx[[col]], na.rm=na.rm),#
	          sd   = sd     (xx[[col]], na.rm=na.rm)#
	        )#
	      },#
	      measVar#
	    )#
#
	#Standard error of the mean#
	data2$se <- data2$sd / sqrt(data2$N)#
#
	#Confidence interval of 95%#
	data2$CI95 = qnorm(0.975)*data2$se#
#
	return(data2)#
}
printSummary(df, measVar,groupVars)
help(wilcox.test)
wilcox.test(myFormula,df)
printSummary <- function(data=NULL, measVar, groupVars=NULL, na.rm=FALSE,#
	                      conf.interval=.95, .drop=TRUE) {#
		require(plyr)#
		#Handle NAs#
		len <- function (x, na.rm=FALSE) {#
		    if (na.rm) sum(!is.na(x))#
		    else       length(x)#
		}#
		#Main summary#
		data2 <- ddply(data, groupVars, .drop=.drop,#
		      .fun = function(xx, col) {#
		        c(N    = len(xx[[col]], na.rm=na.rm),#
		          sum = sum   (xx[[col]], na.rm=na.rm),#
		          mean = mean   (xx[[col]], na.rm=na.rm),#
		          sd   = sd     (xx[[col]], na.rm=na.rm)#
		        )#
		      },#
		      measVar#
		    )#
		#Standard error of the mean#
		data2$se <- data2$sd / sqrt(data2$N)#
		#Confidence interval of 95%#
		data2$CI95 = qnorm(0.975)*data2$se#
		return(data2)#
	}
printSummary(df, measVar,groupVars)
quit()
require(ggplot2)
mag <- c(5.0,2.5,2.0,1.0,0.6)#
umperpx <- c(2.270,4.525,5.850,11.350,19.178)#
data <- data.frame(mag,umperpx)
data
require(ggplot2)#
#The following shows that the relationship between objective magnification and pixel dimensions clearly follows a power equation relationship, a log-log relationship which follows the form y = bx^m#
qplot(log(mag),log(umperpx),data=data)#
qplot(mag,log(umperpx),data=data)#
qplot(mag,umperpx,data=data)
qplot(log(mag),log(umperpx),data=data)
qplot(mag,log(umperpx),data=data)
qplot(mag,umperpx,data=data)
qplot(mag,log(umperpx),data=data)
qplot(log(mag),log(umperpx),data=data)
lm.fit <- lm(umperpx ~ mag,data=data)
lm.fit <- lm(log(umperpx) ~ log(mag),data=data)
lm.fit
summary(lm.fit)
with(data, plot(mag, umperpx, log="xy"))#
lines(seq(0.5,5,by=0.5), exp(predict(lm.fit, data.frame(mag = seq(0.5,5,by=0.5)))))
with(data, plot(mag, umperpx))#
lines(seq(0.5,5,by=0.5), predict(lm.fit, data.frame(mag = seq(0.5,5,by=0.5))))
x=seq(0.5,5,by=0.01)#
#plot(data,pch=22)#
b = exp(2.439593)   #the intercept estimate was predicted based on the log-transformed data, this must be converted back to the normal space with natural exponent exp()#
m = -1.004096
fitted.data <- data.frame(x = x, y = b*x^m)
ggplot(data, aes(x = mag, y = umperpx)) + geom_line(data = fitted.data, aes(x = x, y = y), colour = "red") + geom_point() + title(main = "y=11.46837*x^-1.004096")
quartz();#
p <- ggplot(data, aes(x = mag, y = umperpx))#
p + geom_line(data = fitted.data, aes(x = x, y = y), colour = "red") + geom_point()
p + geom_line(data = fitted.data, aes(x = x, y = y), colour = "red") + geom_point() + title(main = "y=11.46837*x^-1.004096")
help.start()
p + geom_line(data = fitted.data, aes(x = x, y = y), colour = "red") + geom_point() + ggtitle(main = "y=11.46837*x^-1.004096")
ggtitle("y=11.46837*x^-1.004096")
p + geom_line(data = fitted.data, aes(x = x, y = y), colour = "red") + geom_point() + ggtitle("y=11.46837*x^-1.004096")
x
cbind(x,y)
fitted.data
ggsave(file="obj-pixeldim-power-law-regression")
ggsave(file="obj-pixeldim-power-law-regression.pdf")
quit()
time <- c(0, 0.200, 1.2, 2.2, 5.6)
frame <-  c(1,2,7,12,29)
marker <- c(0,0.001,1.207,2.212,5.628)
plot(time,frame)
help(diff)
marker - time
marker <- c(0,0.201,1.207,2.212,5.628)
marker - time
plot(time,marker)
plot(frame,time-marker)
plot(frame,marker-time)
data = data.frame(time,frame,marker)
lm.fit <- lm(marker-time ~ frame,data=data)
lm.fit
summary(lm.fit)
lm.fit
b = -0.0005725; m = 0.0009973;
b
m
(0.200 - b)/m
(0.200)/m
200*0.2
time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)
plot(frame,marker-time)
data = data.frame(time,frame,marker)#
	lm.fit <- lm(marker-time ~ frame,data=data)#
	lm.fit
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	lm.fit <- lm(marker-time ~ frame,data=data)#
 	coef(lm.fit)
coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	data = data.frame(time,frame,marker)#
	lm.fit <- lm(marker-time ~ frame,data=data)#
 	coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
data <- data.frame(time,frame,marker)#
	lm.fit <- lm(marker-time ~ frame,data=data)#
 	coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)#
	data <- data.frame(time,frame,marker)#
 	lm.fit <- lm(marker-time ~ frame,data=data)#
 	lm.fit
coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	data <- data.frame(time,frame,marker)#
	data$diff <- marker-time#
	lm.fit <- lm(diff ~ frame,data=data)#
 	coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
quartz()
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(frame,marker-time)
plot(time,marker-time)
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time
21*0.2
21*5
20.2*5
time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)#
	marker - time
time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)#
	marker - time#
	df <- data.frame(time,frame,marker)#
 	lm.fit <- lm(marker-time ~ frame,data=df)#
 	lm.fit#
	coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
quartz();#
	frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(time,marker-time)#
	df <- data.frame(time,frame,marker)#
	data$diff <- marker-time#
	lm.fit <- lm(diff ~ frame,data=df)#
 	coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
quartz();#
	frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(time,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lm.fit <- lm(diff ~ frame,data=df)#
 	coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(time,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	x <- coef(lmfit)["frame"] #
	(0.200-b)/m
quartz();#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)#
	marker - time#
	plot(frame, marker-time)#
	df <- data.frame(time,frame,marker)#
 	lm.fit <- lm(marker-time ~ frame,data=df)#
 	lm.fit#
	coef(lm.fit) #
	b = coef(lm.fit)["(Intercept)"] #
	x = coef(lm.fit)["frame"] #
	(0.200-b)/m
quartz();#
	frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	x <- coef(lmfit)["frame"] #
	(0.200-b)/m
rm(lm.fit)
rm(lmfit)
time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)#
	marker - time#
	plot(frame, marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
 	lm.fit <- lm(diff ~ frame,data=df)#
 	lm.fit#
	coef(lm.fit) #
	b <- coef(lm.fit)["(Intercept)"] #
	x <- coef(lm.fit)["frame"] #
	(0.200-b)/m
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	x <- coef(lmfit)["frame"] #
	(0.200-b)/m
df
time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)#
	marker - time#
	plot(frame, marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time
df
(0.200-b)/m
b
m
(0.200)/m
b
is.numeric(b)
quartz();#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	frame <-  c(1,2,6,26,101)#
	marker <- c(0.0005,0.2015,1.0055,5.0255,20.1004)#
	marker - time#
	plot(frame, marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
 	lm.fit <- lm(diff ~ frame,data=df)#
 	lm.fit#
	coef(lm.fit) #
	b <- coef(lm.fit)["(Intercept)"] #
	m <- coef(lm.fit)["frame"] #
	(0.200-b)/m
quartz();#
	frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	m <- coef(lmfit)["frame"] #
	(0.200-b)/m
10000/200
coef(lmfit)
1/600
1/6100
1/100
1/30
1/60
1/300
1/600
quartz();#
	frame <-  c(1,2,6,26,101)#
	time <- c(0.0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.0,0.2012,1.0016,5.0034,20.0102)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	m <- coef(lmfit)["frame"] #
	(0.200-b)/m
frame <-  c(2,6,26,101)#
	time <- c(0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.2012,1.0016,5.0034,20.0102)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	m <- coef(lmfit)["frame"] #
	(0.200)/m
(0.2-b)/m
1100*0.2
quartz();#
	frame <-  c(2,6,26,101)#
	time <- c(0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.2012,1.0016,5.0034,20.0102)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	m <- coef(lmfit)["frame"] #
	(0.200-b)/m
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	m <- coef(lmfit)["frame"] #
	(0.200-b)/m
m*3000 + b
m*3000
b
m
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.00002,0.20004,1.00012,5.00052,20.00202)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	m <- coef(lmfit)["frame"] #
	(0.200-b)/m
m
m*3000
m*5000
m*3000 + b
frame <-  c(1,2,6,26,101)#
	time <- c(0, 0.200, 1.0, 5.0, 20.0)#
	marker <- c(0.0,0.200,0.999960,4.9998,19.999240)#
	marker-time#
	plot(frame,marker-time)#
	df <- data.frame(time,frame,marker)#
	df$diff <- marker-time#
	lmfit <- lm(diff ~ frame,data=df)#
 	coef(lmfit) #
	b <- coef(lmfit)["(Intercept)"] #
	m <- coef(lmfit)["frame"] #
	(0.200-b)/m
m*3000 + b
quit();
help(pearson)
??pearson
help(cor)
quit()
edgelist<-read.delim('/Users/ackman/Data/2photon/131208/2014-01-07-003602/dCorr.txt')
edgelist
require(pylr)#
df <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue = mean(rvalue),#
sd = sd(rvalue), #
N = length(rvalue), #
se = sd/sqrt(N))
require(plyr)#
df <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue = mean(rvalue),#
sd = sd(rvalue), #
N = length(rvalue), #
se = sd/sqrt(N))
df
df <- ddply(edgelist, c("node1","node2"), summarize,#
rvalueMean = mean(rvalue),#
sd = sd(rvalue), #
N = length(rvalue), #
se = sd/sqrt(N))
df
df <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = sd/sqrt(N))
df <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))
df
colnames(df)
colnames(df)['rvalue.mean']
colnames(df) == 'rvalue.mean'
colnames(df)[['rvalue.mean']]
colnames(df)
colnames(df) == 'rvalue.mean'
colnames(df)[colnames(df) == 'rvalue.mean']
colnames(df)[colnames(df) == 'rvalue.mean'] <- 'rvalue'
colname(df)
colnames(df)
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(df,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
load(plyr)
help(require)
library(plyr)
library(igraph)#
library(RColorBrewer)
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(df,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
rthresh <- 0.2#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(df,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
edgelist<-read.delim('/Users/ackman/Data/2photon/120518i/2014-01-03-231550/dCorr.txt'
)
edgelist<-read.delim('/Users/ackman/Data/2photon/120518i/2014-01-03-231550/dCorr.txt')
df <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(df)[colnames(df) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.2#
fnm <- '120518'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(df,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
rthresh <- 0.1#
fnm <- '120518'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(df,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
summary(g)
print(g)
print(fastGreeyCom)
print(fastGreedyCom)
print(fastgreedyCom)
degree(g)#
degree.distribution(g)#
degree.distribution(g,cumulative = TRUE)
average.path.length(g)
library(ggplot2)
df <- data.frame(degree(g))#
colnames(df) <- c("degree")#
p <- ggplot(df, aes(x=degree)) + xlab("degree") + theme_bw()#
p + geom_histogram(binwidth = 2) + scale_colour_brewer(palette="Set1") + opts(aspect.ratio=1)  #raw counts#
ggsave(file=paste("120518_07-degreeDist", format(Sys.time(),"%y%m%d-%H%M%S"), ".pdf",sep=""))
g <- barabasi.game(1000) # increase this number to have a better estimate#
d <- degree(g, mode="in")#
fit1 <- power.law.fit(d+1, 10)#
fit2 <- power.law.fit(d+1, 10, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
g
d
df <- data.frame(degree(g))#
colnames(df) <- c("degree")#
p <- ggplot(df, aes(x=degree)) + xlab("degree") + theme_bw()#
p + geom_histogram(binwidth = 2) + scale_colour_brewer(palette="Set1") + opts(aspect.ratio=1)  #raw counts
This should approximately yield the correct exponent 3#
g <- barabasi.game(1000) # increase this number to have a better estimate#
d <- degree(g, mode="in")#
fit1 <- power.law.fit(d+1, 10)#
fit2 <- power.law.fit(d+1, 10, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)#
#
df <- data.frame(degree(g))#
colnames(df) <- c("degree")#
p <- ggplot(df, aes(x=degree)) + xlab("degree") + theme_bw()#
p + geom_histogram(binwidth = 2) + scale_colour_brewer(palette="Set1") + opts(aspect.ratio=1)  #raw counts#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
title("barabasi.game(1000), powerlaw")#
quartz.save(file=paste(dateStr, "-degreeDist-", "barabasiGame-powerlaw",  ".pdf",sep=""), type = "pdf")
ggsave(file=paste(dateStr, "-degreeDist-", "barabasiGame-powerlaw",  ".pdf",sep=""))
This should approximately yield the correct exponent 3#
g <- barabasi.game(33) # increase this number to have a better estimate#
d <- degree(g, mode="in")#
fit1 <- power.law.fit(d+1, 10)#
fit2 <- power.law.fit(d+1, 10, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)#
#
df <- data.frame(degree(g))#
colnames(df) <- c("degree")#
p <- ggplot(df, aes(x=degree)) + xlab("degree") + theme_bw()#
p + geom_histogram(binwidth = 2) + scale_colour_brewer(palette="Set1") + opts(aspect.ratio=1)  #raw counts#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
title("barabasi.game(33), powerlaw")#
ggsave(file=paste(dateStr, "-degreeDist-", "barabasiGame-powerlaw",  ".pdf",sep=""))
fit1
fit2
d
This should approximately yield the correct exponent 3#
g <- barabasi.game(1000) # increase this number to have a better estimate#
d <- degree(g, mode="in")#
fit1 <- power.law.fit(d+1, 10)#
fit2 <- power.law.fit(d+1, 10, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
fit2
fit2
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
print(fastgreedyCom)#
degree(g)#
degree.distribution(g)#
degree.distribution(g,cumulative = TRUE)#
average.path.length(g) #
diameter(g)
centrality(g)
hub.score(g)
hub.score(g)$vector
------Histogram of degree distribution-------------------------------------------------------------#
df <- data.frame(degree(g))#
colnames(df) <- c("degree")#
p <- ggplot(df, aes(x=degree)) + xlab("degree") + theme_bw()#
p + geom_histogram(binwidth = 2) + scale_colour_brewer(palette="Set1") + opts(aspect.ratio=1)  #raw counts#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
ggsave(file=paste(dateStr, "-degreeDist-", fnm,  ".pdf",sep=""))
edgelist<-read.delim('/Users/ackman/Data/2photon/131208/2014-01-07-003602/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
print(fastgreedyCom)#
degree(g)#
degree.distribution(g)#
degree.distribution(g,cumulative = TRUE)#
average.path.length(g) #
diameter(g)#
hub.score(g)$vector
hub.score(g)
------Histogram of degree distribution-------------------------------------------------------------#
df <- data.frame(degree(g))#
colnames(df) <- c("degree")#
p <- ggplot(df, aes(x=degree)) + xlab("degree") + theme_bw()#
p + geom_histogram(binwidth = 2) + scale_colour_brewer(palette="Set1") + opts(aspect.ratio=1)  #raw counts#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
ggsave(file=paste(dateStr, "-degreeDist-", fnm,  ".pdf",sep=""))
help(degree)
g <- graph.ring(10)#
degree(g)#
g2 <- erdos.renyi.game(1000, 10/1000)#
degree.distribution(g2)
plot(degree.distribution(g2))
plot(degree.distribution(g2))
quartz;plot(degree.distribution(g2))
quartz; plot(degree.distribution(g2))
edgelist<-read.delim('/Users/ackman/Data/2photon/131208/2014-01-07-003602/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
quartz()
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
print(fastgreedyCom)#
degree(g)#
degree.distribution(g)#
degree.distribution(g,cumulative = TRUE)#
average.path.length(g) #
diameter(g)#
hub.score(g)$vector
mean(degree(g))
d <- degree(g)#
fit1 <- power.law.fit(d+1, 10)#
fit2 <- power.law.fit(d+1, 10, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
fit2
hist(degree(g))
help(power.law.fit)
d <- degree(g)#
fit1 <- power.law.fit(d+1)#
fit2 <- power.law.fit(d+1, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d)#
fit2 <- power.law.fit(d, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d
d <- degree(g)#
fit1 <- power.law.fit(d,2)#
fit2 <- power.law.fit(d,2, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d,3)#
fit2 <- power.law.fit(d,3, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d,4)#
fit2 <- power.law.fit(d,4, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d,5)#
fit2 <- power.law.fit(d,5, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d,8)#
fit2 <- power.law.fit(d,8, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d,7)#
fit2 <- power.law.fit(d,7, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d,6)#
fit2 <- power.law.fit(d,6, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d)
fit1$alpha
fit1$xmin
fit1
edgelist<-read.delim('/Users/ackman/Data/2photon/120518i/2014-01-03-231550/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
edgelist<-read.delim('/Users/ackman/Data/2photon/120518i/2014-01-03-231550/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
print(fastgreedyCom)#
degree(g)#
degree.distribution(g)#
degree.distribution(g,cumulative = TRUE)#
average.path.length(g) #
diameter(g)#
hub.score(g)$vector
mean(degree(g))
d <- degree(g)#
fit1 <- power.law.fit(d)
fit1
d <- degree(g)#
fit1 <- power.law.fit(d,7)#
fit2 <- power.law.fit(d,7, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
d <- degree(g)#
fit1 <- power.law.fit(d,3)#
fit2 <- power.law.fit(d,3, implementation="R.mle")#
#
fit1$alpha#
coef(fit2)#
fit1$logLik#
logLik(fit2)
edgelist<-read.delim('/Users/ackman/Data/2photon/131208/2014-01-07-003602/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
quartz()
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.1#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
rthresh <- 0.15#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
edgelist<-read.delim('/Users/ackman/Data/2photon/120518i/2014-01-03-231550/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.15#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
rthresh <- 0.15#
fnm <- 'younger'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
edgelist<-read.delim('/Users/ackman/Data/2photon/131208/2014-01-07-003602/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.15#
fnm <- 'P8'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")
edgelist<-read.delim('/Users/ackman/Data/2photon/131208/2014-01-07-003602/dCorr.txt')
help(for)
help()
rthresh <- 0.2#
fnm <- '131208_01'#
fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
d3 <- subset(edgelist,filename==fnm2)#
d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d4,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
# quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
for(i in c('131208_01','131208_03','131208_04','131208_05')) {#
for(j in c(0.1,0.15,0.2)) {#
rthresh <- j#
fnm <- '131208_03'#
fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
d3 <- subset(edgelist,filename==fnm2)#
d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d4,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
# quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")#
}#
}
for(j in c(0.1,0.15,0.2)) {#
for(i in c('131208_01','131208_03','131208_04','131208_05')) {#
rthresh <- j#
fnm <- i#
fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
d3 <- subset(edgelist,filename==fnm2)#
d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d4,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
# quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
# quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")#
}#
}
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.15#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
edgelist<-read.delim('/Users/ackman/Data/2photon/120518i/2014-01-03-231550/dCorr.txt')
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.15#
fnm <- '131208'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
d2 <- ddply(edgelist, c("node1","node2"), summarize,#
rvalue.mean = mean(rvalue),#
rvalue.sd = sd(rvalue), #
N = length(rvalue), #
rvalue.sem = rvalue.sd/sqrt(N))#
colnames(d2)[colnames(d2) == 'rvalue.mean'] <- 'rvalue'#
#
rthresh <- 0.15#
fnm <- '120518'#
# fnm2 <- paste(fnm,".tif",sep="")#
lo <- 'layout.fruchterman.reingold'#
# lo <- 'layout.kamada.kawai'#
# lo <- 'layout.lgl'#
# d3 <- subset(edgelist,filename==fnm2)#
# d4 <- with(d3,data.frame(node1,node2,rvalue))#
edgelist2<-subset(d2,rvalue > rthresh)#
g <- graph.data.frame(edgelist2, directed=FALSE)#
E(g)$weight <- E(g)$rvalue#
E(g)$width <- 1#
E(g)[ weight >= 0.3 ]$width <- 3#
E(g)[ weight >= 0.5 ]$width <- 5#
fastgreedyCom<-fastgreedy.community(g,weights=E(g)$weight)#
V(g)$color <- fastgreedyCom$membership#
quartz();#
# palette(rainbow(max(V(g)$color),alpha=0.5))#
mypalette <- adjustcolor(brewer.pal(max(V(g)$color),"Set1"),0.6)#
palette(mypalette)#
plot(g, layout=eval(parse(text=lo)), edge.width=E(g)$width, edge.color="black", vertex.label.color="black")#
# palette("default")#
title(paste(fnm,', fastgreedy default, ', lo, 'r>', rthresh))#
dateStr=format(Sys.time(),"%y%m%d-%H%M%S")#
quartz.save(file=paste(dateStr, "-", fnm,  ".png",sep=""), type = "png", dpi=150)#
quartz.save(file=paste(dateStr, "-", fnm,  ".pdf",sep=""), type = "pdf")
quit()