## Supplemental file S1 for article:
## Title: "Association studies between chromosomal regions 1q21.3, 5q21.3, 14q21.2 and 17q21.31 and numbers of children in Poland."

## R coding for weighted zero-inflated negative binomial regression, with non-events set as 0 to 3 children.
## INSTRUCTIONS: set choices between lines below, then run

ipak <- function(pkg){ 
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE, repos="http://cran.r-project.org")
sapply(pkg, require, character.only = TRUE)
} ## end function
packages <- c("boot", "Exact", "glmmTMB", "nloptr", "data.table", "survival", "MASS", "rcompanion", "Matrix", "performance", "lme4", "pscl", "car", "dplyr", "stringr") 
ipak(packages) ## References at the end of this file

## _____________________________________________________
DDD_FILE <- "" 			## path + chosen .csv file here (wihtout headers)  
nickname <- "NegativeSNPs" ## "H1_H2" or "NegativeSNPs" or "Barban"
sexOnly <- "" ## "MALE" or "FEMALE" or "" for both
## possible override selecting some SNPs by column number:
SNPoverride <- FALSE  ## TRUE or FALSE
SNPoverrideVec <- NA		## NA or vector of column numbers
## _____________________________________________________

if (sexOnly == "MALE") {
	sexOnly2 <- "Male_Only"
}
if (sexOnly == "FEMALE") {
	sexOnly2 <- "Female_Only"
}
if (sexOnly == "") {
	sexOnly2 <- "Male + Female"
}

WORKING_FOLDER<-dirname(DDD_FILE)
WORKING_FOLDER
setwd(WORKING_FOLDER)
BASE_NAME<-basename(WORKING_FOLDER)
NAME_FILE<-BASE_NAME
NAME_FILE
PHENO_GENO<-read.csv(DDD_FILE)

PHENO_GENO <- fread(DDD_FILE, encoding = "Latin-1")
PHENO_GENO <- as.data.frame(PHENO_GENO)

head(PHENO_GENO)

PHENO_GENO_SEL<-PHENO_GENO
PHENO_GENO_SEL<-PHENO_GENO_SEL[PHENO_GENO_SEL$Age!=-9,]
PHENO_GENO_SEL<-PHENO_GENO_SEL[PHENO_GENO_SEL$NCI!=-9,]

PHENO_GENO_SEL2<-PHENO_GENO_SEL[(PHENO_GENO_SEL$Age>44 & PHENO_GENO_SEL$SEX==1) | (PHENO_GENO_SEL$Age>54 & PHENO_GENO_SEL$SEX==2),]

if (sexOnly == "FEMALE") { ## override
PHENO_GENO_SEL2<-PHENO_GENO_SEL[(PHENO_GENO_SEL$Age>44 & PHENO_GENO_SEL$SEX==1), ]
} ## from if (sexOnly == 

if (sexOnly == "MALE") { ## override
PHENO_GENO_SEL2<-PHENO_GENO_SEL[(PHENO_GENO_SEL$Age>54 & PHENO_GENO_SEL$SEX==2), ]
} ## from if (sexOnly == 

PHENO_GENO_SEL1 <- PHENO_GENO_SEL2

## Setting non-events:
for(n in 1:nrow(PHENO_GENO_SEL1)){
  if(PHENO_GENO_SEL1$NCI[n]<=3){
    PHENO_GENO_SEL1$NCI[n]<-0
  }
}

## creation of weights using age bins. 
Age <- PHENO_GENO_SEL1$Age
tableAge <- table(Age)
numSubjects <- length(Age)
wts_table <- 1/tableAge
## scale so that sum = number of subjects:
## create wts_scaled vector:
wts2 <- list(); 
for (i in 1:nrow(PHENO_GENO_SEL1)) {	
	wts2[[i]] <- wts_table[as.character(PHENO_GENO_SEL1[i , "Age"])]
	}
nameswts2 <- list(); 
for (i in 1:length(wts2)) {
nameswts2[[i]] <- names(wts2[[i]])
}
wts <- as.vector(unlist(wts2))
names(wts) <- nameswts2
length(wts)
wts_scaled <- (wts * numSubjects) / sum(wts)
sum(wts_scaled)
PHENO_GENO_SEL1$wts_scaled <- wts_scaled

