## ------------------------------------------------------------------------ ## ## R scripts for "Uncertainty on LTC survival probabilities and SCR" ## ## ------------------------------------------------------------------------ ## ## Description Script to produce quantites and figures for the article ## ## ------------------------------------------------------------------------ ## ## Author Frédéric Plachet & Tomas Julien ## ## ------------------------------------------------------------------------ ## ## Version 01 - 15/03/2013 ## ## 02 - 19/07/2013 ## ## 03 - 07/02/2014 ## ## ------------------------------------------------------------------------ ## ## ---------------- If necessary, rm(list=ls(all=TRUE)) ## ------------------------------------------------------------------------ ## ## Load functions ## ## ------------------------------------------------------------------------ ## ## ---------------- Graphic packages require(lattice); require(grid); require(latticeExtra) ## ---------------- ISFA couleur: isfa <- rgb(121,165,0, max=225) ## ---------------- Define logit function f.logit = function(x) { log(x/(1-x)) } f.expit = function(x) { exp(x)/(1 + exp(x)) } ## ------------------------------------------------------------------------ ## ## Load data ## ## ------------------------------------------------------------------------ ## ## --- Original data from previous article load("Mu_i.j_hat_bf.RData") q_uv <- 1-exp(-Mu_i.j_hat_bf) ## --- Disability law loi_incidence_Lx <- as.vector(as.matrix(read.table(file = "loiincidence.txt", header=T))) QxInc <- (loi_incidence_Lx[1:120] - loi_incidence_Lx[-1])/loi_incidence_Lx[1:120] ## --- Mortality law TH_Lx <- as.vector(as.matrix(read.table(file = "TH00-02.txt", header=T))) TF_Lx <- as.vector(as.matrix(read.table(file = "TF00-02.txt", header=T))) QxTab <- (((TH_Lx[1:120] - TH_Lx[-1])/TH_Lx[1:120])+((TF_Lx[1:120] - TF_Lx[-1])/TF_Lx[1:120]))/2 piInc <- pp <- vector(,120) for(i in 1:120){ piInc[i] <- QxInc[i] * prod((1-QxTab[1:i]) * (1-QxInc[1:i])) } ## ------------------------------------------------------------------------ ## ## Uncertainty ## ## ------------------------------------------------------------------------ ## ## ---------------- Function to obtain expectation and variance of remaining ## ---------------- lifetime conditionally to the disturbance MomentsLife <- function(lgt,disturbance){ lgt_disturbed <- (lgt + disturbance) q_disturbed <- f.expit(lgt_disturbed) q_disturbed[is.na(q_disturbed)] <- 1 LT <- apply(1 - q_disturbed, 2, cumprod) L0 <- pmin(1, c(1 : ncol(LT))) LT <- rbind(L0, LT) TT <- matrix(nrow=nrow(LT),ncol=ncol(LT),data=t(c(1:nrow(LT)))) - 1 mu <- colSums(LT[1:nrow(LT),]) - 1 sigma <- (2*colSums(TT * LT[1:nrow(LT),]) - mu^2)^.5 cbind(mu,sigma) } f.1 <- function(c){ MomentsLife(lgt,c)[,1] } f.2 <- function(c){ MomentsLife(lgt,c)[,2] } f.3 <- function(x){ (quantile(x,.95) - mean(x)) / mean(x) } ## ---------------- Simulation of empirical distribution of remaning lifetime ## ---------------- expectancy volatility <- seq(.01,.5,.01) pIncVec <- piInc[71:91] lgt <- f.logit(q_uv) delta <- matrix(, length(volatility), ncol(lgt)) esp_c <- vol_c <- vector("list",length(volatility)) for (k in 1:length(volatility)){ disturbance <- rnorm(n=5000, mean=0, sd=volatility[k]) temp_1 <- apply(as.matrix(disturbance),1,f.1) temp_2 <- apply(as.matrix(disturbance),1,f.2) temp_1[is.na(temp_1)] <- 0 temp_2[is.na(temp_2)] <- 0 esp_c[[k]] <- vol_c[[k]] <- matrix(,ncol(lgt),ncol(temp_1)) for(j in 1:ncol(temp_1)){ for(i in 1: ncol(lgt)){ esp_c[[k]][i,j] <- sum(temp_1[i:ncol(lgt),j]*pIncVec[i:ncol(lgt)]) vol_c[[k]][i,j] <- sum(temp_2[i:ncol(lgt),j]*pIncVec[i:ncol(lgt)]) } } delta[k,] <- apply(esp_c[[k]],1,f.3) } esp_m <- matrix(,21,20) for(i in 1:20){ esp_m[,i] <- apply(esp_c[[i]],1,mean) } ## ---------------- Figure 1 surface_plot_color = function(xx, zexpr, axis, cc){ jet.colors <- colorRampPalette(c("white", cc)) color <- jet.colors(100) wireframe(matrix(xx,nrow(xx),ncol(xx)), zlab=list(zexpr, rot=90), ylab=list("Age"), xlab=list("Volatility"), aspect=c(2,1), drape = T, colorkey = F, screen = list(z = 20, x = -75), scales = list(arrows = F, y = list(at = c(seq(1,ncol(xx), by = 2)), labels = c(seq(axis[3],axis[4],by=2)), cex=.8), x = list( at = c(seq(1,21, by = 2)), labels = c(seq(axis[1],axis[2], by = .02)),cex=.75)), col.regions = color) } surface_plot_color(delta[1:21,], expression(delta), c(.01,.21,70,90),isfa) ## ---------------- More accurate volatility <- seq(7,11,.01)/100 pIncVec <- piInc[71:91] lgt <- f.logit(q_uv) delta.5 <- matrix(, length(volatility), ncol(lgt)) esp_c.5 <- vol_c.5 <- vector("list",length(volatility)) for (k in 1:length(volatility)){ disturbance <- rnorm(n=5000, mean=0, sd=volatility[k]) temp_1 <- apply(as.matrix(disturbance),1,f.1) temp_2 <- apply(as.matrix(disturbance),1,f.2) temp_1[is.na(temp_1)] <- 0 temp_2[is.na(temp_2)] <- 0 esp_c.5[[k]] <- vol_c.5[[k]] <- matrix(,ncol(lgt),ncol(temp_1)) for(j in 1:ncol(temp_1)){ for(i in 1: ncol(lgt)){ esp_c.5[[k]][i,j] <- sum(temp_1[i:ncol(lgt),j]*pIncVec[i:ncol(lgt)]) vol_c.5[[k]][i,j] <- sum(temp_2[i:ncol(lgt),j]*pIncVec[i:ncol(lgt)]) } } delta.5[k,] <- apply(esp_c.5[[k]],1,f.3) print(volatility[k])} ## ---------------- Obtain mu mu <- vector(,ncol(lgt)) for(i in 1: ncol(lgt)){ mu[i] <- sum((MomentsLife(lgt,0)[,1])[i:ncol(lgt)] * pIncVec[i:ncol(lgt)]) } ## ---------------- Obtain volatility expert vol.expert <- ind.exp <- vector(,21) for (i in 1:21){ mu5 <- sum(5 * pIncVec[i:ncol(lgt)]) ind.exp[i] <- max(which(delta.5[,i] <= (mu5/mu)[i])) vol.expert[i] <- volatility[ind.exp[i]] } ## ---------------- Table 1 cbind(round(as.matrix(mu),2), round(mu5/mu * 100,2), as.matrix(vol.expert*100)) ## ---------------- Obtain quantile 99,5% volatility <- seq(.01,.5,.01) q995 = function(v){ quantile(v, .995) } esp_995 <- matrix(,length(volatility),21) mu <- vector(,ncol(lgt)) for(i in 1: ncol(lgt)){ mu[i] <- sum((MomentsLife(lgt,0)[,1])[i:ncol(lgt)] *pIncVec[i:ncol(lgt)]) } for(i in 1:length(volatility)){ esp_995[i,] <- apply(esp_c[[i]],1,q995) / mu } volatility <- seq(7,11,.01)/100 esp_995.exp <- vector(,21) for(i in 1:21){ esp_995.exp[i] <- (apply(esp_c.5[[ind.exp[i]]],1,q995) / mu)[i] } ## ---------------- Comparison with QIS5 tt <- vector(,ncol(lgt)) for(i in 1: ncol(lgt)){ tt[i] <- sum((MomentsLife(f.logit(q_uv * .8),0)[,1])[i:ncol(lgt)] *pIncVec[i:ncol(lgt)]) } esp_80 <- tt / mu ## ---------------- Figure 2 volatility <- seq(.01,.5,.01) plot(c(70:90), esp_80,ylim=c(1,1.6),type="b",pch=16, xlim=c(70,92.3), xlab="Age of occurrence",ylab=expression(paste(rho[99.5], " / BEL"))) for (i in 1:20){ points(c(70:90),esp_995[i,],col=rainbow(20)[i],type="l") } for (i in 1:13){ text(90.8,esp_995[i,21],labels = expression(paste(sigma[epsilon],"=")),cex=.7, col=rainbow(20)[i] ) text(91.9,esp_995[i,21],labels = volatility[i],cex=.7, col=rainbow(20)[i] ) } text(90.8,esp_995[14,21]+0.005,labels = expression(paste(sigma[epsilon],"=")),cex=.7, col=rainbow(20)[14] ) text(91.9,esp_995[14,21]+0.005,labels = volatility[14],cex=.7, col=rainbow(20)[14] ) for (i in 15:16){ text(90.8,esp_995[i,21],labels = expression(paste(sigma[epsilon],"=")),cex=.7, col=rainbow(20)[i] ) text(91.9,esp_995[i,21],labels = volatility[i],cex=.7, col=rainbow(20)[i] ) } text(90.8,esp_995[17,21]+0.006,labels = expression(paste(sigma[epsilon],"=")),cex=.7, col=rainbow(20)[17] ) text(91.9,esp_995[17,21]+0.01,labels = volatility[17],cex=.7, col=rainbow(20)[17] ) text(90.8,esp_995[18,21]-0.006,labels = expression(paste(sigma[epsilon],"=")),cex=.7, col=rainbow(20)[18] ) text(91.9,esp_995[18,21]-0.006,labels = volatility[18],cex=.7, col=rainbow(20)[18] ) text(90.8,esp_995[19,21]+0.006,labels = expression(paste(sigma[epsilon],"=")),cex=.7, col=rainbow(20)[19] ) text(91.9,esp_995[19,21]+0.006,labels = volatility[19],cex=.7, col=rainbow(20)[19] ) text(90.8,esp_995[20,21],labels = expression(paste(sigma[epsilon],"=")),cex=.7, col=rainbow(20)[20] ) text(91.9,esp_995[20,21],labels = volatility[20],cex=.7, col=rainbow(20)[20] ) points(c(70:90),esp_995.exp,type="b",col=2, pch=16) ## ---------------- Application to SCR FDR = function(x, m, v){ y <- (x-m)/sqrt(v) sum(apply(as.matrix(y),1,pnorm)) / length(m) } InverseFDR = function(q, m, v){ epsilon=10^-2; xR <- 10 * max(m); xL <- 0 while ((xR - xL)>epsilon) { xM <- (xR + xL) / 2 zL <- FDR(xL, m, v) - q zM <- FDR(xM, m, v) - q if (zL * zM > 0) { xL = xM } else { xR = xM } } xR } SCR = function(nIndi,e,v){ q <- InverseFDR(.995, nIndi * e, nIndi * (v^2)) bel <- nIndi * mean(e) khi <- (q - bel) / bel alpha <- .06 d <- mean(e) khi / (1 - alpha * d * khi) } nn <- c(100,250,500,750,1000,2500,5000,7500,10000) esp_SCR <- matrix(,length(nn),21) for(j in 1:21){ for (n in 1:length(nn)){ esp_SCR[n,j] <- SCR(nn[n], esp_c[[9]][j,], vol_c[[9]][j,]) } } surface_plot_color = function(xx, zexpr, axis, cc){ jet.colors <- colorRampPalette(c("white", cc)) color <- jet.colors(100) wireframe(matrix(xx,nrow(xx),ncol(xx)), zlab=list(zexpr, rot=90), ylab=list("Age"), xlab=list("Number of individuals"), aspect=c(2,1), drape = T, colorkey = F, screen = list(z = 20, x = -75), scales = list(arrows = F, y = list(at = c(seq(1,ncol(xx), by = 5)), labels = c(seq(axis[3],axis[4],by=5)), cex=.7), x = list( at = c(seq(1,9, by = 1)), labels = c(nn),cex=.8)), col.regions = color) } surface_plot_color(esp_SCR, "SCR / BEL", c(.01,.2,70,90),isfa)