all-analyses.Rmd

Analysis for Cornish et al. (submitted)
============================================================================================

"In order to determine whether string sets are being acquired more faithfully over time, we calculated error between generations in terms of the normalized edit distance"

```{r}
# setwd('/Users/rickdale/Dropbox/projects/studies/rickmorten/cornishetal/paper/newstuff/newdata/')
set.seed(2) # stable randomization seed for statistics in paper; patterns hold in all seeds explored
library(lme4)
mindists = matrix(nrow=8,ncol=10)
numzeros = matrix(nrow=8,ncol=10)
avglen = matrix(nrow=8,ncol=10)
for (chain in 1:8) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')
  for (generation in 2:11) {
    dists = adist(strs[,generation-1],strs[,generation])
    lens = matrix(nchar(as.character(strs[,generation-1])),ncol=15,nrow=15)
    dists = dists/lens # normalized by length
    mindists[chain,generation-1] = mean(apply(dists,1,min)) # simply takes first minimum
    numzeros[chain,generation-1] = sum(apply(dists,1,function(x) {x==0})) # count # of totally correct ones
    avglen[chain,generation-1] = mean(nchar(as.character(strs[,generation])))
  }
}
shft = .1
ms = colMeans(mindists)
sds = apply(mindists,2,sd)
plot(ms,type='b',xlab='Generation',
     ylab='Minimal normalized Levensthein (error)',
     ylim=c(0,.45),main='Error over generations')
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
regdata = cbind(which(mindists>-1,arr.ind=T),mindists[which(mindists>-1,arr.ind=T)])
lmo = lmer(regdata[,3]~regdata[,2]+(1+regdata[,2]|regdata[,1]))
coefs = data.frame(summary(lmo)$coefficients)
coefs$p = 2*(1-pnorm(abs(coefs$t.value)));coefs
t.test(mindists[,1]-mindists[,10])
mean(mindists[,1]);sd(mindists[,1])
mean(mindists[,10]);sd(mindists[,10])
```

"The boost in overall accuracy translates into a significant increase in the number of correctly recalled items"

```{r}
ms = colMeans(numzeros)
sds = apply(numzeros,2,sd)
plot(ms,type='b',xlab='Generation',ylab='Number of exact reproductions',
     ylim=c(0,9),main='Correct strings over generations')
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
regdata = cbind(which(numzeros>-1,arr.ind=T),numzeros[which(numzeros>-1,arr.ind=T)])
lmo = lmer(regdata[,3]~regdata[,2]+(1+regdata[,2]|regdata[,1]))
coefs = data.frame(summary(lmo)$coefficients)
coefs$p = 2*(1-pnorm(abs(coefs$t.value)));coefs
t.test(numzeros[,1]-numzeros[,10])
mean(numzeros[,1]);sd(numzeros[,1])
mean(numzeros[,10]);sd(numzeros[,10])
```

"Importantly, the improved learnability did not come at the cost of a collapse of the string sets into very short sequences"

```{r}
ms = colMeans(avglen)
sds = apply(avglen,2,sd)
plot(ms,type='b',xlab='Generation',ylab='Average string length',
     ylim=c(0,4.5),main='Average length of strings')
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
regdata = cbind(which(avglen>-1,arr.ind=T),avglen[which(avglen>-1,arr.ind=T)])
lmer(regdata[,3]~regdata[,2]+(1+regdata[,2]|regdata[,1]))
lmo = lmer(regdata[,3]~regdata[,2]+(1+regdata[,2]|regdata[,1]))
coefs = data.frame(summary(lmo)$coefficients)
coefs$p = 2*(1-pnorm(abs(coefs$t.value)));coefs
t.test(avglen[,1]-avglen[,10])
mean(avglen[,1]);sd(avglen[,1])
mean(avglen[,10]);sd(avglen[,10])
```

"Figure 2 indicates that the amount of reuse of chunks (structure) does increase over time"

