--- imach/src/simulation.R	2023/07/14 11:25:16	1.1
+++ imach/src/simulation.R	2023/07/19 15:21:39	1.2
@@ -1,3 +1,6 @@
+# $Id: simulation.R,v 1.2 2023/07/19 15:21:39 brouard Exp $
+#  $State: Exp $
+#  
 # Simulating data for IMaCh tests
 # Feinuo Sun and Nicolas Brouard (INED, July 2023)
 # 
@@ -9,7 +12,8 @@ if(!require("tidyverse")){
   }
 # LIBRARIES
 library(haven)
-samplesize<- 10000
+#samplesize<- 10000
+samplesize<- 50000
 
 # Which Life table?
 # Simulating a gompertz law of mortality: slope of line mu_m (about 9% increase per year)
@@ -26,7 +30,8 @@ exp(exp(mum*(65-am)))/mum*expint(exp(mum
 e65 <- exp(exp(mum*(65-am)))/mum*expint(exp(mum*(65-am)))
 
 e65 <- exp(exp(mum*(65-am)))/mum*expint(exp(mum*(65-am)))
-# Life table from 0
+e65
+                                        # Life table from 0
 l <- function(a,am,mum){
     exp(-exp(mum*(a-am)))/exp(-exp(-mum*am))
 }
@@ -59,18 +64,25 @@ mum
 l65bisinv <- function(la){
     84.8506+log(exp(0.09*(65-84.8506)) -log(la))/0.09
 }
-l65bisinv(0.001)
+l65bisinv(0.0001)
 l65bisinv(1)
 
 set.seed(128)
+zeroone<-runif(samplesize,min=0.00000001,max=1)
 zeroone<-runif(samplesize,min=0.001,max=1)
 lifelength <- lapply(zeroone,l65bisinv)
+#slife<- sort(unlist(lifelength),decreasing=TRUE)
+#?sort
+#head(slife)
+#mean(unlist(lifelength)-65) #17.99482
+#head(lifelength,10)
+#class(lifelength)
 
 #lifelength <- rnorm(n=samplesize, mean=85, sd=16)
 
 set.seed(124)
 # First interview 4th semester
-monthinterview <- runif(samplesize, min=9,max=12)
+monthinterview <- runif(samplesize, min=10,max=13)
 st_98 <- rbinom(n=samplesize, size=1, prob=0.5)+1
 # state 2
 st_00 <- rbinom(n=samplesize, size=1, prob=0.2)+1
@@ -79,28 +91,34 @@ popage65110in1998<- runif(samplesize, mi
 #gender
 ragender <- rbinom(n=samplesize, size=1, prob=0.56)
 
+r98iwmid <- 1998 + monthinterview/12
 yrinterview1 <- floor(r98iwmid)
-monthinw1 <- floor((monthinterview - floor(monthinterview))*12)+1
+#monthinw1 <- floor((monthinterview - floor(monthinterview))*12)+1
+monthinw1 <- floor(monthinterview)
+
 
 int_98 <- paste0(monthinw1,"/",yrinterview1)
 
-r98iwmid <- 1998 + monthinterview/12
 rabdate <- r98iwmid - popage65110in1998 
 birthyr <- floor(rabdate)
 monthdb <- floor((rabdate - floor(rabdate))*12)+1
 brt <- paste0(monthdb,"/",birthyr)
 # date of death for dropping cases
-raddate <- rabdate + lifelength
+#raddate <- rabdate + lifelength
+raddate <- rabdate + unlist(lifelength)
 dateinterview2 <- 1998 + 2 + monthinterview/12
 #head(cbind(r98iwmid,dateinterview2,raddate),70)
 # people whose death occured before the interview will be dropped (radnate)
 radndate <- if_else(raddate < r98iwmid,  NA, raddate )
+head(cbind(r98iwmid,dateinterview2,raddate,radndate,lifelength),n=5299)
 #head(cbind(r98iwmid,dateinterview2,raddate,radndate,lifelength),70)
 # in order to avoid date of death known after last wave and potential bias
 # people whose death will occur after last interview will have an unkwonn (99/9999) date of death
 lastinterview<- dateinterview2
-radldate<- if_else(radndate > dateinterview2,  9999, radndate )
-head(cbind(r98iwmid,dateinterview2,raddate,radndate,lifelength,radldate),70)
+radldate<- if_else(radndate > dateinterview2+1/12,  9999, radndate )
+head(cbind(r98iwmid,dateinterview2,raddate,radndate,lifelength,radldate),n=5299)
+radldate<- if_else(dateinterview2+1/12 > radndate & radndate > dateinterview2, radndate , radldate )
+head(cbind(r98iwmid,dateinterview2,raddate,radndate,lifelength,radldate),n=5299)
 ddtyr <- if_else((!is.na(radldate) & radldate ==9999),  9999, floor(radldate))
 monthdd <- if_else((!is.na(radldate) & radldate ==9999),  99,floor((radldate - floor(radldate))*12)+1)
 #head(cbind(r98iwmid,dateinterview2,raddate,lifelength,radldate,ddtyr,monthdd),70)
@@ -114,6 +132,7 @@ ageatinterview1 <- r98iwmid - rabdate
 yrinterview2 <- floor(dateinterview2)
 monthinw2 <- floor((dateinterview2 - floor(dateinterview2))*12)+1
 int_00 <- paste0(monthinw2,"/",yrinterview2)
+#head(cbind(r98iwmid,dateinterview2,int_00,raddate,radndate,lifelength,radldate),70)
 ageatinterview2 <- dateinterview2 - rabdate
 # state 2
 st_00 <- if_else(raddate < dateinterview2, 3, st_00)