## "MALE" creation of weights using age bins. 
if (sexOnly == "MALE") { ## override
Age <- PHENO_GENO_SEL1$Age
tableAge <- table(Age)
numSubjects <- length(Age)
wts_table <- 1/tableAge
## scale so that sum = number of subjects:
## create wts_scaled vector:
wts2 <- list(); 
for (i in 1:nrow(PHENO_GENO_SEL1)) {	
	wts2[[i]] <- wts_table[as.character(PHENO_GENO_SEL1[i , "Age"])]
	}
nameswts2 <- list(); 
for (i in 1:length(wts2)) {
nameswts2[[i]] <- names(wts2[[i]])
}
wts <- as.vector(unlist(wts2))
names(wts) <- nameswts2
length(wts)
wts_scaled <- (wts * numSubjects) / sum(wts)
sum(wts_scaled)
PHENO_GENO_SEL1$wts_scaled <- wts_scaled
} ## if (sexOnly == "MALE")

## "FEMALE" creation of weights using age bins. 
if (sexOnly == "FEMALE") { ## override
Age <- PHENO_GENO_SEL1$Age
tableAge <- table(Age)
numSubjects <- length(Age)
wts_table <- 1/tableAge
## scale so that sum = number of subjects:
## create wts_scaled vector:
wts2 <- list(); 
for (i in 1:nrow(PHENO_GENO_SEL1)) {	
	wts2[[i]] <- wts_table[as.character(PHENO_GENO_SEL1[i , "Age"])]
	}
nameswts2 <- list(); 
for (i in 1:length(wts2)) {
nameswts2[[i]] <- names(wts2[[i]])
}
wts <- as.vector(unlist(wts2))
names(wts) <- nameswts2
length(wts)
wts_scaled <- (wts * numSubjects) / sum(wts)
sum(wts_scaled)
PHENO_GENO_SEL1$wts_scaled <- wts_scaled
} ## if (sexOnly == "FEMALE")

## rs12373123 = representative SNP - number 42
if (nickname == "H1_H2") {
PHENO_GENO_SEL2[ , "rs12373123"]
PHENO_GENO_SEL2[ , "NCI"]
} ## from if (nickname == "H1_H2")

if (nickname == "H1_H2") {
if (sexOnly == "") {
mydf <- data.frame(PHENO_GENO_SEL2[ , "rs12373123"], PHENO_GENO_SEL2[ , "NCI"])
mytable2 <- table(mydf)
mytable1 <- as.matrix(mytable2, nrow = 3)
mytable1 <- rbind(mytable1[1, ], mytable1[2, ], mytable1[3, ]) 
mytable <- as.data.frame(mytable1)
rownames(mytable) <- rownames(mytable2)
HetplusHomoH1 <- mytable[2, ] + mytable[3, ]
HetplusHomoH2 <- mytable[1, ] + mytable[2, ]
Allele1st <- (mytable[1, ] * 2) + mytable[2, ]
Allele2nd <- mytable[2, ] + (mytable[3, ] * 2)
mytable <- rbind(mytable, HetplusHomoH1, HetplusHomoH2, Allele1st, Allele2nd)
rownames(mytable)[4] <- paste0(rownames(mytable)[2], "+", rownames(mytable)[3])
rownames(mytable)[5] <- paste0(rownames(mytable)[1], "+", rownames(mytable)[2])
rownames(mytable)[6] <- paste0("allele_", substr(rownames(mytable)[1],1,1))
rownames(mytable)[7] <-  paste0("allele_", substr(rownames(mytable)[3],1,1))
mytable$from0to13 <- c(sum(mytable[1, ]), sum(mytable[2, ]), sum(mytable[3, ]), sum(mytable[4, ]), sum(mytable[5, ]), sum(mytable[6, ]), sum(mytable[7, ]))
mytable$mean <- NA
colnameNums <- as.numeric(c(colnames(mytable)[1:14]))
mytable[1, "mean"] <- signif(sum(colnameNums * mytable[1, c(1:14)]) / mytable[1, "from0to13"], 4)
mytable[2, "mean"] <- signif(sum(colnameNums * mytable[2, c(1:14)]) / mytable[2, "from0to13"], 4)
mytable[3, "mean"] <- signif(sum(colnameNums * mytable[3, c(1:14)]) / mytable[3, "from0to13"], 4)
mytable[4, "mean"] <- signif(sum(colnameNums * mytable[4, c(1:14)]) / mytable[4, "from0to13"], 4)
mytable[5, "mean"] <- signif(sum(colnameNums * mytable[5, c(1:14)]) / mytable[5, "from0to13"], 4)
mytable[6, "mean"] <- signif(sum(colnameNums * mytable[6, c(1:14)]) / mytable[6, "from0to13"], 4)
mytable[7, "mean"] <- signif(sum(colnameNums * mytable[7, c(1:14)]) / mytable[7, "from0to13"], 4)
mytable
} ## from if (sexOnly == "") 
} ## from if (nickname == "H1_H2

