Cacher is the code snippet organizer for pro developers

#Sample size
N=10000
#Simulate a population sample, ages ranging from 20-40
age1 <- runif(N,20,40)
#Where trauma follows an exponential distribution (most people have little trauma)
trauma1 <- rexp(N,rate=.5)

#Model methylation age as real age + random noise (3 years deviation on average) + Acceleration due to trauma (.1 years for every unit increase in trauma)
age1_dnam <- age1 + rnorm(N,0,3) + .5*trauma1

#Model methylation age as a function of age, get the residuals
age1_resids <- lm(age1_dnam ~ age1 )$residual
summary(lm(age1_dnam ~ age1 ))$coefficients[2,1] #Beta should be close to 1

#Association the residuals to trauma
summary(lm(age1_resids ~ trauma1))

##Traumatized sample

#Take a sample of data with the same ages
age2 <-  age1

#But who have been sampled on the upper range of the trauma distribution (Median or above)
trauma2 <- runif(N,median(trauma1),max(trauma1))

#Model acceleration in the same way
age2_dnam <- age2 + rnorm(N,0,3) + .5*trauma2

#Model methylation age as a function of age, get the residuals
age2_resids <- lm(age2_dnam ~ age2 )$residual
summary(lm(age2_dnam ~ age2 ))$coefficients[2,1] #Beta should be close to 1


#Association the residuals to trauma
summary(lm(age2_resids ~ trauma2))




pdf('age1plots.pdf',7,7)
plot(age1,age1_dnam)
plot(age2,age2_dnam)

dev.off()