```{r}
bitrifound = function(s1) { 
  bigs = unlist(strsplit(s1,""))
  bigs = unique(cbind(bigs[1:(length(bigs)-1)],bigs[2:length(bigs)]))
  nbigs = 0;totbigs = data.frame()
  for (s2 in strings) { # loop through all strings, comparing bigrams... determine # of shared bigrams across all
    s2 = unlist(strsplit(s2,""))
    s2 = unique(cbind(s2[1:(length(s2)-1)],s2[2:length(s2)]))
    all = rbind(bigs,s2)
    nbigs = nbigs + dim(all)[1] - dim(unique(all))[1] # total bigrams shared
    totbigs = rbind(totbigs,s2)
  }
  ntris = 0;tottris = data.frame()
  if (nchar(s1)>=3) {
    tris = unlist(strsplit(s1,""))
    tris = unique(cbind(tris[1:(length(tris)-2)],tris[2:(length(tris)-1)],tris[3:length(tris)]))
    for (s2 in strings[nchar(strings)>=3]) { # loop through only those that had trigrams 
      s2 = unlist(strsplit(s2,""))
      s2 = unique(cbind(s2[1:(length(s2)-2)],s2[2:(length(s2)-1)],s2[3:length(s2)]))
      all = rbind(tris,s2)
      ntris = ntris + dim(all)[1] - dim(unique(all))[1] # total trigrams shared
      tottris = rbind(tottris,s2)
    }
  }  
  return(ntris/dim(unique(tottris))[1]+nbigs/dim(unique(totbigs))[1])
  #return(ntris+nbigs)
}

ACS = matrix(nrow=8,ncol=10)
for (chain in 1:8) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')
  for (generation in 2:11) {
    strings <<- as.character(strs[,generation-1])
    thisgen = as.character(strs[,generation])
    total_shared = mean(apply(data.frame(thisgen),1,bitrifound))
    ACS[chain,generation-1] = total_shared
  }
}

shft = .1
ms = colMeans(ACS)
sds = apply(ACS,2,sd)
plot(ms,type='b',xlab='Generation',
     ylab='Associated chunk strength',
     ylim=c(0,0.8),main='ACS over generations')
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
regdata = cbind(which(ACS>-1,arr.ind=T),ACS[which(ACS>-1,arr.ind=T)])
lmo = lmer(regdata[,3]~regdata[,2]+(1+regdata[,2]|regdata[,1]))
coefs = data.frame(summary(lmo)$coefficients)
coefs$p = 2*(1-pnorm(abs(coefs$t.value)));coefs
t.test(ACS[,1]-ACS[,10])
mean(ACS[,1]);sd(ACS[,1])
mean(ACS[,10]);sd(ACS[,10])
```

[Here we plot the 2x2 Figure 1]

```{r}

# dev.new()
# pdf(file="Figure1.pdf",width=10,height=10)
par(mfrow = c(2, 2))
shft = .1
ms = colMeans(mindists)
sds = apply(mindists,2,sd)
plot(ms,type='b',xlab='',
     ylab='Minimal normalized Levensthein (error)',
     ylim=c(0,.45),main='Error over generations',cex.lab=1.5)
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
ms = colMeans(numzeros)
sds = apply(numzeros,2,sd)
plot(ms,type='b',xlab='',ylab='Number of exact reproductions',
     ylim=c(0,9),main='Correct strings over generations',cex.lab=1.5)
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
ms = colMeans(avglen)
sds = apply(avglen,2,sd)
plot(ms,type='b',xlab='',ylab='Average string length',
     ylim=c(0,4.5),main='Average length of strings',cex.lab=1.5)
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
ms = colMeans(ACS)
sds = apply(ACS,2,sd)
plot(ms,type='b',xlab='',
     ylab='Associative chunk strength',
     ylim=c(0,0.8),main='ACS over generations',cex.lab=1.5)
segments(1:10,ms+sds/sqrt(10),1:10,ms-sds/sqrt(10))
segments(1:10-shft,ms+sds/sqrt(10),1:10+shft,ms+sds/sqrt(10))
segments(1:10-shft,ms-sds/sqrt(10),1:10+shft,ms-sds/sqrt(10))
mtext("Generation", side=1, outer=TRUE, line=-1.5, cex=1.5)
# dev.off()

```


"Figure 2: Generations 0 (left) and 10 (right) of chain 8. These network diagrams link strings which share at least one bigram sequence"


