R/MFDM.R

In ftsa: Functional Time Series Analysis

MFDM <- function(mort_female, mort_male, mort_ave, percent_1 = 0.95, percent_2 = 0.95, fh, level = 80,
                 alpha = 0.2, MCMCiter = 100, fmethod = c("auto_arima", "ets"), BC = c(FALSE, TRUE), lambda)
{
   	fmethod = match.arg(fmethod)
  	if(BC == TRUE)
  	{
    		mort_female = BoxCox(mort_female, lambda)
    		mort_male = BoxCox(mort_male, lambda)
    		mort_ave = BoxCox(mort_ave, lambda)
    }
  	
  	# mean
  	
  	mort_femalemean = rowMeans(mort_female)
  	mort_malemean   = rowMeans(mort_male)
  	mort_avemean    = rowMeans(mort_ave)
  	
  	# de-mean
  	
  	mort_avedemean = mort_ave - matrix(rep(mort_avemean, ncol(mort_female)), ncol = ncol(mort_female))
  	mort_femaledemean = mort_female - matrix(rep(mort_femalemean, ncol(mort_female)), ncol = ncol(mort_female))
  	mort_maledemean = mort_male - matrix(rep(mort_malemean, ncol(mort_female)), ncol = ncol(mort_female))
  	
  	# svd
  	
  	W_1 = scale(t(mort_female), scale = FALSE)
  	W_2 = scale(t(mort_male), scale = FALSE)
  	mort_avesvd = svd(t(mort_avedemean))
  	ncomp_1 = which.min(abs(cumsum(mort_avesvd$d^2/sum(mort_avesvd$d^2)) - percent_1))
  	first_percent = (mort_avesvd$d[1])^2/sum(mort_avesvd$d^2)
  	basis_ave = mort_avesvd$v[,1:ncomp_1]
  	score_ave = t(mort_avedemean) %*% basis_ave
  	decomp_ave = basis_ave %*% t(score_ave)
  	
  	mort_femalediff = t(W_1) - decomp_ave
  	mort_malediff = t(W_2) - decomp_ave
  	
  	# svd (2nd time) (female)
  	
  	mort_femalesvd = svd(t(mort_femalediff))
  	ncomp_female = which.min(abs(cumsum(mort_femalesvd$d^2/sum(mort_femalesvd$d^2)) - percent_2))
  	female_percent = (mort_femalesvd$d[1])^2/sum(mort_femalesvd$d^2)
  	basis_female = mort_femalesvd$v[,1:ncomp_female]
  
  	mort_malesvd = svd(t(mort_malediff))
  	ncomp_male = which.min(abs(cumsum(mort_malesvd$d^2/sum(mort_malesvd$d^2)) - percent_2))
      male_percent = (mort_malesvd$d[1])^2/sum(mort_malesvd$d^2)
  	basis_male = mort_malesvd$v[,1:ncomp_male]
  
  	psi_1 = as.matrix(basis_ave)
  	psi_2 = as.matrix(basis_female)
  	psi_3 = as.matrix(basis_male)
  
  	N_subj = ncol(mort_female)
  	N_obs  = nrow(mort_female)
  	dim_space_b = ncomp_1
  	dim_space_w = ncomp_female
  	dim_space_f = ncomp_male
  	data_A = list("N_subj", "N_obs", "dim_space_b", "dim_space_w", "dim_space_f", "W_1", "W_2", "psi_1", "psi_2", "psi_3")
  
  	inits_A = function()
  	{
    		list(ll_w = c(rep(1, ncomp_female)), ll_f = c(rep(1, ncomp_male)), ll_b = c(rep(1, ncomp_1)), 
       			taueps_1 = 1, taueps_2 = 1)
  	}
  		
  	zi = matrix(NA, N_subj, dim_space_w)
  	xi = matrix(NA, N_subj, dim_space_b)
  	fi = matrix(NA, N_subj, dim_space_f)
  	ll_b = vector("numeric", dim_space_b)
  	ll_w = vector("numeric", dim_space_w)
  	ll_f = vector("numeric", dim_space_f)	
  	
  	inprod = function(a, b)
  	{
  		  return(matrix(a, nrow = 1) %*% matrix(b, ncol = 1))
  	}
  		
  	popmodel = function()
  	{
  		for(i in 1:N_subj)
  		{
  			for(t in 1:N_obs)
  			{
  				W_1[i,t] ~ dnorm(m_1[i,t], taueps_1)
  				W_2[i,t] ~ dnorm(m_2[i,t], taueps_2)
  
  				m_1[i,t] <- X[i,t] + U_1[i,t]
  				m_2[i,t] <- X[i,t] + U_2[i,t]
  
  				X[i,t] <- inprod(xi[i,], psi_1[t,])
  				U_1[i,t] <- inprod(zi[i,], psi_2[t,])
  				U_2[i,t] <- inprod(fi[i,], psi_3[t,])				
  			}
  			for(k in 1:dim_space_b)
  			{
  				xi[i,k] ~ dnorm(0.0, ll_b[k])
  			}
  			for(l in 1:dim_space_w)
  			{
  				zi[i,l]  ~ dnorm(0.0, ll_w[l])
  			}
  			for(j in 1:dim_space_f)
  			{	
  				fi[i,j] ~ dnorm(0.0, ll_f[j])
  			}
  		}	
  		for(k in 1:dim_space_b)
  		{
  			ll_b[k] ~ dgamma(1.0E-3, 1.0E-3)
  			lambda_b[k] <- 1/ll_b[k]
  		}
  		for(l in 1:dim_space_w)
  		{
  			ll_w[l] ~ dgamma(1.0E-3, 1.0E-3)
  			lambda_w[l] <- 1/ll_w[l]
  		}
  		for(j in 1:dim_space_f)
  		{
  			ll_f[j] ~ dgamma(1.0E-3, 1.0E-3)
  			lambda_f[j] <- 1/ll_f[j]
  		}
  		
  		# prior
  
  		taueps_1 ~ dgamma(1.0E-3, 1.0E-3)
  		taueps_2 ~ dgamma(1.0E-3, 1.0E-3)
  	}	
  	if(requireNamespace("R2jags", quietly = TRUE)) 
  	{
  		bugs_chose = R2jags::jags(data = data_A, inits = inits_A, model.file = popmodel,
  		                    parameters.to.save = c("xi", "zi", "fi", "taueps_1", "taueps_2", "lambda_b", "lambda_f", "lambda_w"),
  		                    n.chains = 1, n.burnin = 5000, n.iter = 6000, DIC = TRUE)
  	}
  	else
  	{
  		stop("Please install JAGS")
  	}
  	mortality_female_coda = bugs_chose$BUGSoutput$sims.list
  		
      taueps_female = mortality_female_coda$taueps_1
  	taueps_male   = mortality_female_coda$taueps_2
  		
      lambda_coda_common = mortality_female_coda$lambda_b
  	lambda_coda_male   = mortality_female_coda$lambda_f
  	lambda_coda_female = mortality_female_coda$lambda_w
  
  	score_common = array(NA,c(MCMCiter,ncol(mort_female),ncomp_1))
  	score_female = array(NA,c(MCMCiter,ncol(mort_female),ncomp_female))
  	score_male = array(NA,c(MCMCiter,ncol(mort_male),ncomp_male))
  	
  	if(ncomp_1 == 1)
      {
         for(i in 1:MCMCiter)
         {
             score_common[i,,1] = mortality_female_coda$xi[i,,1]
         }
      }
      else
      {
         for(i in 1:MCMCiter)
         {
             score_common[i,,] = mortality_female_coda$xi[i,,]
         }
  	}
          
      if(ncomp_female == 1)
      {
         for(i in 1:MCMCiter)
         {
             score_female[i,,1] = mortality_female_coda$zi[i,]
         }
      }
      else
      {
         for(i in 1:MCMCiter)
         {
             score_female[i,,] = mortality_female_coda$zi[i,,]
         }
  	}
          
      if(ncomp_male == 1)
      {
         for(i in 1:MCMCiter)
         {
             score_male[i,,1] = mortality_female_coda$fi[i,]
         }
      }
      else
      {
         for(i in 1:MCMCiter)
         {
             score_male[i,,] = mortality_female_coda$fi[i,,]
         }
      }
  	
  	qconf = qnorm(.5 + level/200)
  	score_common_fore = score_common_varfcast = array(, dim = c(MCMCiter, fh, ncomp_1))		
  	for(i in 1:MCMCiter)
  	{
  		if(fmethod == "auto_arima")
  		{
  			for(j in 1:ncomp_1)
  			{
  				dum = forecast(auto.arima(score_common[i,,j]), h = fh, level = level)
  				score_common_varfcast[i,,j] = ((dum$upper - dum$lower)/(2*qconf))^2
  				score_common_fore[i,,j] = dum$mean
  			}
  		}
  		if(fmethod == "ets")
  		{
  			for(j in 1:ncomp_1)
  			{
  				dum = forecast(ets(score_common[i,,j]), h = fh, level = level)
  				score_common_varfcast[i,,j] = ((dum$upper - dum$lower)/(2*qconf))^2
  				score_common_fore[i,,j] = dum$mean
  			}			
  		}
  	}
  
  	score_female_fore = score_female_varfcast = array(,dim=c(MCMCiter, fh, ncomp_female))
  	for(i in 1:MCMCiter)
  	{
  		if(fmethod == "auto_arima")
  		{
  			for(j in 1:ncomp_female)
  			{
  				dum = forecast(auto.arima(score_female[i,,j]), h = fh, level = level)
  				score_female_varfcast[i,,j] = ((dum$upper - dum$lower)/(2*qconf))^2
  				score_female_fore[i,,j] = dum$mean		
  			}
  		}
  		if(fmethod == "ets")
  		{
  			for(j in 1:ncomp_female)
  			{
  				dum = forecast(ets(score_female[i,,j]), h = fh, level = level)
  				score_female_varfcast[i,,j] = ((dum$upper - dum$lower)/(2*qconf))^2
  				score_female_fore[i,,j] = dum$mean		
  			}			
  		}
  	}
  		
  	score_male_fore = score_male_varfcast = array(,dim = c(MCMCiter, fh, ncomp_male))
  	for(i in 1:MCMCiter)
  	{
  		if(fmethod == "auto_arima")
  		{
  			for(j in 1:ncomp_male)
  			{	
  				dum = forecast(auto.arima(score_male[i,,j]), h = fh, level = level)
  				score_male_varfcast[i,,j] = ((dum$upper - dum$lower)/(2*qconf))^2
  				score_male_fore[i,,j] = dum$mean
  			}
  		}
  		if(fmethod == "ets")
  		{
  			for(j in 1:ncomp_male)
  			{	
  				dum = forecast(ets(score_male[i,,j]), h = fh, level = level)
  				score_male_varfcast[i,,j] = ((dum$upper - dum$lower)/(2*qconf))^2
  				score_male_fore[i,,j] = dum$mean
  			}
  		}		
  	}
  	
  	forescore_common = array(, dim = c(MCMCiter, fh, ncomp_1))
  	for(i in 1:MCMCiter)
  	{
  		for(j in 1:ncomp_1)
  		{	
  			for(h in 1:fh)
  			{
  				forescore_common[i,h,j] = rnorm(1, mean = score_common_fore[i,h,j], 
  												sd = sqrt(score_common_varfcast[i,h,j]))
  			}
  		}			
  	}
  			
  	forescore_female = array(, dim = c(MCMCiter, fh, ncomp_female))
  	for(i in 1:MCMCiter)
  	{
  		for(j in 1:ncomp_female)
  		{
  			for(h in 1:fh)
  			{
  				forescore_female[i,h,j] = rnorm(1, mean = score_female_fore[i,h,j], 
  												sd = sqrt(score_female_varfcast[i,h,j]))
  			}
  		}
  	}
  
  	forescore_male = array(, dim = c(MCMCiter, fh, ncomp_male))
  	for(i in 1:MCMCiter)
  	{
  		for(j in 1:ncomp_male)
  		{
  			for(h in 1:fh)
  			{
  				forescore_male[i,h,j] = rnorm(1, mean = score_male_fore[i,h,j],
  												sd = sqrt(score_male_varfcast[i,h,j]))
  			}
  		}
  	}
  	
  	ave_boot_fore = array(, dim = c(MCMCiter, fh, nrow(mort_female)))
  	for(i in 1:MCMCiter)
  	{
  		for(h in 1:fh)
  		{
  			ave_boot_fore[i,h,] = t((forescore_common[i,h,] %*% t(basis_ave)))
  		}
  	}
  			
  	female_boot_fore = array(, dim = c(MCMCiter, fh, nrow(mort_female)))
  	for(i in 1:MCMCiter)
  	{
  		for(h in 1:fh)
  		{
  			female_boot_fore[i,h,] = t((forescore_female[i,h,] %*% t(basis_female))) 
  		}
  	}
  			
  	male_boot_fore = array(, dim = c(MCMCiter, fh, nrow(mort_female)))
  	for(i in 1:MCMCiter)
  	{
  		for(h in 1:fh)
  		{
  			male_boot_fore[i,h,] = t((forescore_male[i,h,] %*% t(basis_male)))
  		}
  	}	
  	
  	female_fore = array(, dim = c(MCMCiter, fh, nrow(mort_female)))
  	for(i in 1:MCMCiter)
  	{
  		for(h in 1:fh)
  		{
  			female_fore[i,h,] = mort_femalemean + ave_boot_fore[i,h,] + female_boot_fore[i,h,]
  		}
  	}
  
  	male_fore = array(, dim = c(MCMCiter, fh, nrow(mort_male)))
  	for(i in 1:MCMCiter)
  	{
  		for(h in 1:fh)
  		{
  			male_fore[i,h,] = mort_malemean + ave_boot_fore[i,h,] + male_boot_fore[i,h,]
  		}
  	}	
  	
  	female_fore_variance = array(, dim = c(MCMCiter, fh, nrow(mort_female)))
  	for(i in 1:MCMCiter)
  	{
  		for(h in 1:fh)
  		{
  			for(j in 1:nrow(mort_female))
  			{
  				female_fore_variance[i,h,j] = rnorm(1, mean = female_fore[i,h,j], 
  													sd = sqrt(1/taueps_female[i]))			
  			}
  		}	
  	}
  			
  	male_fore_variance = array(, dim = c(MCMCiter, fh, nrow(mort_male)))
  	for(i in 1:MCMCiter)
  	{
  		for(h in 1:fh)
  		{
  			for(j in 1:nrow(mort_male))
  			{
  				male_fore_variance[i,h,j] = rnorm(1, mean = male_fore[i,h,j],
  													sd = sqrt(1/taueps_male[i]))		
  			}
  		}
  	}
  	
  	if(BC == TRUE)
  	{
  		mort_female_fore = InvBoxCox(female_fore_variance, lambda)
  		mort_male_fore   = InvBoxCox(male_fore_variance, lambda)
  	}
  	else
  	{
  		mort_female_fore = female_fore_variance
  		mort_male_fore   = male_fore_variance
  	}
  	return(list(first_percent = first_percent, female_percent = female_percent, male_percent = male_percent,
                      mort_female_fore = mort_female_fore, mort_male_fore = mort_male_fore))
}