timriffe

# install.packages("devtools", dependencies = TRUE)
# library(devtools)
# install_github("DemogBerkeley", subdir = "DemogBerkeley", username = "UCBdemography")

# read in data:
library(DemogBerkeley)

mltper <- readHMDweb("USA","mltper_1x1",password=pw,username=us)
Exp    <- readHMDweb("USA","Exposures_1x1",password=pw,username=us)

timriffe

ternary2xy <- function(v1,v2,v3){
    list(x = (2 * v2 + v3) / ((v1+v2+v3) * 2),
         y = (sqrt(3) * v3)/ ((v1+v2+v3) * 2)
 )
}
v  <- seq(0,1,by=.05) 
V  <- expand.grid(v1=v,v2=v,v3=v)
V  <- V / rowSums(V)
xy <- ternary2xy(V[,1],V[,2],V[,3])

timriffe

# Da Yu, Ding's data reshaped in the spreadsheet and copied it into R.
# tip: if you have an object in R, you can make it print conveniently to
# the console using:
# dput(data) # data = name of object
# result:
data <- structure(list(Age = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 
                        10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 
                        23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 
                        36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 
                        49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 

timriffe

# trick is to identify all unique point dyads. 
# One could calculate everything in each direction and take the average as well...
# here we just go in one direction for each dydad and assume it's valid.

# example data:

DAT <- read.table(text=
"2.018233 41.31934 '08301-01-015' 2.01823300593231 41.3193364170593
2.010912 41.31946 '08301-02-021' 2.01091150024005 41.3194594380338
2.026979 41.31763 '08301-01-001' 2.02697878055703 41.3176301571496

timriffe

'%==%' <- function(x,y){
  x == y & !is.na(x) & ! is.na(y)
}

'%>%' <- function(x,y){
  x > y & !is.na(x) & ! is.na(y)
}

'%<%' <- function(x,y){
  x < y & !is.na(x) & ! is.na(y)

timriffe

# I usually generate death counts for each age using rpois() and m(x) * Exposure(x) as the mean
# then for each of the N simulation runs, I calculate e0 from the lifetable, and you can look at the
# properties of the distribution of e0. This script takes a different approach, and I'm not sure how 
# the two approaches compare yet. In the present simulation, I take a vector of m(x) as canonical from the outset,
# and instead what we simulate are lifespans under this set of m(x). To do this, assume that the rate in
# each single age is the hazard, and that it is constant within each single age. We can then generate 
# random draws from the exponential distribution. These are theoretical survival times assuming you're
# exposed only to this rate. If the draw is > 1, then we say the person survived the age interval. The age
# at death is determined by the first interval in which the person didn't make it to the end, so draw < 1.

timriffe

# this approach is asymptotically identical to e0sim1.R. It generates lifespans via sampling ages, using d(x) as weights
# this is different from the other lifespan script, where individuals need to
# survive through age classes sequentially.

mxi <- c(0.00339193594631364, 0.000218176772567409, 0.000112225116810933, 
  6.28291177575766e-05, 0.000101363542293164, 6.51192928664136e-05, 
  0.000149171546489114, 9.1324237373929e-05, 2.37282051423152e-05, 
  0.000147879337023085, 9.65951840540929e-05, 0.000117234660071609, 
  0.000148808909648065, 0.000163764625306839, 0.000149412048660226, 
  0.000213094661860419, 0.00023027681399857, 0.000282336781997834, 

timriffe

ES <- local(get(load(url("http://biogeo.ucdavis.edu/data/gadm2/R/ESP_adm2.RData"))))

library(sp)
plot(ES, main = "antes")
Cani <- which(ES@data$NAME_1 == "Islas Canarias")
L <- length(ES@polygons[[45]]@Polygons)
for (i in 1:L){
  ES@polygons[[Cani]]@Polygons[[i]]@coords <- cbind(ES@polygons[[Cani]]@Polygons[[i]]@coords[,1]+6,ES@polygons[[Cani]]@Polygons[[i]]@coords[,2]+5)
}
plot(ES, main = "despues")

timriffe

# R code for 4th of July jet lag blog post
# Tim Riffe
############################################

# library(devtools)
# install_github("timriffe/TR1/TR1/HMDHFDplus")
library(HMDHFDplus)

# it'll ask you to give your HMD password and username in the console (no quotes)
Mx <- readHMDweb("USA", "Mx_1x1")

timriffe

# matrix of marriage counts with female ages in columns and male ages in rows
MAR <- matrix(c(4145,24435,8140,1865,1655,54515,45010,15030,80,6735,20870,19530,5,920,5435,42470),ncol=4)
rownames(MAR) <- colnames(MAR) 	<- c("15-19","20-24","25-29","30-60")
# unmarried males and females at start of the year
unMARf 	<- c(254876,147705,61804,415497)
unMARm 	<- c(265755,199437,114251,429655)
PanmicticRates <- function(MAR,unMARf,unMARm){
	# allocate P, array of stacked component counts
	P <- array(0,dim=c(dim(MAR),min(dim(MAR))))
	Pi <- Ri <- MAR