```{r}
library(igraph)

bigramfound = function(s1,s2) {
  s1 = unlist(strsplit(s1,""))
  s2 = unlist(strsplit(s2,""))
  s1 = unique(cbind(s1[1:(length(s1)-1)],s1[2:length(s1)]))
  s2 = unique(cbind(s2[1:(length(s2)-1)],s2[2:length(s2)]))
  all = rbind(s1,s2)
  if (dim(all)[1]>dim(unique(all))[1]) {
    bigramfound = 1
  } else {
    bigramfound = 0
  }  
  return(bigramfound)
}

buildconxs = function(ix) {
    A = as.character(strings[ix[1]])
    B = as.character(strings[ix[2]])
    if (ix[1]>ix[2]) { return(bigramfound(A,B)) }
    else { return(0) }
}   

for (chain in 8) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')  
  for (generation in c(1,11)) {
    nodes = as.character(strs[,generation])
    cmb = expand.grid(i=1:15, j=1:15)
    strings <<- nodes
    conxs = matrix(apply(cmb,1,buildconxs),15,15) 
    edges = which(conxs==1,arr.ind=T)
    edges = data.frame(from=nodes[edges[,1]],to=nodes[edges[,2]])
    network=graph.data.frame(edges,directed=F,vertices=as.character(nodes))
    plot(network,layout=layout.circle,vertex.size=0,
         vertex.color=rgb(1,1,1),vertex.shape="circle",
         edge.color=rgb(0.25,0.25,0.25),
         vertex.label=as.character(nodes),
         vertex.label.color=rgb(0,0,0),vertex.label.dist=c(1:7*0+1,1:8*0-1))    
  }
}

```

[Building Figure 2]

```{r}

# cTo = c('S','R','Z','M','L','P') # see notes about translation in "prior_code"
# cFrom = c('c','b','a','d','e','f')
# dev.new()
# pdf(file="Figure2.pdf",width=10,height=10)
par(mfrow=c(1,2))
for (chain in 8) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')  
  for (generation in c(1,11)) {
    nodes = as.character(strs[,generation])
    #for (charNum in 1:6) { nodes = gsub(cFrom[charNum],cTo[charNum],nodes) }      
    cmb = expand.grid(i=1:15, j=1:15)
    strings <<- nodes
    conxs = matrix(apply(cmb,1,buildconxs),15,15) 
    edges = which(conxs==1,arr.ind=T)
    edges = data.frame(from=nodes[edges[,1]],to=nodes[edges[,2]])
    network=graph.data.frame(edges,directed=F,vertices=as.character(nodes))
    plot(network,layout=layout.circle,vertex.size=0,vertex.label='',
         vertex.color=rgb(1,1,1),vertex.shape="circle",
         edge.color=rgb(0.25,0.25,0.25),
         vertex.label.color=rgb(0,0,0),vertex.label.dist=0*cos(seq(from=0,to=2*pi,length.out=15))) 
    text(x=layout.circle(network)[,1]*1.2,y=layout.circle(network)[,2]*1.1,labels=nodes)
  }
}
# dev.off()

```

"For each set of networks, both experimental and natural-language (and their baselines), we extracted (1) string length, and (2) the proportion of other elements to which a given string is connected. The relationship between these variables is shown as blue and red lines in Figure 4"

```{r}

buildconxs_full = function(ix) {
    A = as.character(strings[ix[1]])
    B = as.character(strings[ix[2]])
    if (ix[1]!=ix[2]) { return(bigramfound(A,B)) }
    else { return(0) }
}   

shuffle_nodes_internally = function(node) {
  node = unlist(strsplit(as.character(node),""))
  node = paste(node[sample(length(node))],collapse="")
  return(node)
}

gens=c();props=c();lens=c();chains=c();shuff=c()
for (chain in 1:8) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')  
  for (generation in c(1:11)) {
    nodes = as.character(strs[,generation])
    cmb = expand.grid(i=1:15, j=1:15)
    strings <<- nodes
    conxs = matrix(apply(cmb,1,buildconxs_full),15,15) # include some connection matrix
    lens = c(lens,nchar(as.character(nodes))) # length of the strings in this generation
    props = c(props,colMeans(conxs)) # proportion to other nodes connected
    gens = c(gens,1:15*0+generation)
    chains = c(chains,1:15*0+chain)
    shuff = c(shuff,1:15*0+0)
    # time to do the shuffled versions
    nodes = apply(data.frame(nodes),1,shuffle_nodes_internally)
    cmb = expand.grid(i=1:15, j=1:15)
    strings <<- nodes
    conxs = matrix(apply(cmb,1,buildconxs_full),15,15) # include some connection matrix
    lens = c(lens,nchar(as.character(nodes))) # length of the strings in this generation
    props = c(props,colMeans(conxs)) # proportion to other nodes connected
    gens = c(gens,1:15*0+generation)
    chains = c(chains,1:15*0+chain)    
    shuff = c(shuff,1:15*0+1)
  }
}

net_data = data.frame(lens,props,gens,chains,shuff)

```

