R/familyStructure2.R

In FamEvent: Family Age-at-Onset Data Simulation and Penetrance Estimation

# only if data with pedigrees and covariates is provided 
familyStructure2 <-
  function (i, cumind=cumind, depend=1, 
            base.dist="Weibull", base.parms=c(0.016, 3),
            data, corr_type = "kinship", IBD = NULL, 
            var_names = NULL, vbeta=c(-1.13, 2.35, 0.5, 0.5), agemin=20)
{
    
    #to samo, ale od razu dla calej rodziny
    tmpdata <- numeric()
    data <- data[order(data$indID),]
    IBD <- IBD[order(data$indID), order(data$indID)]
    
    famID <- data$famID
    indID <- data$indID
    motherID <- data$motherID#c(0,0) 
    fatherID <- data$fatherID#c(0,0) 
    gender <- data$gender  # 0--Female, 1--Male
    proband <- data$proband  
#    proband <- rep(0,dim(data)[1])#c(0,0) ## 0--Non-proband,  1-- proband
    generation <- data$generation#c(1,1) ## generation number, 0 in second generation means married from outside

      nIndi <- length(indID)     
      famID <- data$famID#rep(i, nIndi)    
      ageonset <- rep(0, nIndi)      ## age at onset
      censorage <- rep(0, nIndi)     ## current age
      status <- rep(0, nIndi)        ## disease staus
      affected <- rep(0, nIndi)	
      
      pos <- c(1:nIndi) 
      if(is.null(data$proband)){ # Randomly choose a family member to be a proband
        if(length(which(data$generation==0))>0) potential_probands <- which(data$generation!=0)
        else potential_probands <- 1:nIndi
        pos_proband <- sample(potential_probands, 1)  # Randomly choose a family member to be a proband
        proband[pos_proband] <- 1
        data$proband <- proband
      }
      #------------------------------------------------------------------------------
      # Familial correlation will be generated by gamma frailty, log-normal frailty, or secondgene
      

      ## Log-normal frailty (logX~Normal with mean=0, variance=1/depend)
      if(length(corr_type)==1){
        if(corr_type == "kinship"){
          # different coding of gender in kinship
          gender2 <- ifelse(gender == 1, 1, 2)
          fam <- data.frame(cbind(indID, fatherID, motherID, gender2, famID))
          mped <- with(fam, pedigree(indID, fatherID, motherID, gender2, famid = famID))
          kmat <- as.matrix(kinship(mped))
          alpha <- mvrnorm(1, mu = rep(0, nIndi), Sigma = (1/depend) * 2*kmat)#
        }else if(corr_type == "IBD") alpha <- mvrnorm(1, mu = rep(0, nIndi), Sigma = (1/depend) * IBD) #diag(1, nIndi))#
      }
      else{ #corr_type = c("kinship", "IBD")
          gender2 <- ifelse(gender == 1, 1, 2)
          fam <- data.frame(cbind(indID, fatherID, motherID, gender2, famID))
          mped <- with(fam, pedigree(indID, fatherID, motherID, gender2, famid = famID))
          kmat <- kinship(mped)
          alpha <- mvrnorm(1, mu = rep(0, nIndi), Sigma = 1/depend[1]*2*kmat + 1/depend[2]* IBD)
        }
      
      
      # Individual random effects correlated via the IBD matrix
      #  IBD <- read.table("F:/Stat/PostDoc/FamEvent/IBD_matrix.txt", header = F, sep = " ")
      #  IBD <- as.matrix(IBD)
      # alpha <- mvrnorm(1, mu = rep(0, nIndi), Sigma = depend * IBD[(cumind+1):(cumind + nIndi), (cumind+1):(cumind+nIndi)])
      #print(kmat)
      
      #------------------------------------------------------------------------------
      
      xmat <- data[, var_names]
      xvbeta <- c(t(as.matrix(xmat)%*%vbeta))
      uni <- runif(nIndi)
      
      ageonset <- apply(cbind(xvbeta,alpha,uni), 1, inv2.surv, base.dist=base.dist, parms=base.parms)
      ageonset <- ageonset + agemin
      currentage <- data$currentage  
      
      time <- pmin(currentage, ageonset)  
      status <- ifelse(currentage >= ageonset, 1, 0)
      fsize <- length(famID) # family size
      naff <- sum(status) # number of affected in the family
      tmpdata <- cbind(data, ageonset, time, status, fsize, naff)
      
    return(tmpdata)
    
  }