if (nickname == "H1_H2") {
if (sexOnly == "" | sexOnly == "FEMALE") {
mydf_F <- data.frame(PHENO_GENO_SEL2[PHENO_GENO_SEL2$SEX %in% 1, "rs12373123"], PHENO_GENO_SEL2[PHENO_GENO_SEL2$SEX %in% 1, "NCI"])
mytable2_F <- table(mydf_F)
mytable1_F <- as.matrix(mytable2_F, nrow = 3)
mytable1_F <- rbind(mytable1_F[1, ], mytable1_F[2, ], mytable1_F[3, ]) 
mytable_F <- as.data.frame(mytable1_F)
## add missing columns - ONLY for females !
mytable_F["11"] <- c(0, 0, 0)
mytable_F["12"] <- c(0, 0, 0)
mytable_F["13"] <- c(0, 0, 0)
rownames(mytable_F) <- rownames(mytable2_F)
HetplusHomoH1 <- mytable_F[2, ] + mytable_F[3, ]
HetplusHomoH2 <- mytable_F[1, ] + mytable_F[2, ]
Allele1st <- (mytable_F[1, ] * 2) + mytable_F[2, ]
Allele2nd <- mytable_F[2, ] + (mytable_F[3, ] * 2)
mytable_F <- rbind(mytable_F, HetplusHomoH1, HetplusHomoH2, Allele1st, Allele2nd)
rownames(mytable_F)[4] <- paste0(rownames(mytable_F)[2], "+", rownames(mytable_F)[3])
rownames(mytable_F)[5] <- paste0(rownames(mytable_F)[1], "+", rownames(mytable_F)[2])
rownames(mytable_F)[6] <- paste0("allele_", substr(rownames(mytable_F)[1],1,1))
rownames(mytable_F)[7] <-  paste0("allele_", substr(rownames(mytable_F)[3],1,1))
mytable_F$from0to13 <- c(sum(mytable_F[1, ]), sum(mytable_F[2, ]), sum(mytable_F[3, ]), sum(mytable_F[4, ]), sum(mytable_F[5, ]), sum(mytable_F[6, ]), sum(mytable_F[7, ]))
mytable_F$mean <- NA
colnameNums <- as.numeric(c(colnames(mytable_F)[1:14]))
mytable_F[1, "mean"] <- signif(sum(colnameNums * mytable_F[1, c(1:14)]) / mytable_F[1, "from0to13"], 4)
mytable_F[2, "mean"] <- signif(sum(colnameNums * mytable_F[2, c(1:14)]) / mytable_F[2, "from0to13"], 4)
mytable_F[3, "mean"] <- signif(sum(colnameNums * mytable_F[3, c(1:14)]) / mytable_F[3, "from0to13"], 4)
mytable_F[4, "mean"] <- signif(sum(colnameNums * mytable_F[4, c(1:14)]) / mytable_F[4, "from0to13"], 4)
mytable_F[5, "mean"] <- signif(sum(colnameNums * mytable_F[5, c(1:14)]) / mytable_F[5, "from0to13"], 4)
mytable_F[6, "mean"] <- signif(sum(colnameNums * mytable_F[6, c(1:14)]) / mytable_F[6, "from0to13"], 4)
mytable_F[7, "mean"] <- signif(sum(colnameNums * mytable_F[7, c(1:14)]) / mytable_F[7, "from0to13"], 4)
mytable_F
} ## from if (sexOnly == "" | sexOn
} ## from if (nickname == "H1_H2") 