"In the first panel, Generation 0, the shuffled strings (red) are in fact significantly greater in their overall connectivity, t = 7.5, p < .0001. This gradually changes, and by the final three generations (8, 9, 10) the original data are more greatly connected as a function of string length, t's > 2.5,  p's < .005."

```{r}
subd = net_data[net_data$gens==1,]
subd$lensC = subd$lens-mean(subd$lens)
summary(lm(props~shuff*lensC,data=subd))

subd = net_data[net_data$gens==9,]
subd$lensC = subd$lens-mean(subd$lens)
summary(lm(props~shuff*lensC,data=subd))

subd = net_data[net_data$gens==10,]
subd$lensC = subd$lens-mean(subd$lens)
summary(lm(props~shuff*lensC,data=subd))

subd = net_data[net_data$gens==11,]
subd$lensC = subd$lens-mean(subd$lens)
summary(lm(props~shuff*lensC,data=subd))
```

"We did this same analysis across our generations of the experiment, shown in Figure 3."

```{r}
subd = net_data[net_data$gens==1,]
plot(props~lens,data=subd[subd$shuff==0,],
     col=rgb(.55,.55,.9),pch=1,cex=1.25,ylim=c(0,1),xlim=c(3,6))
points(props~lens,data=subd[subd$shuff==1,],col=rgb(.9,.55,.55),pch=15,cex=.5)
lmo = lm(props~lens,data=subd[subd$shuff==0,]);
  abline(lmo,col='blue',lwd=2,xlim=c(3,6),ylim=c(0,1),xpd=F)
lmo = lm(props~lens,data=subd[subd$shuff==1,]);
  abline(lmo,col='red',lwd=2,xpd=F)

subd = net_data[net_data$gens==11,]
plot(props~lens,data=subd[subd$shuff==0,],
     col=rgb(.55,.55,.9),pch=1,cex=1.25,ylim=c(0,1),xlim=c(3,6))
points(props~lens,data=subd[subd$shuff==1,],col=rgb(.9,.55,.55),pch=15,cex=.5)
lmo = lm(props~lens,data=subd[subd$shuff==0,]);
  abline(lmo,col='blue',lwd=2,xlim=c(3,6),ylim=c(0,1),xpd=F)
lmo = lm(props~lens,data=subd[subd$shuff==1,]);
  abline(lmo,col='red',lwd=2,xpd=F)
```

"For the natural-language (CHILDES) network, the original data (unshuffled) have overall greater connectivity than the shuffled data by (on average) 10%, t = 47.6, p < .0001, and the interaction in Fig. 4 (bottom right) is significant, t = 20.7,  p < .0001. Importantly, these effects are still present when just focusing on strings of length 3 and 4 alone: It is not driven by the longer string sequences alone (p's < .0001)."

```{r}
pos = read.table('data/human_POS_for_R.xls')
colnames(pos) = list('shuff','lens','props')
subd = pos # just to match the code above
plot(props~lens,data=subd[subd$shuff==0,],
     col=rgb(.55,.55,.9),pch=1,cex=1.25,ylim=c(0,1),xlim=c(3,6))
points(props~lens,data=subd[subd$shuff==1,],col=rgb(.9,.55,.55),pch=15,cex=.5)
lmo = lm(props~lens,data=subd[subd$shuff==0,]);
  abline(lmo,col='blue',lwd=2,xlim=c(3,6),ylim=c(0,1),xpd=F)
lmo = lm(props~lens,data=subd[subd$shuff==1,]);
  abline(lmo,col='red',lwd=2,xpd=F)

subd$lensC = subd$lens-mean(subd$lens)
summary(lm(props~shuff*lensC,data=subd))
summary(lm(props~shuff*lensC,data=subd[subd$lens<=4,]))
```

