# "SIZE GROWTH." This program is designed to compare small, fragmentary (i.e. fossil) samples with more complete ones, to test the null hypothesis that patterns of size change (i.e. growth) are the same between the two groups. The resampling procedure runs a few minutes for 10,000 iterations, and around 30 minutes for 50,000. The number of iterations you want to run will depend on your sample sizes.

# This is set up to read a dataset in .csv format, with individuals/observations' names in the first column, individuals' ages (or dental stages) in the second column, and the variables used to calculate individuals' size (calculated as a geometric mean) occupying subsequent columns

setwd("")	# set the working directory, from which data can be read and to which output is saved
  data = data.frame(read.csv(".csv")) 	# enter the "filename.csv"
  	row.names(data)=data[,1]; data=data[,2:ncol(data)]
    lib=data[x:x,]; rob = data[y:y,]	# 'lib' is your large comparative dataset, 'rob' the small (fossil) dataset. Enter the rows occupied by each sample x and y. 
    p.lib = data.frame(t(lib)); p.rob=data.frame(t(rob))	# transposes the 2 groups or subsamples, to be manipulated below

# Geometric mean f(x)
GM = function(r1)  {
  prod(r1)^(1/nrow(r1))
	}

#### randomly grab 2 diff-aged robusts, and 2 humans of same ages
#	 First declare what variables will be extracted in each resampling
D = NA    # test statistic, measuring species difference in size change between dental stages
	# if D = 0 then this means no difference, positive/negative values mean more size change in the small-sample group than the other comparative sample
	# when resampled, if 0 is within 95% of the resampled distribution, the null hypothesis of no change cannot be rejected
k = NA    # number of traits shared by subsampled individuals
a1 = NA; a2 = NA   # ages in the comparison. Later can be used to examine differences between specific age groups. a1 is older and a2 is the younger
f1 = NA; f2 = NA		# these are the individuals in the small (i.e. fossil) sample. 1-2 correspond to a1-2 above
h1 = NA; h2 = NA 		# these are the individuals in the larger comparitive sample.

#    Now resample relative growth stats 'D,' describing species' differences in the amount size change (i.e. growth) between dental stages
for (i in 1:n)  {	# n being the number of resamplings
r1 = sample(p.rob,1); r2 = sample(p.rob[,p.rob[1,]!=r1[1,1]],1) # randomly select 2 differently-aged individuals of one group
    l1 = sample(p.lib[,p.lib[1,]==r1[1,1]],1)   # then 2 from the second group of same ages as the pair in the first group
    l2 = sample(p.lib[,p.lib[1,]==r2[1,1]],1)
  sub = data.frame(na.omit(cbind(r1,r2,l1,l2)))		# creates matrix of resampled subset, omits traits not shared by all 4
	k[i] = nrow(sub)-1;   # number of traits all share
  s2 = sub[2:nrow(sub),]  # get rid of Age for subsequent GM calculations
    # calculate test statistic
  if (k[i] == 0) {D[i] = NA}  else {	#no traits means D cannot be calculated, this is stored as NA to be omitted later
    if (sub[1,1]>sub[1,2])  {D[i] = (GM(s2[1])/GM(s2[2])) - (GM(s2[3])/GM(s2[4]))}	# calculate ratio of oldest / youngest
      else {D[i] = (GM(s2[2])/GM(s2[1])) - (GM(s2[4])/GM(s2[3]))}
  # extract information about subsample which can be explored later; redundant resamplings can be omitted
    if (sub[1,1]<sub[1,2])  {a1[i] = sub[1,1]; a2[i]=sub[1,2]} else {a1[i]=sub[1,2]; a2[i]=sub[1,1]}
    if (sub[1,1]<sub[1,2])  {f1[i]=names(sub[1]); f2[i]=names(sub[2])} else {f1[i]=names(sub[2]); f2[i]=names(sub[1])}
    if (sub[1,1]<sub[1,2])  {h1[i]=names(sub[3]); h2[i]=names(sub[4])} else {h1[i]=names(sub[4]); h2[i]=names(sub[3])}
    }
}#                                  FIN!

###### OUTPUT: all of the resampled info ######
da = a2-a1  # diff in age groups b/w pairs; a2 should always be greater than a1
total = data.frame(cbind(D,k,da,a1,a2,f1,f2,h1,h2))  # all the resampled data, including NAs
  tot = data.frame(unique(na.omit(total)))		# this gets rid of all NAs and redundant resamplings. Can remove unique() if desired
  write.csv(tot,file="insertname.csv") 			# save a copy of the output for a rainy day; if you do not save these output data you will need to run the analysis each time you restart R.
  
#R sometimes doesn't create the 'total' matrix above so that each column is in the right right class, but will only treat data easily if re-read from a .csv file, hence...
d2 = data.frame(read.csv("insertname.csv"));d2 = d2[,2:5] # this is the same as tot, above

# this output can be examined however you want from here. Here's how you get a p value for if D is significantly greater than 0
length(d2$D[d2$D<=0])/length(d2$D)	# neat-o!

#You may be interested in how D behaves across all age groups. Here's a pic:
hist(d2$D, col=4,main = "All possible resampled D", xlab = "D",breaks=100,border=4)
  abline(v=0,lty=4,lwd=4)
  
# Or you could look only at those D calculated between individuals in adjacent age groups - is sort of like comparing the deviation (or residuals) of each group from a linear model of growth rates across ontogeny
one = d2[d2$da==1,]
	plot(one$D~one$a1)	# this plots the D statistic from each age group comparison (ages or stages 1-2, 2-3...)