if (nickname == "H1_H2") {
if (sexOnly == "" | sexOnly == "MALE") {
mydf_M <- data.frame(PHENO_GENO_SEL2[PHENO_GENO_SEL2$SEX %in% 2, "rs12373123"], PHENO_GENO_SEL2[PHENO_GENO_SEL2$SEX %in% 2, "NCI"])
mytable2_M <- table(mydf_M)
mytable1_M <- as.matrix(mytable2_M, nrow = 3)
mytable1_M <- rbind(mytable1_M[1, ], mytable1_M[2, ], mytable1_M[3, ]) 
mytable_M <- as.data.frame(mytable1_M)
rownames(mytable_M) <- rownames(mytable2_M)
HetplusHomoH1 <- mytable_M[2, ] + mytable_M[3, ]
HetplusHomoH2 <- mytable_M[1, ] + mytable_M[2, ]
Allele1st <- (mytable_M[1, ] * 2) + mytable_M[2, ]
Allele2nd <- mytable_M[2, ] + (mytable_M[3, ] * 2)
mytable_M <- rbind(mytable_M, HetplusHomoH1, HetplusHomoH2, Allele1st, Allele2nd)
rownames(mytable_M)[4] <- paste0(rownames(mytable_M)[2], "+", rownames(mytable_M)[3])
rownames(mytable_M)[5] <- paste0(rownames(mytable_M)[1], "+", rownames(mytable_M)[2])
rownames(mytable_M)[6] <- paste0("allele_", substr(rownames(mytable_M)[1],1,1))
rownames(mytable_M)[7] <-  paste0("allele_", substr(rownames(mytable_M)[3],1,1))
mytable_M$from0to13 <- c(sum(mytable_M[1, ]), sum(mytable_M[2, ]), sum(mytable_M[3, ]), sum(mytable_M[4, ]), sum(mytable_M[5, ]), sum(mytable_M[6, ]), sum(mytable_M[7, ]))
mytable_M$mean <- NA
colnameNums <- as.numeric(c(colnames(mytable_M)[1:14]))
mytable_M[1, "mean"] <- signif(sum(colnameNums * mytable_M[1, c(1:14)]) / mytable_M[1, "from0to13"], 4)
mytable_M[2, "mean"] <- signif(sum(colnameNums * mytable_M[2, c(1:14)]) / mytable_M[2, "from0to13"], 4)
mytable_M[3, "mean"] <- signif(sum(colnameNums * mytable_M[3, c(1:14)]) / mytable_M[3, "from0to13"], 4)
mytable_M[4, "mean"] <- signif(sum(colnameNums * mytable_M[4, c(1:14)]) / mytable_M[4, "from0to13"], 4)
mytable_M[5, "mean"] <- signif(sum(colnameNums * mytable_M[5, c(1:14)]) / mytable_M[5, "from0to13"], 4)
mytable_M[6, "mean"] <- signif(sum(colnameNums * mytable_M[6, c(1:14)]) / mytable_M[6, "from0to13"], 4)
mytable_M[7, "mean"] <- signif(sum(colnameNums * mytable_M[7, c(1:14)]) / mytable_M[7, "from0to13"], 4)
mytable_M
} ## from if (sexOnly == "" | sexOn
} ## from if (nickname == "H1_H2") 

INTERACTION_TABELKA<-as.data.frame(matrix(nrow=1, ncol=2))
colnames(INTERACTION_TABELKA)<-c("SNP","ZERO_INFLA")
INTERACTION_TABELKA
PROBLEM_SNP<-c()
				
TABLIST <- list(); TABLIST_INTERACTION <- list(); WEK_COUNT <- list(); INTERACTION_COUNT <- list(); INTERACTION_COUNT_1 <- list(); INTERACTION_COUNT_2 <- list(); TAB_COUNT <- list(); genotypelowp <- list(); genotype_1lowp <- list(); genotype_2lowp <- list(); TAB_PART1_no_interaction <- list(); TAB_PART2_interaction <- list(); TAB_ANOVA  <- list(); TAB_ANOVA_INT  <- list(); ANOVA_INT  <- list(); TAB_INTERACTION <- list(); 
if (sexOnly == "") {
TABELKA<-as.data.frame(matrix(nrow=8, ncol = 5))
}
if ((sexOnly == "MALE") || (sexOnly == "FEMALE")) {
TABELKA<-as.data.frame(matrix(nrow=6, ncol = 5))
}

colnames(TABELKA)<-c("Estimate", "Std. Error or df", "p-value", "low p-value ?", "Interaction with age ?")            	
## 14:ncol(PHENO_GENO_SEL1)
whichSNPs <- c(14:(ncol(PHENO_GENO_SEL1)-1))

if (SNPoverride == TRUE) {
whichSNPs <- SNPoverrideVec
}