"We can now compare the human part-of-speech data to the experimental data directly, because they can be compared on the same scale (proportion of connectivity). In the final three generations (8,9,10) the human data does statistically differ from the experimental data in extent of connectivity. In particular, the experimental data are more connected, by about 7% (p < .0001). This is likely due to the fact that the POS CHILDES data involve more categories (parts of speech), and thus more bigram types, and lower probability of drawing edges between sequences. Importantly, the interaction term in this analysis is not significant (p = .72): There is no difference in slope for how proportion connectivity scales with string length. However, the human data do differ from the first three experimental generations considered together (0,1,2). The human data show considerably more connectivity, and the interaction term is significant (p < .0001), suggesting that natural-language connectivity scales more robustly with length than the first few generations of the experiment, but more similarly to the final three generations."

```{r}
subd1 = net_data[net_data$gens>=9,];subd1$human=0
subd1 = subset(subd1,subd1$shuff==0,c('human','lens','props'))
subd2 = pos;subd2$human=1
subd2 = subset(subd2,subd2$shuff==0,c('human','lens','props'))
compare_data = rbind(subd1,subd2)
compare_data$lensC = compare_data$lens-mean(compare_data$lens)
summary(lm(props~human*lensC,data=compare_data))

subd1 = net_data[net_data$gens %in% c(2,3,4),];subd1$human=0
subd1 = subset(subd1,subd1$shuff==0,c('human','lens','props'))
subd2 = pos;subd2$human=1
subd2 = subset(subd2,subd2$shuff==0,c('human','lens','props'))
compare_data = rbind(subd1,subd2)
compare_data$lensC = compare_data$lens-mean(compare_data$lens)
summary(lm(props~human*lensC,data=compare_data))
```

[Building Figure 3]

```{r}

# dev.new()
# pdf(file="Figure3.pdf",width=10,height=10)
par(mfrow=c(3,4),oma = c(5,4,0,0) + 0.7,mar = c(0,0,1,1) + 0.7)
for (generation in c(1:11)) {
  subd = net_data[net_data$gens==generation,]
  xset = "n";yset = "n"
  if (generation<=8) { xset="n" }
  if (generation %in% c(1,5,9)) { yset="s" }  
  plot(props~lens,data=subd[subd$shuff==0,],
       col=rgb(.55,.55,.9),pch=1,cex=1.25,ylim=c(0,1),xlim=c(3,6),
       main=paste('Generation',generation-1),xaxt=xset,yaxt=yset)
  if (generation>8) { axis(1,at=c(3:6)) }
  points(props~lens,data=subd[subd$shuff==1,],col=rgb(.9,.55,.55),pch=15,cex=.5)
  lmo = lm(props~lens,data=subd[subd$shuff==0,]);
    abline(lmo,col='blue',lwd=2,xlim=c(3,6),ylim=c(0,1),xpd=F)
  lmo = lm(props~lens,data=subd[subd$shuff==1,]);
    abline(lmo,col='red',lwd=2,xpd=F)
}

title(ylab='Proportion strings connected',xlab='String size',outer=T,cex.lab=1.5)

pos = read.table('data/human_POS_for_R.xls')
colnames(pos) = list('shuff','lens','props')
subd = pos # just to match the code above
plot(props~lens,data=subd[subd$shuff==0,],
     col=rgb(.55,.55,.9),pch=1,cex=1.25,ylim=c(0,1),xlim=c(3,6),
     main="Human CHILDES POS",yaxt="n",xaxt="n")
axis(1,at=c(3:6))
points(props~lens,data=subd[subd$shuff==1,],col=rgb(.9,.55,.55),pch=15,cex=.5)
lmo = lm(props~lens,data=subd[subd$shuff==0,]);
  abline(lmo,col='blue',lwd=2,xlim=c(3,6),ylim=c(0,1),xpd=F)
lmo = lm(props~lens,data=subd[subd$shuff==1,]);
  abline(lmo,col='red',lwd=2,xpd=F)

# dev.off()

```

Checking how many elements CHILDES POS has, hence higher connectivity for experiments vs. human data.

```{r}
pos_list = read.table('data/Mega-Eng-POS-CHILDES.txt',sep='\t')
items = unlist(strsplit(paste(pos_list$V2,sep=' '),' '))
unique(items)
length(unique(items))-1
```

[Fig. 4 PDF is generated in the supplementary-analysis folder.]

[Fig. 5 showing contrast in generation 1 and 10 in one chain]

