[BACK]Return to simulation.R CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/simulation.R between versions 1.1 and 1.2

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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>