## adjusted with interaction:
for (n in c(whichSNPs)) {

problem <- "no" ## overridden later

SNP<-colnames(PHENO_GENO_SEL1)[n]
  print(n)    
 if(SNP == "rs4792887" || ##  H1_H2 can't form model 
SNP == "rs75638861" || ## H1_H2 only two genotypes 
SNP == "rs144845296" || ## H1_H2 can't form model 
SNP == "rs63750072" || ## H1_H2 can't form model 
SNP == "rs2732674"|| ## H1_H2 only two genotypes 
SNP == "rs12947457" || ## H1_H2 can't form model

SNP == "rs11012" || ## H1_H2 can't form interaction model 
SNP == "rs17631676" || ## H1_H2 can't form interaction model 
SNP == "rs17763658" || ## H1_H2 can't form interaction model
SNP == "rs242928" || ## H1_H2 can't form interaction model
SNP == "rs2532419" || ## H1_H2 can't form interaction model
SNP == "rs1378358" || ## H1_H2 can't form interaction model
SNP == "rs199533" ) ## H1_H2 can't form interaction model
{
problem <- "yes"
 PROBLEM_SNP<-c(PROBLEM_SNP, SNP)
 }
 
if (sexOnly == "MALE") {
	if(SNP == "rs12373123" ) ## MALE - H1_H2 representatibve SNP rs12373123 can't form model
{
problem <- "yes"
PROBLEM_SNP<-c(PROBLEM_SNP, SNP)
 }
 } ## from if (sexOnly == "MALE")

if (sexOnly == "FEMALE") {
	if(SNP == "rs12373123" || ## FEMALE - H1_H2 representative SNP rs12373123 can't form model
SNP == "rs10505385") ## FEMALE negative SNPs 
{
problem <- "yes"
PROBLEM_SNP<-c(PROBLEM_SNP, SNP)
 }
 } ## from if (sexOnly == "FEMALE")

if (SNP == "rs354082" || ## Negative - can't form model
SNP == "rs7358238" || ##  Negative - can't form model
SNP == "rs418889" || ##  Negative - too many NAs
SNP == "rs6121444" ||  ##  Negative - only two genotypes ?
SNP == "rs113831133" || ##  Negative - only two genotypes 

SNP == "rs10068788" || ##  Negative - can't form interaction model
SNP == "rs145522270" || ##  Negative - can't form model
SNP == "rs2242312" || ##  Negative - can't form interaction model
SNP == "rs9939343" || ##  Negative - can't form interaction model
SNP == "rs1507732" || ##  Negative - can't form interaction model
SNP == "rs619865" ) ##  Negative - can't form interaction model
{
problem <- "yes"
PROBLEM_SNP<-c(PROBLEM_SNP, SNP)
 }
 
if (problem == "no")  {
    p<-na.omit(PHENO_GENO_SEL1[ , c(2, 5, 10, 8, n, ncol(PHENO_GENO_SEL1))])
    p[,5]<-as.factor(p[,5])

myweights <- p[ , 6]

if (sexOnly == "") {
MODEL <- zeroinfl(NCI ~ p[,5] + Age + SEX, data = p, dist = "negbin", weights = myweights)
MODEL_INT <- zeroinfl(NCI ~ p[,5] * Age * SEX, data = p, dist = "negbin", weights = myweights)
} ## from if (sexOnly == ""
if ((sexOnly == "MALE") || (sexOnly == "FEMALE")) {
MODEL <- zeroinfl(NCI ~ p[,5] + Age, data = p, dist = "negbin", weights = myweights)
MODEL_INT <- zeroinfl(NCI ~ p[,5] * Age, data = p, dist = "negbin", weights = myweights)
} ## from if (sexOnly == "MALE"

TAB_INTERACTION[[n]] <- summary(MODEL_INT)
TABLIST_INTERACTION[[n]] <- TAB_INTERACTION[[n]]
names(TABLIST_INTERACTION)[n] <- SNP
  
TAB_COUNT[[n]]<-as.data.frame(TAB_INTERACTION[[n]]$coefficients$count)
 
print(SNP)

               if (sexOnly == "") {
            WEK_COUNT[[n]]<-TAB_COUNT[[n]][-(1:5), 4]<=0.05                
      if(is.na(WEK_COUNT[[n]][1])==TRUE){
        PROBLEM_SNP<-c(PROBLEM_SNP, SNP)
      } else {
        if(length(which(WEK_COUNT[[n]]==FALSE)) == 7){ ## No interactions
          INTERACTION_COUNT[[n]]<-c("No")
        }else{
          INTERACTION_COUNT[[n]]<-c("Yes")
          }
          }
          } ## from  if (sexOnly == ""
          
               if ((sexOnly == "MALE") || (sexOnly == "FEMALE")) {
            WEK_COUNT[[n]]<-TAB_COUNT[[n]][-(1:4), 4]<=0.05
                  if(is.na(WEK_COUNT[[n]][1])==TRUE){
        PROBLEM_SNP<-c(PROBLEM_SNP, SNP)
      } else {
        if(length(which(WEK_COUNT[[n]]==FALSE)) == 2){ ## No interactions
          INTERACTION_COUNT[[n]]<-c("No")
        }else{
          INTERACTION_COUNT[[n]]<-c("Yes")
          }
          }
          } ## from  if (sexOnly ==
            
          ## the following only considers first two interactions: with Age:
        if (length(WEK_COUNT[[n]][1]) != 0) {
        	if (!is.na(WEK_COUNT[[n]][1])) {
        if(WEK_COUNT[[n]][1] ==FALSE) { 
          INTERACTION_COUNT_1[[n]]<-c("No")
        } else {
          INTERACTION_COUNT_1[[n]]<-c("Yes")
          }
        } else {
        	 INTERACTION_COUNT_1[[n]] <- NA
        	} ## fron if (!is.na(W
        	} else {
          	 INTERACTION_COUNT_1[[n]] <- NA
        	} ## from if (length(W_ 		
        if (length(WEK_COUNT[[n]][2]) != 0) {
        	if (!is.na(WEK_COUNT[[n]][2])) {
        	  if(WEK_COUNT[[n]][2] ==FALSE) {
          INTERACTION_COUNT_2[[n]]<-c("No")
        }else{
          INTERACTION_COUNT_2[[n]]<-c("Yes")
          }
          } else {
        	 INTERACTION_COUNT_2[[n]] <- NA
        	} ## fron if (!is.na(W
        	} else {
          	 INTERACTION_COUNT_2[[n]] <- NA
        	} ## from if (length(W_ 
        	
 ## Not adjusted with interaction; Anova used identically as in Mangiafico, S.S. (2015):
 TAB_PART1_no_interaction[[n]] <- Anova(MODEL, type="II", test="Chisq")
TAB_ANOVA[[n]] <- as.data.frame(as.matrix(TAB_PART1_no_interaction[[n]]))
TAB_ANOVA[[n]] <-  TAB_ANOVA[[n]] %>% relocate("Df", .after = last_col())
TAB_ANOVA[[n]] <-  TAB_ANOVA[[n]] %>% relocate("Pr(>Chisq)", .after = last_col())
rownames( TAB_ANOVA[[n]]) <-str_replace(rownames(TAB_ANOVA[[n]]), "p\\[, 5\\]", SNP)
if (as.numeric(as.character(TAB_ANOVA[[n]][1, 3])) < 0.05) {
	genotypelowp[[n]] <- "yes"
} else {
	genotypelowp[[n]] <- "No"
	}
names(TAB_PART1_no_interaction)[n] <- SNP

## Adjusted with interaction:

TAB_PART2_interaction[[n]] <- Anova(MODEL_INT, type="II", test="Chisq")
 TAB_ANOVA_INT[[n]] <- as.data.frame(as.matrix(TAB_PART2_interaction[[n]]))
 TAB_ANOVA_INT[[n]] <-  TAB_ANOVA_INT[[n]] %>% relocate("Df", .after = last_col())
 TAB_ANOVA_INT[[n]] <-  TAB_ANOVA_INT[[n]] %>% relocate("Pr(>Chisq)", .after = last_col())
 
 
if (sexOnly == "") {
if (TAB_ANOVA_INT[[n]][4, 3] < 0.05) { ## only considers p[, 5]:Age
	ANOVA_INT[[n]] <- "Yes"
} else {
	ANOVA_INT[[n]] <- "No"
}
} ## from if (sexOnly == "") 

if ((sexOnly == "MALE") || (sexOnly == "FEMALE")) {
if (TAB_ANOVA_INT[[n]][3, 3] < 0.05) { 
	ANOVA_INT[[n]] <- "Yes"
} else {
	ANOVA_INT[[n]] <- "No"
}
} ## from if (sexOnly == 

names(TAB_PART2_interaction)[[n]] <- SNP

TAB_FULL <- summary(MODEL)
TABLIST[[n]] <- TAB_FULL
names(TABLIST)[n] <- SNP
          TAB_COUNT[[n]]<-as.data.frame(TAB_FULL$coefficients$count)
          TAB_COUNT[[n]]<- TAB_COUNT[[n]][-1, -3]
if (as.numeric(as.character(TAB_COUNT[[n]][1, 3])) < 0.05) {
	genotype_1lowp[[n]] <- "yes"
} else {
	genotype_1lowp[[n]] <- "No"
	}
if (as.numeric(as.character(TAB_COUNT[[n]][2, 3]) < 0.05)) {
	genotype_2lowp[[n]] <- "yes"
} else {
	genotype_2lowp[[n]] <- "No"
	}

            if (sexOnly == "") {
            TEMP_TABELKA<-as.data.frame(matrix(nrow = 8, ncol = 5))            
            colnames(TEMP_TABELKA) <- c("Estimate", "Std. Error or df", "p-value", "low p-value ?", "Interaction with age ?")   
            rownames(TEMP_TABELKA) <- c(paste(SNP, rownames(TAB_COUNT[[n]])[1], sep = ":"), paste(SNP, rownames(TAB_COUNT[[n]])[2], sep = ":"),  paste(SNP,rownames(TAB_COUNT[[n]])[3], sep = ":"), paste(SNP,rownames(TAB_COUNT[[n]])[4], sep = ":"), paste(SNP,rownames(TAB_COUNT[[n]])[5], sep = ":"),        
paste(SNP, "genotype", sep = ":"), paste(SNP, "ANOVA", rownames(TAB_ANOVA[[n]])[2], sep = ":"), paste(SNP, "ANOVA", rownames(TAB_ANOVA[[n]])[3], sep = ":")) 
            TEMP_TABELKA[1:5, 1:3] <- TAB_COUNT[[n]]
            TEMP_TABELKA[6:8, ] <- TAB_ANOVA[[n]]
     } ## from   if (sexOnly == "")
     
 if ((sexOnly == "MALE") || (sexOnly == "FEMALE")) {
            TEMP_TABELKA<-as.data.frame(matrix(nrow = 6, ncol = 5))            
            colnames(TEMP_TABELKA) <- c("Estimate", "Std. Error or df", "p-value", "low p-value ?", "Interaction with age ?")   
           rownames(TEMP_TABELKA) <- c(paste(SNP,rownames(TAB_COUNT[[n]])[1], sep = ":"), paste(SNP,rownames(TAB_COUNT[[n]])[2], sep = ":"),  paste(SNP,rownames(TAB_COUNT[[n]])[3], sep = ":"), paste(SNP,rownames(TAB_COUNT[[n]])[4], sep = ":"), paste(SNP, "genotype", sep = ":"), paste(SNP, "ANOVA", rownames(TAB_ANOVA[[n]])[2], sep = ":")) 
           
            TEMP_TABELKA[1:4 , 1:3] <- TAB_COUNT[[n]]
            TEMP_TABELKA[5:6, ] <- TAB_ANOVA[[n]]
                 } ## from   if (sexOnly == 
     
 if (sexOnly == "") {
 if ((length(genotype_1lowp) != 0) & (length(genotype_2lowp) != 0)) {
TEMP_TABELKA[ , 4] <- c(genotype_1lowp[[n]], genotype_2lowp[[n]], "", "", "", genotypelowp[[n]], "", "")   
}
} ## from  if (sexOnly == "")

if ((sexOnly == "MALE") || (sexOnly == "FEMALE")) {
 if ((length(genotype_1lowp) != 0) & (length(genotype_2lowp) != 0)) {
TEMP_TABELKA[ , 4] <- c(genotype_1lowp[[n]], genotype_2lowp[[n]], "", "", genotypelowp[[n]], "")   
}
} ## from if ((sexOnly == 

if (length(INTERACTION_COUNT_1[[n]]) != 0) {
if (!is.na(INTERACTION_COUNT_1[[n]])) {
     	if (INTERACTION_COUNT_1[[n]] == "Yes") {
			TEMP_TABELKA[1, 5] <- c("Yes")
		} else {
			TEMP_TABELKA[1, 5] <- c("No")
		}
		} else {
			TEMP_TABELKA[1, 5] <- NA			
		} ## from if (!is.na(INTE
		} else {
			TEMP_TABELKA[1, 5] <- NA			
		} ## from if (length(INT

if (length(INTERACTION_COUNT_2[[n]]) != 0) {
if (!is.na(INTERACTION_COUNT_2[[n]])) {
	    if (INTERACTION_COUNT_2[[n]] == "Yes") {
			TEMP_TABELKA[2, 5] <- c("Yes")
		} else {
			TEMP_TABELKA[2, 5] <- c("No")
		}
		} else {
			TEMP_TABELKA[2, 5] <- NA			
		} ## from if (!is.na(INTE
		} else {
			TEMP_TABELKA[2, 5] <- NA			
		} ## from if (length(INT		

if (sexOnly == "") {	
		if (ANOVA_INT[[n]] == "Yes") {
			TEMP_TABELKA[6, 5] <- c("Yes")
		} else {
			TEMP_TABELKA[6, 5] <- c("No")
		}		
		TEMP_TABELKA[c(3:5, 7:8), 5] <- c("")
		} ## from if (sexOnly == "")

if ((sexOnly == "MALE") || (sexOnly == "FEMALE")) {
		if (ANOVA_INT[[n]] == "Yes") {
			TEMP_TABELKA[5, 5] <- c("Yes")
		} else {
			TEMP_TABELKA[5, 5] <- c("No")
		}		
		TEMP_TABELKA[c(3:4, 6), 5] <- c("")
		} ## from if (sexOnly ==	
		
		    TABELKA <- rbind(TABELKA, TEMP_TABELKA)
            TEMP_TABELKA<-c()

TABELKA2 <- TABELKA
TABELKA2$rownames <- rownames(TABELKA2)
TABELKA2 <- TABELKA2 %>% select("rownames", everything())
colnames(TABELKA2)[1] <- sexOnly2

    write.csv(TABELKA2, paste0(WORKING_FOLDER,"/", NAME_FILE, nickname, sexOnly, "_ZERO_INFL.csv"), row.names = FALSE)
    write.csv(INTERACTION_TABELKA, paste0(WORKING_FOLDER,"/", NAME_FILE, nickname, sexOnly, "_INTERACTION.csv"), row.names = FALSE)
} ## from for(n in c(whichSNPs)
} ## from for(n in c(wh

