breuer.R

richard — Jan 9, 2014, 8:41 AM

## This script: http://rpubs.com/gill1109/11774

## Reference: "Empirical Patterns in Google Scholar Citation Counts"
## Peter T. Breuer and Jonathan P. Bowen
## Preprint available from Peter Breuer (University of Birmingham)
## through ResearchGate.net and/or Academia.edu

## Breuer plot: log positive citation count against square root of reverse rank number
## (most cited has reverse rank 0, next most cited 1, ...)
## Sometimes need a power different from 1/2

## My data (extracted from Google Scholar) 
## http://www.math.leidenuniv.nl/~gill/citations.txt

## Previous analyses: http://rpubs.com/gill1109/11755

## The Breuer plot corresponds to the following model.
## Suppose U ~ Unif(0,1). Define X = c(1 - U^a), in my case, c = 7.5, a = 1/2.
## Then my positive citation counts are distributed as a sample of size N = 127
## from the distribution of exp(X), rounded up to positive whole numbers.

## After fitting this straight line to my data I simulate 9 times from the estimated model
## and show 9 Breuer plots (which of course give beautiful straight lines)
## and 9 Normal probability plots (which show a good fit to the log normal).

data <- read.csv("http://www.math.leidenuniv.nl/~gill/citations.txt", 
                 sep = "$", 
                 stringsAsFactors = FALSE, 
                 header = FALSE, 
                 blank.lines.skip = FALSE)
V1 <- data$V1
counts <- as.numeric(V1[seq(from = 4, by = 7, to = 1100)])
counts[is.na(counts)] <- 0


## log positive citation counts, decreasing order

lpcounts <- sort(log(counts[counts>0]), decreasing = TRUE)

## create Breuer plot

plot(sqrt(ppoints(lpcounts)), lpcounts,
    main = "Breuer plot, my data",
    xlab = "sqrt(1-r/(n+1))",
    ylab = "log positive counts",
    sub = "r  runs from 1 to n (reverse rank)",
    xlim = c(0, 1))

## Fit a straight line to body of data - determined by first, third quartile

## I'd like the straight line to go through the point (1, 0)
## So if it's y = a + bx I want a = -b
## Adjusting the quartiles line a bit to get equality I settle on a = 7.5

quartiles <- quantile(lpcounts, probs = c(0.75, 0.25))
quartiles

  75%   25% 
3.726 1.099 

Uquartiles <- sqrt(c(0.25,0.75))
Uquartiles

[1] 0.500 0.866

coefs <- lm(quartiles ~ Uquartiles)$coef
coefs

(Intercept)  Uquartiles 
      7.314      -7.177 

abline(coefs, col = "blue")
abline(7.5, -7.5, col = "red")
abline(7.2, -7.2, col = "green")

legend(0.1, 2, 
    c("quartile fit", "both coeffs 7.2", "both coeffs 7.5"), 
    col = c("blue", "green", "red"), 
    lty = 1)

qqnorm(lpcounts,
    main = "Normal probability plot, my data (log positive counts)",
    xlab = "standard normal quantiles",
    ylab = "log positive counts")
qqline(lpcounts)

## Now simulate (9 times) from the estimated model according to the Breuer plot:
## log positive count = lpc = 7.5 ( 1 - sqrt U) where U ~ Unif(0,1)
## Equivalently: (1 - (log pos cnt)/7.5)^2 ~ Unif(0,1)
## However to add realism at low counts I'll transform the continuous lpc to discrete
## by the transformation discrete lpc = log(ceiling(exp(continuous lpc)))

set.seed(1234)
N <- length(lpcounts)

par(mfrow = c(3, 3))
for (i in 1: 9) {
    sim <- log(ceiling(exp(7.5 * (1 - sqrt(sort(runif(N)))))))
    plot(sqrt(ppoints(sim)), sim,
        main = "Breuer plot, simulated data",
        xlab = "sqrt(1-r/(n+1))",
        ylab = "log positive counts",
        sub = "r  runs from 1 to n (reverse rank)",
        xlim = c(0, 1))
    abline(7.5, -7.5)
}

par(mfrow = c(3, 3))
set.seed(1234)
for (i in 1: 9) {
    sim <- log(ceiling(exp(7.5 * (1 - sqrt(sort(runif(N)))))))
    qqnorm(sim)
    qqline(sim)
}

par(mfrow = c(1, 1))

## We noticed that the mean and standard devation of the normal seemed about equal.
## Is this a coincidence?
## Below I plot the QQ plot of the Breuer distribution (y-axis)
## against the standard normal (x-axis).
## The straight line is the tangent line at the medians.
## It is approximately y = 2.2 + 2.1 x
## So with this normal approximation, mu = 2.2 and sigma = 2.1
## Changing the power "1/2" a little, I could get them equal
## but changing the upper bound "7.5" does not help.


qBreuer <- 7.5 * (1 - sqrt(1-ppoints(100)))

qNorm <- qnorm(ppoints(100))


plot(qNorm, qBreuer, type="b", pch="*",
    main = "QQ plot of log-Breuer against log-normal",
    xlab = "standard normal quantiles",
    ylab = "Breuer(7.5, 1/2) quantiles",
    sub = "straight line = tangent at median: y = 2.2 + 2.1 x")

points(qNorm[50], qBreuer[50], col="red")

a <- 7.5*(1-sqrt(1/2))
b <- 7.5/(2*sqrt(pi))

a

[1] 2.197

b

[1] 2.116

abline(a, b)