Reshaping data

Preparing and reshaping data is the ever continuing task of a data analyst. Luckily we have many tools for it. The default tool in R would be reshape(), although this is so user friendly that a reshape package has been added too. I try to use reshape() (the function) because I feel it is a good tool, though with a somewhat cryptical manual. The latter may be because it is written in terms of longitudal data, whereas my experience is  conversion of data from easy to enter in Excel to suitable for analysis in R.
To exercise myself a bit more I have taken all examples from the SAS transpose procedure and implemented them in R.

Examples 1 to 3

These examples are so simple, the best tool is the t() function.

score <- read.table(textConnection('
Student StudentID  Section  Test1 Test2 Final
Capalleti 0545 1  94 91 87
Dubose    1252 2  51 65 91
Engles    1167 1  95 97 97
Grant     1230 2  63 75 80
Krupski   2527 2  80 76 71
Lundsford 4860 1  92 40 86
McBane    0674 1  75 78 72'),header=TRUE,
  colClasses=rep(c('character','numeric'),each=3))
# id, only numerical data
# this case - t(score[,-1:-3])
# general
t(score[,sapply(score,class)=='numeric'])

      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
Test1   94   51   95   63   80   92   75
Test2   91   65   97   75   76   40   78
Final   87   91   97   80   71   86   72

#example 2: studentid
score2 <- score
rownames(score2) <- paste('sn',score2$StudentID,sep='')
t(score2[,-1:-3])

      sn0545 sn1252 sn1167 sn1230 sn2527 sn4860 sn0674
Test1     94     51     95     63     80     92     75
Test2     91     65     97     75     76     40     78
Final     87     91     97     80     71     86     72

#example 3 student names

rownames(score2) <- score2$Student
t(score2[,-1:-3])

      Capalleti Dubose Engles Grant Krupski Lundsford McBane
Test1        94     51     95    63      80        92     75
Test2        91     65     97    75      76        40     78
Final        87     91     97    80      71        86     72

Example 4 by groups

The first 'real' example. Some columns are used to transpose on. In addition the data is ragged, not everywhere there is data. Two unexpected problems appeared here. Location contains a space, which I solved by separating this part and adding it later. Date is in a 7 digit format, but not my locale and a non default layout. Since I need to sort in order to get the data ordered like SAS does, I converted it to a proper date variable.
The transpose itself is not so difficult. I chose to extract the variables to transpose via grep(), so I could reuse this part in the times variable. To make tracking between call and result easy, I used lower and upper case in 'length'. Length is transposed, so weight is dropped from the data. 

# example 4 by groups

fishdata1 <- readLines(textConnection(
'Cole Pond   2JUN95 31 .25 32 .3  32 .25 33 .3
Cole Pond   3JUL95 33 .32 34 .41 37 .48 32 .28
Cole Pond   4AUG95 29 .23 30 .25 34 .47 32 .3
Eagle Lake  2JUN95 32 .35 32 .25 33 .30
Eagle Lake  3JUL95 30 .20 36 .45
Eagle Lake  4AUG95 33 .30 33 .28 34 .42'))
fishdata <- read.table(text=c('Date Length1 Weight1 Length2 Weight2 Length3 Weight3 Length4 Weight4',
           substring(fishdata1,11)),
       fill=TRUE,
       header=TRUE)
fishdata$Location <- gsub(' $','',substring(fishdata1,1,10))
Sys.setlocale(category = "LC_TIME", locale = "C")
fishdata$Date <- as.Date(fishdata$Date,'%d%b%y')
lengthfishdata <- reshape(fishdata,
    varying=list(grep('Length',names(fishdata),value=TRUE)),
    direction='long',
    idvar=c('Date','Location'),
    drop=grep('Weight',names(fishdata),value=TRUE),
    v.names='LENGTH',
    timevar='NAME',
    times=tolower(grep('Length',names(fishdata),value=TRUE)),
  )
rownames(lengthfishdata)  <- 1:nrow(lengthfishdata)
lengthfishdata[order(lengthfishdata$Location,
        lengthfishdata$Date,
        lengthfishdata$NAME),]

         Date   Location    NAME LENGTH
1  1995-06-02  Cole Pond length1     31
7  1995-06-02  Cole Pond length2     32
13 1995-06-02  Cole Pond length3     32
19 1995-06-02  Cole Pond length4     33
2  1995-07-03  Cole Pond length1     33
8  1995-07-03  Cole Pond length2     34
14 1995-07-03  Cole Pond length3     37
20 1995-07-03  Cole Pond length4     32
3  1995-08-04  Cole Pond length1     29
9  1995-08-04  Cole Pond length2     30
15 1995-08-04  Cole Pond length3     34
21 1995-08-04  Cole Pond length4     32
4  1995-06-02 Eagle Lake length1     32
10 1995-06-02 Eagle Lake length2     32
16 1995-06-02 Eagle Lake length3     33
22 1995-06-02 Eagle Lake length4     NA
5  1995-07-03 Eagle Lake length1     30
11 1995-07-03 Eagle Lake length2     36
17 1995-07-03 Eagle Lake length3     NA
23 1995-07-03 Eagle Lake length4     NA
6  1995-08-04 Eagle Lake length1     33
12 1995-08-04 Eagle Lake length2     33
18 1995-08-04 Eagle Lake length3     34
24 1995-08-04 Eagle Lake length4     NA
In practice my call would be different, I would keep both height and weight, and stick the number in a more logical variable than NAME. I think SAS can only achieve this by two transpose calls and a subsequent merge, though I may be mistaken.
reshape(fishdata,
    varying=list(grep('Length',names(fishdata),value=TRUE),
        grep('Weight',names(fishdata),value=TRUE)),
    direction='long',
    idvar=c('Date','Location'),
    v.names=c('Length','Weight'),
    timevar='Fish number',
    new.row.names=1:24
)

         Date   Location Fish number Length Weight
1  1995-06-02  Cole Pond           1     31   0.25
2  1995-07-03  Cole Pond           1     33   0.32
3  1995-08-04  Cole Pond           1     29   0.23
4  1995-06-02 Eagle Lake           1     32   0.35
5  1995-07-03 Eagle Lake           1     30   0.20
6  1995-08-04 Eagle Lake           1     33   0.30
7  1995-06-02  Cole Pond           2     32   0.30
8  1995-07-03  Cole Pond           2     34   0.41
9  1995-08-04  Cole Pond           2     30   0.25
10 1995-06-02 Eagle Lake           2     32   0.25
11 1995-07-03 Eagle Lake           2     36   0.45
12 1995-08-04 Eagle Lake           2     33   0.28
13 1995-06-02  Cole Pond           3     32   0.25
14 1995-07-03  Cole Pond           3     37   0.48
15 1995-08-04  Cole Pond           3     34   0.47
16 1995-06-02 Eagle Lake           3     33   0.30
17 1995-07-03 Eagle Lake           3     NA     NA
18 1995-08-04 Eagle Lake           3     34   0.42
19 1995-06-02  Cole Pond           4     33   0.30
20 1995-07-03  Cole Pond           4     32   0.28
21 1995-08-04  Cole Pond           4     32   0.30
22 1995-06-02 Eagle Lake           4     NA     NA
23 1995-07-03 Eagle Lake           4     NA     NA
24 1995-08-04 Eagle Lake           4     NA     NA

Example 5

Example 5 is named: Naming Transposed Variables When the ID Variable Has Duplicate Values. I am not sure what the naming part is, it seems that only the closing call price is needed, which boils down to taking only the last observation in a category. The data here comes in a fixed format array, I have used (first time) read.fwf() and some calls to strip the spaces from the factors.
To get the transpose I borrowed the idea of a last function from SAS, which is basically a function which indicates that the current record in a variable is different from the next record. Which is very important in SAS because it is fundamentally row (record) organized, whereas R is column (variable) organized. Proper R would just processing dependent on Time='closing' but what is used is functionally closer to the SAS call.

stocks <- read.fwf(textConnection(
'Horizon Kites jun11 opening 29
Horizon Kites jun11 noon    27
Horizon Kites jun11 closing 27
Horizon Kites jun12 opening 27
Horizon Kites jun12 noon    28
Horizon Kites jun12 closing 30
SkyHi Kites   jun11 opening 43
SkyHi Kites   jun11 noon    43
SkyHi Kites   jun11 closing 44
SkyHi Kites   jun12 opening 44
SkyHi Kites   jun12 noon    45
SkyHi Kites   jun12 closing 45'
),col.names=c('Company','Date','Time','Price'),
widths=c(14,5,8,3))
levels(stocks$Company) <- gsub('(^ +)|( +$)','',levels(stocks$Company))
levels(stocks$Time) <- gsub('(^ +)|( +$)','',levels(stocks$Time))
# only last observation
islast <- function(x)  c(x[-length(x)]!=x[-1],TRUE)
reshape(stocks[islast(stocks$Date),],direction='wide',
    timevar='Date',idvar=c('Company'),drop='Time',times=Time)

        Company Price.jun11 Price.jun12
3 Horizon Kites          27          30
9   SkyHi Kites          44          45

Example 6 

Transposing for statistical analysis. We have 5 subjects, who did 3 programs and 7 strength assessments. Why subject is not in the original data baffles me. In SAS this is not run through proc transpose, but rather a data step. But that doesn't stop me from using reshape(). The subject variable is added. In the next step SAS rebuilds to get one file with both formats, seems silly in R context, I am just transforming back, now using subject. 

weights <- read.table(textConnection(
'Program s1 s2 s3 s4 s5 s6 s7
CONT  85 85 86 85 87 86 87
CONT  80 79 79 78 78 79 78
CONT  78 77 77 77 76 76 77
CONT  84 84 85 84 83 84 85
CONT  80 81 80 80 79 79 80
RI    79 79 79 80 80 78 80
RI    83 83 85 85 86 87 87
RI    81 83 82 82 83 83 82
RI    81 81 81 82 82 83 81
RI    80 81 82 82 82 84 86
WI    84 85 84 83 83 83 84
WI    74 75 75 76 75 76 76
WI    83 84 82 81 83 83 82
WI    86 87 87 87 87 87 86
WI    82 83 84 85 84 85 86
'),header=TRUE)
weights1 <- reshape(weights,direction='long',timevar='time',idvar='row',
    varying=list(paste('s',1:7,sep='')),v.names='Strength')
weights1$subject <- 1+ ((weights1$row-1) %% 5)
weights1 <- weights1[order(weights1$Program,weights1$subject,weights1$time),]
(Weights1 <- weights1[,c(1,5,2,3)])[1:15,]

    Program subject time Strength
1.1    CONT       1    1       85
1.2    CONT       1    2       85
1.3    CONT       1    3       86
1.4    CONT       1    4       85
1.5    CONT       1    5       87
1.6    CONT       1    6       86
1.7    CONT       1    7       87
2.1    CONT       2    1       80
2.2    CONT       2    2       79
2.3    CONT       2    3       79
2.4    CONT       2    4       78
2.5    CONT       2    5       78
2.6    CONT       2    6       79
2.7    CONT       2    7       78
3.1    CONT       3    1       78

reshape(Weights1,direction='wide',timevar='time',

    idvar=c('Program','subject'))

     Program subject Strength.1 Strength.2 Strength.3 Strength.4 Strength.5
1.1     CONT       1         85         85         86         85         87
2.1     CONT       2         80         79         79         78         78
3.1     CONT       3         78         77         77         77         76
4.1     CONT       4         84         84         85         84         83
5.1     CONT       5         80         81         80         80         79
6.1       RI       1         79         79         79         80         80
7.1       RI       2         83         83         85         85         86
8.1       RI       3         81         83         82         82         83
9.1       RI       4         81         81         81         82         82
10.1      RI       5         80         81         82         82         82
11.1      WI       1         84         85         84         83         83
12.1      WI       2         74         75         75         76         75
13.1      WI       3         83         84         82         81         83
14.1      WI       4         86         87         87         87         87
15.1      WI       5         82         83         84         85         84
     Strength.6 Strength.7
1.1          86         87
2.1          79         78
3.1          76         77
4.1          84         85
5.1          79         80
6.1          78         80
7.1          87         87
8.1          83         82
9.1          83         81
10.1         84         86
11.1         83         84
12.1         76         76
13.1         83         82
14.1         87         86
15.1         85         86

Simulation shows gain of clmm over ANOVA is small

After last post's setting up for a simulation, it is now time to look how the models compare. To my disappointment with my simple simulations of assessors behavior the gain is minimal. Unfortunately, the simulation took much more time than I expected, so I will not expand it.

Background

I have been looking at the ordered logistic model in a number of postings. The reason is that in general I find people use ANOVA for analyzing data on a nine point scale, whereas you would think an ordered logistic model works better. Two posts (first and second) showed the current methods in R, and and JAGS are easy to use and with some tweaking provide suitable output to present. Subsequently I made the simulated data generator. Now it is time to make the final comparison.  

Simulations

The core of the simulator is explained elsewhere so I won't explain here again. I did however notice a small error, so the corrected code is given here. Some new parts are added, wrappers around the data generator and the analysis. And to my big disappointment I could not even build that as desired. The call anova(Res.clmm2,Res.clmm) with subsequent extraction has been replaced by the ugly pchisq(2*(Res.clmm$logLik-Res.clmm2$logLik),2,lower.tail=FALSE ). 2 represents the degrees of freedom for products in my simulations. Somehow that call to ANOVA did not run within a function, after trying too many variations I choose the short cut.

library(ordinal) 
library(ggplot2)

num2scale <- function(x,scale) {
  cut(x,c(-Inf,scale,Inf),1:(length(scale)+1))
}
pop.limits2ind.limits <- function(scale,sd=.5) {
  newscale <- scale+rnorm(length(scale),sd=sd)
  newscale[order(newscale)]
}

obs.score <- function(obs,pop.score,pop.limits,
    sensitivity.sd=.1, precision.sd=1,
    additive.sd=2,center.scale=5,
    labels=LETTERS[1:length(pop.score)]) {
  # individual sensitivity (multiplicative)
  obs.errorfreeintensity <- center.scale + 
      (pop.score-center.scale)*rlnorm(1,sd=sensitivity.sd)
  #individual (additive) 
  obs.errorfreeintensity <- obs.errorfreeintensity +
      rnorm(1,mean=0,sd=additive.sd)
  # individual observation error 
  obs.intensity <- obs.errorfreeintensity+
      rnorm(length(pop.score),sd=precision.sd)
  
  # individual cut offs between categories  
  obs.limits <- pop.limits2ind.limits(pop.limits)
  obs.score <- num2scale(obs.intensity,obs.limits)
  data.frame(obs=obs,
      score = obs.score,
      product=labels
  )
}

panel.score <- function(n=100,pop.score,pop.limits,
    sensitivity.sd=.1,precision.sd=1,
    additive.sd=2,center.scale=5,
    labels=LETTERS[1:length(pop.score)]) {
  la <- lapply(1:n,function(x) {
        obs.score(x,pop.score,pop.limits,sensitivity.sd=sensitivity.sd,
            precision.sd=precision.sd,additive.sd=additive.sd,labels=labels)
      })
  dc <- do.call(rbind,la)
  dc$obs <- factor(dc$obs)
  dc
}

overallP <- function(scores) {
  Res.aov <- aov( numresponse ~ obs +  product , data=scores)
  paov <- anova(Res.aov)['product','Pr(>F)']
  Res.clmm <- clmm(score  ~ product + (1|obs),data=scores)
  Res.clmm2 <- clmm(score ~ 1 + (1|obs),data=scores)
  c(ANOVA=paov,
    clmm = pchisq(2*(Res.clmm$logLik-Res.clmm2$logLik),2,lower.tail=FALSE ) #.1687
   )
}

onesim <- function(prodmeans,pop.limits,center.scale,additive.sd=2) {
  scores <- panel.score(40,prodmeans,pop.limits,
      center.scale=center.scale,additive.sd=additive.sd)
  scores$numresponse <- as.numeric(levels(scores$score))[scores$score]
  overallP(scores)
}

Simulation I, 5 categories

The first simulation is with 5 categories. This represents the Just About Right (JAR) scale. The final plot shows the difference between ANOVA and clmm is minimal.

pop.limits <- c(1,2.5,4.5,6)

nsim <- 250
sim5cat1 <- lapply(seq(0,.6,.05),function(dif) {
    prodmeans=rep(3.5,3)+c(-dif,0,dif)
    sa <- sapply(1:nsim,function(x) onesim(prodmeans,pop.limits))
    data.frame(dif=dif,nreject=as.numeric(rowSums(sa<0.05)),
        method=rownames(sa))
  }
)

sim5cat1tb <- do.call(rbind,sim5cat1)
ggplot(sim5cat1tb, aes(dif,nreject/nsim ,colour=method)) +
     geom_line() + xlab('Difference between products') + 
     ylab('Proportion significant (at 5% Test)') +
     theme(legend.position = "bottom") + ylim(c(0,1)) +
     guides(colour = guide_legend('Analysis Method'))

Simulation 2, 9 categories

This simulation represents the intensity and liking scales. Again the difference between ANOVA and clmm are minimal.

pop.limits <- c(1,3:8,10)
prodmeans <- c(7,7,7)
scores <- panel.score(40,prodmeans,pop.limits,additive.sd=1)
scores$numresponse <- as.numeric(levels(scores$score))[scores$score]
nsim <- 250
sim9cat <- lapply(seq(0.,.6,.05),function(dif) {
      prodmeans=c(7,7,7)+c(-dif,0,dif)
      sa <- sapply(1:nsim,function(x) onesim(prodmeans,pop.limits,
                additive.sd=1))
      data.frame(dif=dif,nsign=as.numeric(rowSums(sa>0.05)),
          method=rownames(sa))
    }
)
sim9cattb <- do.call(rbind,sim9cat)
ggplot(sim9cattb, aes(dif,(nsim-nsign)/nsim ,colour=method)) +
    geom_line() + xlab('Difference between products') + 
    ylab('Proportion significant (at 5% Test)') +
    theme(legend.position = "bottom") + ylim(c(0,1)) +
    guides(colour = guide_legend('Analysis Method'))

Conclusion

The differences between clmm and ANOVA seem to be minimal. I had not expected large differences, but especially at 5 categories my expectation were to find real differences as the continuous data is more violated. Obviously, a load more simulations would be needed to draw final conclusions. This is beyond the scope of my blog. 

To conclude, in theory clmm is much more suitable than ANOVA for ordinal data. There are no reasons in terms of presentation to prefer ANOVA over clmm. But in practice the difference may be minimal.

Ordinal data, models with observers

I recently made three posts regarding analysis of ordinal data. A post looking at all methods I could find in R, a post with an additional method and a post using JAGS. Common in all three was using the cheese data, a data set where only product data was available, observers were not there. The real world is different, there are observers and they are quite different. So, in this post I have added observers and look at two models again; Simple ANOVA and, complex, an ordinal model with random observers.

Data

To spare myself looking for data, and keeping in mind I may want to run some simulations, I have created a data generator. In this generator I entered some of the things I could imagine regarding the observers. First of all they all have different sensitivities. In other words, what is very salt for some, is not at all salt for others. This is expressed in an additive and a multiplicative effect. Second, they have observation error. The same product does not always give the same response. Third is the categories, observers use them slightly different. I am absolutely aware these are much more properties than I might estimate. After all, a typical consumer takes two or three products in a test, while I have many more parameters. Finally, the outer limits (categories 1 and 9) are placed a bit outward. This is because there is an aversion to using these categories.
num2scale <- function(x,scale) {
  cut(x,c(-Inf,scale,Inf),1:(length(scale)+1))
}
pop.limits2ind.limits <- function(scale,sd=.5) {
  newscale <- scale+rnorm(length(scale),sd=sd)
  newscale[order(newscale)]
}

obs.score <- function(obs,pop.score,pop.limits,scale.sd=0.5,
    sensitivity.sd=.1, precision.sd=1,
    additive.sd=2,center.scale=5,
    labels=LETTERS[1:length(pop.score)]) {
  # individual sensitivity (multiplicative)
  obs.errorfreeintensity <- center.scale + 
      (pop.score-center.scale)*rlnorm(1,sd=sensitivity.sd)
  #individual (additive) 
  obs.errorfreeintensity <- obs.errorfreeintensity +
      rnorm(1,mean=0,sd=additive.sd)
  # individual observation error 
  obs.intensity <- obs.errorfreeintensity+
      rnorm(length(pop.score),sd=precision.sd)
  # individual cut offs between categories  
  obs.limits <- pop.limits2ind.limits(pop.limits)
  obs.score <- num2scale(obs.intensity,obs.limits)
  data.frame(obs=obs,
      score = obs.score,
      product=labels
  )
}

panel.score <- function(n=100,pop.score,pop.limits,scale.sd=0.5,
    sensitivity.sd=.1,precision.sd=1,
    additive.sd=2,center.scale=5,
    labels=LETTERS[1:length(pop.score)]) {
  la <- lapply(1:n,function(x) {
        obs.score(x,pop.score,pop.limits,scale.sd,sensitivity.sd,
            precision.sd,labels=labels)
      })
  dc <- do.call(rbind,la)
  dc$obs <- factor(dc$obs)
  dc
}

pop.limits <- c(1,3:8,10)
prodmeans <- c(7,7,8)
scores <- panel.score(40,prodmeans,pop.limits)
scores$numresponse <- as.numeric(levels(scores$score))[scores$score]

Analysis - ANOVA

The first analysis method is ANOVA. Plain ANOVA. Not because I dislike variance components or such, but because this is what is done most commonly. On top of that, I cannot imagine somebody taking these data, pulling it through a mixed model, while forgetting it is not from a continuous scale.
library(multcomp)
Res.aov <- aov( numresponse ~ obs + product , data=scores)
anova(Res.aov)

Analysis of Variance Table

Response: numresponse
          Df  Sum Sq Mean Sq F value    Pr(>F)    
obs       39 239.300  6.1359  6.3686 2.321e-12 ***
product    2   9.517  4.7583  4.9388   0.00956 ** 
Residuals 78  75.150  0.9635                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

mt <- model.tables(Res.aov,type='means',cterms='product') 

data.frame(dimnames(mt$tables$product)[1],
    Response=as.numeric(mt$tables$product),
    LSD=cld(glht(Res.aov,linfct=mcp(product='Tukey'))
    )$mcletters$monospacedLetters
)    
  product Response LSD
A       A    6.725  a 
B       B    6.850  a 
C       C    7.375   b

Analysis - cumulative link mixed model (clmm)

The alternative, is the other extreme. Do it as correct as available, using both a cumulative link and making observers random. This model is standard available in the ordinal package. This means two intermediate models are not shown. Both of these models lack the simplicity of ANOVA while being less correct than clmm, hence inferior to both clmm and ANOVA.
library(ordinal) 
Res.clmm <- clmm(score  ~ product + (1|obs),data=scores)
summary(Res.clmm)

Cumulative Link Mixed Model fitted with the Laplace approximation

formula: score ~ product + (1 | obs)
data:    scores

 link  threshold nobs logLik  AIC    niter    max.grad cond.H
 logit flexible  120  -179.57 379.14 37(5046) 8.77e-06 Inf   

Random effects:
      Var Std.Dev
obs 8.425   2.903
Number of groups:  obs 40 

Coefficients:
         Estimate Std. Error z value Pr(>|z|)   
productB   0.1309     0.4225   0.310  0.75680   
productC   1.4640     0.4718   3.103  0.00192 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Threshold coefficients:
    Estimate Std. Error z value
2|3  -6.0078         NA      NA
3|4  -5.6340         NA      NA
4|5  -4.1059     0.4528  -9.069
5|6  -3.0521     0.4769  -6.400
6|7  -1.0248     0.4970  -2.062
7|8   0.5499     0.5303   1.037
8|9   4.0583     0.7429   5.463

Res.clmm2 <- clmm(score ~ 1 + (1|obs),data=scores)

anova(Res.clmm2,Res.clmm)

Likelihood ratio tests of cumulative link models:

          formula:                    link: threshold:
Res.clmm2 score ~ 1 + (1 | obs)       logit flexible  
Res.clmm  score ~ product + (1 | obs) logit flexible  

          no.par    AIC  logLik LR.stat df Pr(>Chisq)   
Res.clmm2      8 387.41 -185.70                         
Res.clmm      10 379.14 -179.57  12.266  2    0.00217 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

co <- coef(Res.clmm)[c('7|8','productB','productC')]

vc <- vcov(Res.clmm)[c('7|8','productB','productC'),
    c('7|8','productB','productC')]
names(co) <- levels(scores$product)
sd <- sqrt(c(vc[1,1],diag(vc)[-1]+vc[1,1]-2*vc[1,-1]))
data.frame(
    `top 2 box`=arm::invlogit(c(-co[1],-co[1]+co[-1])),
    `2.5% limit`=arm::invlogit(c(-co[1],-co[1]+co[-1])+qnorm(.025)*sd),
    `97.5% limit`=arm::invlogit(c(-co[1],-co[1]+co[-1])+qnorm(.975)*sd),
    check.names=FALSE
)

  top 2 box 2.5% limit 97.5% limit
A 0.3658831  0.1694873   0.6199713
B 0.3967401  0.1753347   0.6704337
C 0.7138311  0.4498042   0.8838691

Continuing Sync

I am continuing in Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life by Steven Strogatz. To get a feeling on it, I was building a group of things which have only a minute influence on each other are able to synchronize their behavior. I got some more explanations how things fit together and adapted my calculations. Two things are changed.
You can actually use one of the items as a kind of clock, only register status when it goes to zero. This makes it possible to plot behavior over much longer times. Then, I found my function of last week a bit primitive, by going step-wise through time. I was thinking at first to reduce the step-size  why not do the more exact thing and follow a function? So, I scaled the logit function a bit to provide behavior of an item and build some functions to determine when an item goes over the threshold. Code is at the bottom.

Even with a small additive effect they synchronize, as shown below. They are not completely synchronized, but close enough for me.

R code

library(ggplot2)
library(arm)


plotpart <- function(xmat,fname='een.png',thin=1) {
  nstep <- ncol(xmat)
  coluse <- seq(1,nstep,by=thin)
  xmat <- xmat[,coluse]
  fin <- length(coluse)
  nitems <- nrow(xmat)
  df <- data.frame(score=as.numeric(xmat),
      cycle=rep((coluse),each=nitems),
      item =rep(1:nitems,fin))
  g<- ggplot(df,aes(x=cycle,y=score,group=item,alpha=score)) + 
      geom_line() + 
      scale_alpha_continuous(range=c(.02,.0201) )+
      theme(legend.position='none') +
      xlab('Iteration') + theme(panel.background = element_rect(fill='white'))
  png(fname)
  print(g)
  dev.off()
}

time2score <- function(x) 2*invlogit(x)-1
score2time <- function(x) logit((1+x)/2)

time_delta_score <- function(score1,score2=.99) {
  score2time(score2)-score2time(score1)
}
score_delta_time <- function(scorenow,tdelta) {
  tnow <- score2time(scorenow)
  time2score(tnow+tdelta)
}

onestep <- function(score,limit=0.99,spil=1e-4) {
  mscore <- max(score)
  dtnew <- time_delta_score(mscore,limit)
  scorenew <- score_delta_time(score,dtnew)
  maxed1 <- maxed <- scorenew==max(scorenew)#>limit*.99999999
  while(sum(maxed)>0) {
    scorenew[!maxed] <- scorenew[!maxed] + sum(maxed)*spil
    scorenew[maxed] <- 0
    maxed <- scorenew>limit
    maxed1 <- maxed | maxed1
  }
  scorenew[maxed1] <- 0
  scorenew
}

nitems <- 500
niter <- 10000
xmat <- matrix(0,nrow=nitems,ncol=niter)

score  <- time2score(runif(nitems,0,score2time(.99)))
n <- 0
while (n < niter) {
  score <- onestep(score,spil=1e-7)
  if (score[1]==0) {
    n<- n+1
    xmat[,n] <- score
  }
}

plotpart(xmat[1:250,],'een.png',thin=20)

Sync

I am listening to the audiobook Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life by Steven Strogatz which I got from Audible. Obviously a mathematical book is not ideal to listen to, but lacking illustrations I can make them myself. On top of that, you learn more from programming than listening.

He mentioned a simple system which synchronizes itself. As far as I understand a large number of items are considered. Each  item has more or less the same properties. Small steps are taken. Every step has a small additive effect, a multiplicative decrease. When an item gets over a limit, then it resets to zero and causes a small increase in the other items.

What is supposed to happen is that at low values you get steep slopes and at high values small slopes. As time progresses the small increases due to resets cause a synchronization.

To program this was not difficult, except for the initialization. If you start the items at random levels you have too many at low levels. In the end I decided to add an item at each time point for an init period. What did surprise me, was that calculating what happens actually takes less computer time than plotting it. But maybe that is because I wanted to use ggplot and a low alpha to show the bunching of items. I also discovered it is also not wise to display too complex a figure in ggplot running in Eclipse. It can freeze the whole program. I took to saving all plots rather than displaying within Eclipse.

In the init phase every item has different start and ending times. Note the black line at level 0 which consists of not yet running items

After 10000 steps it is still chaos, with a few somewhat fatter lines.

But at 50000 steps it is clear that the things are synchronizing.

Not that they are perfect after 300000 steps, but it gets more and more closer.

It is really happening. So simple and it works. Having said that, it took some tweaking to get parameters which take a while to result in sync, yet do not degenerate either. I guess it gets more easy to sync when more items are syncing, but that rapidly gets to the end of my hardware. I must say that I look forward to listening to the rest of the book. Hope to make more plots and get a bit of the math which is supposed to be behind all this.

R code

library(ggplot2)

plotpart <- function(xmat,fname='een.png',xstart=1) {
  fin <- ncol(xmat)
  nitems <- nrow(xmat)
  df <- data.frame(score=as.numeric(xmat),
      time=rep((1:fin)+xstart,each=nitems),
      item =rep(1:nitems,fin))
  g<- ggplot(df,aes(x=time,y=score,group=item,alpha=score)) + 
      geom_line() + 
      scale_alpha_continuous(range=c(.05,.0501) )+
      theme(legend.position='none')
  png(fname)
  print(g)
  dev.off()
}

nextgen <- function(x,add,shrink=.99,limit=1,spil=.0105) {
  x <- x*shrink+add
  maxed <- x>limit
  x[maxed]<- .01#add[maxed]
  x[!maxed] <- x[!maxed] + sum(maxed)*spil
  x[x>limit] <- limit
  x
}

ngenmat <- function(xstart,niter,add,spil) {
  nitems <- length(xstart)
  xmat <- matrix(0,nrow=nitems,ncol=niter)
  xmat[,1 ] <- xstart
  for ( i in 2:niter) {
    xmat[,i] <- nextgen(xmat[,i-1],add=add,spil=spil)
  }
  xmat
}

ngenmove <- function(xstart,niter,add,spil) {
  for ( i in 1:niter) {
    xstart <- nextgen(xstart,add=add,spil=spil)
  }
  xstart
}

# init a run
nitems <- 1000
niter <- 1000
xmat <- matrix(0,nrow=nitems,ncol=niter)
for (i in 1:nitems) xmat[i,i] <- runif(1,min=.01,max=.02)
add <- rnorm(nitems,.0103,.00005)
spil=.000025
for ( i in 2:niter) {
  todo <- xmat[,i-1]>0
  xmat[todo,i] <- nextgen(xmat[todo,i-1],add=add[todo],spil=spil)
}
#and show it
show <- sample(1:nrow(xmat),50)
plotpart(xmat[show,],fname='een.png',xstart=1)

# move 9000 steps
x1 <- ngenmove(xmat[,ncol(xmat)],niter=9000,add=add,spil=spil)
xm1 <- ngenmat(x1,niter=1000,add=add,spil=spil)
plotpart(xm1[show,],fname='twee.png',xstart=10000)

x2 <- ngenmove(xm1[,ncol(xm1)],niter=39000,add=add,spil=spil)
xm2 <- ngenmat(x2,niter=1000,add=add,spil=spil)
plotpart(xm2[show,],fname='drie.png',xstart=50000)

x3 <- ngenmove(xm2[,ncol(xm2)],niter=49000,add=add,spil=spil)
xm3 <- ngenmat(x3,niter=1000,add=add,spil=spil)
plotpart(xm3[show,],99000,fname='vier.png')

x4 <- ngenmove(xm3[,ncol(xm3)],niter=49000,add=add,spil=spil)
xm4 <- ngenmat(x4,niter=1000,add=add,spil=spil)
plotpart(xm4[show,],150000,fname='vijf.png')

x5 <- ngenmove(xm4[,ncol(xm4)],niter=49000,add=add,spil=spil)
xm5 <- ngenmat(x5,niter=1000,add=add,spil=spil)
plotpart(xm5[show,],200000,fname='zes.png')

x6 <- ngenmove(xm5[,ncol(xm5)],niter=99000,add=add,spil=spil)
xm6 <- ngenmat(x6,niter=1000,add=add,spil=spil)
plotpart(xm6[show,],300000,fname='zeven.png')