summary(TABELKA)

print(PROBLEM_SNP)

TABLIST_INTERACTION
TABLIST
WEK_COUNT
INTERACTION_COUNT
TAB_COUNT

if (nickname == "H1_H2" & sexOnly == "") {
mytable
}
if (nickname == "H1_H2" & (sexOnly == "" | sexOnly == "FEMALE")) {
mytable_F
}
if (nickname == "H1_H2" & (sexOnly == "" | sexOnly == "MALE")) {
mytable_M
}

TAB_PART1_no_interaction
TAB_PART2_interaction

## References:

## (1) "boot", (2) "Exact", (3) "glmmTMB", (4) "nloptr", (5) "data.table", (6) "survival", (7) "MASS", (8) "rcompanion", (9) "Matrix", (10) "performance", (11) "lme4", (12) "pscl", (13) "car", (14) "dplyr", (15) "stringr" (16) "R companion to the Handbook of Biological Statistics"

## 1.	Achim Z, Christian K, Simon J. Regression Models for Count Data in R. Journal of Statistical Software. 2008;27
## 2.	Angelo C, B. DR. boot: Bootstrap R (S-Plus) Functions. 2021
## 3.	Daniel L, Mattan SB-S, Indrajeet P, Philip W, Dominique M. performance: An R Package for Assessment, Comparison and Testing of Statistical Models. Journal of Open Source Software. 2021;6:3139.
## 4.	Douglas B, Martin M, Ben B, Steve W. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software. 2015;67:1-48.
## 5.	Douglas B, Martin M, Mikael J. Matrix: Sparse and Dense Matrix Classes and Methods. 2022
## 6.	Hadley W. stringr: Simple, Consistent Wrappers for Common String Operations. 2019
## 7.	Hadley W, Romain F, Lionel H, Kirill M. dplyr: A Grammar of Data Manipulation. 2022
## 8.	John F, Sanford W. An R Companion to Applied Regression. Thousand Oaks CA: Sage; 2019
## 9.	Matt D, Arun S. data.table: Extension of `data.frame`. 2021
## 10.	Mollie EB, Kasper K, Koen JVB et al. glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal. 2017;9:378-400.
## 11.	Peter C. Exact: Unconditional Exact Test. 2021
## 12.	Salvatore M. rcompanion: Functions to Support Extension Education Program Evaluation. 2022
## 13.	Steven GJ. The NLopt nonlinear-optimization package. 2022
## 14.	Terry MT, Patricia MG. Modeling Survival Data: Extending the Cox Model. New York: Springer; 2000
## 15.	W. NV, B. DR. Modern Applied Statistics with S. New York: Springer; 2002
## 16. Mangiafico SS. An R Companion for the Handbook of Biological Statistics, version 1.3.2. 2015.  https://rcompanion.org/rcompanion/