```{r}

# cTo = c('S','W','L','C','X','K') # for chain 3 (had to determine by hand)
# cFrom = c('c','b','a','d','e','f') # see notes about translation in prior_code
#cTo = c('c','b','a','d','e','f')
# dev.off();dev.new()
# pdf(file="Figure5.pdf",width=70,height=70)
par(mfrow=c(1,2),oma = c(5,4,0,0) + 1,mar = c(0,0,1,1) + 1)
for (chain in 3) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')  
  for (generation in c(1,11)) {
    nodes = as.character(strs[,generation])
    cmb = expand.grid(i=1:15, j=1:15)
    # for (charNum in 1:6) { nodes = gsub(cFrom[charNum],cTo[charNum],nodes) }   
    strings <<- nodes
    conxs = matrix(apply(cmb,1,buildconxs),15,15) 
    dists = adist(strings,strings)
    edges = which(conxs==1,arr.ind=T)
    weights = 10*1/dists[conxs==1]
    edges = data.frame(from=nodes[edges[,1]],to=nodes[edges[,2]])
    network=graph.data.frame(edges,directed=F,vertices=as.character(nodes))
    E(network)$weights = weights
    plot(network,layout=layout.kamada.kawai,vertex.size=66,vertex.label=nodes,
         vertex.color=rgb(1,1,1),vertex.shape="rectangle",edge.width=weights,
         edge.color=rgb(0.25,0.25,0.25),
         vertex.label.color=rgb(0,0,0),vertex.label.cex=.7) 
    #text(x=layout.circle(network)[,1]*1.2,y=layout.circle(network)[,2]*1.1,labels=nodes)
  }
}
# dev.off()

```

[All networks using Kamada-Kawai force-directed algorithm for visualization; see Fig. S1 PDF]

```{r}

dev.off();dev.new()
pdf(file="FigureS1.pdf",width=70,height=70)
par(mfrow=c(8,11))
for (chain in 1:8) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')  
  for (generation in c(1:11)) {
    nodes = as.character(strs[,generation])
    cmb = expand.grid(i=1:15, j=1:15)
    strings <<- nodes
    conxs = matrix(apply(cmb,1,buildconxs),15,15) 
    dists = adist(strings,strings)
    edges = which(conxs==1,arr.ind=T)
    weights = 10*1/dists[conxs==1]
    edges = data.frame(from=nodes[edges[,1]],to=nodes[edges[,2]])
    network=graph.data.frame(edges,directed=F,vertices=as.character(nodes))
    E(network)$weights = weights
    plot(network,layout=layout.kamada.kawai,vertex.size=50,vertex.label=nodes,
         vertex.color=rgb(1,1,1),vertex.shape="rectangle",edge.width=weights,
         edge.color=rgb(0.25,0.25,0.25),
         vertex.label.color=rgb(0,0,0)) 
    #text(x=layout.circle(network)[,1]*1.2,y=layout.circle(network)[,2]*1.1,labels=nodes)
  }
}
dev.off()

```

[All networks for demonstration; using bigrams; see FigureS2.pdf]

```{r}

dev.off();dev.new()
pdf(file="FigureS2.pdf",width=70,height=70)
par(mfrow=c(8,11))
for (chain in 1:8) {
  # load this chain file
  strs = read.table(paste('data/chain',chain,'.chain',sep=''),sep='\t')  
  for (generation in c(1:11)) {
    nodes = as.character(strs[,generation])
    cmb = expand.grid(i=1:15, j=1:15)
    strings <<- nodes    
    conxs = matrix(apply(cmb,1,buildconxs),15,15) 
    edges = which(conxs==1,arr.ind=T)
    edges = data.frame(from=nodes[edges[,1]],to=nodes[edges[,2]])
    network=graph.data.frame(edges,directed=F,vertices=as.character(nodes))
    plot(network,layout=layout.circle,vertex.size=0,vertex.label='',
         vertex.color=rgb(1,1,1),vertex.shape="circle",
         edge.color=rgb(0.25,0.25,0.25),
         vertex.label.color=rgb(0,0,0),vertex.label.dist=0*cos(seq(from=0,to=2*pi,length.out=15))) 
    text(x=layout.circle(network)[,1]*1.2,y=layout.circle(network)[,2]*1.1,labels=nodes)
  }
}
dev.off()

```