From 4f51569cf615fbe6094b21f317c18601a2c52bc3 Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Tue, 15 Oct 2024 17:09:00 +0700
Subject: [PATCH 01/14] concise the code v1

Concise the code part prior to calculating the weight
---
 R/ipwtm.R | 151 ++++++++++++++++++++++++++++++++----------------------
 1 file changed, 89 insertions(+), 62 deletions(-)

diff --git a/R/ipwtm.R b/R/ipwtm.R
index f1e3dc6..0c08f05 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -1,66 +1,93 @@
 ipwtm <- function(
-	exposure,
-	family,
-	link,
-	numerator = NULL,
-	denominator,
-	id,
-	tstart,
-	timevar,
-	type,
-	data,
-	corstr = "ar1",
-	trunc = NULL,
-	...)
-	{
-		#save input
-			tempcall <- match.call()
-		#some basic input checks
-			if (!("exposure" %in% names(tempcall))) stop("No exposure variable specified")
-			if (!("family" %in% names(tempcall)) | ("family" %in% names(tempcall) & !(tempcall$family %in% c("binomial", "survival", "multinomial", "ordinal", "gaussian")))) stop("No valid family specified (\"binomial\", \"survival\", \"multinomial\", \"ordinal\", \"gaussian\")")
-			if (tempcall$family == "binomial") {if(!(tempcall$link %in% c("logit", "probit", "cauchit", "log", "cloglog"))) stop("No valid link function specified for family = binomial (\"logit\", \"probit\", \"cauchit\", \"log\", \"cloglog\")")}
-			if (tempcall$family == "ordinal" ) {if(!(tempcall$link %in% c("logit", "probit", "cauchit", "cloglog"))) stop("No valid link function specified for family = ordinal (\"logit\", \"probit\", \"cauchit\", \"cloglog\")")}
-			if (!("denominator" %in% names(tempcall))) stop("No denominator model specified")
-			if (!is.null(tempcall$numerator) & !is(eval(tempcall$numerator), "formula")) stop("Invalid numerator formula specified")
-			if (!is.null(tempcall$denominator) & !is(eval(tempcall$denominator), "formula")) stop("Invalid denominator formula specified")
-			if (!("id" %in% names(tempcall))) stop("No patient id specified")
-			if (tempcall$family == "survival" & !("tstart" %in% names(tempcall))) stop("No tstart specified, is necessary for family = \"survival\"")
-			if (!("timevar" %in% names(tempcall))) stop("No timevar specified")
-			if (!("type" %in% names(tempcall))) stop("No type specified (\"first\" or \"all\")")
-			if (!(tempcall$type %in% c("first", "all", "cens"))) stop("No type specified (\"first\", \"all\" or \"cens\")")
-			if (tempcall$family %in% c("survival", "multinomial", "ordinal") & tempcall$type == "all") stop(paste("Type \"all\" not yet implemented for family = ", deparse(tempcall$family, width.cutoff = 500), sep = ""))
-			if (tempcall$family %in% c("multinomial", "ordinal", "gaussian") & tempcall$type == "cens") stop(paste("Type \"cens\" not yet implemented for family = ", deparse(tempcall$family, width.cutoff = 500), sep = ""))
-			if (tempcall$family %in% c("gaussian") & tempcall$type == "first") stop(paste("Type \"first\" not implemented for family = ", deparse(tempcall$family, width.cutoff = 500), sep = ""))
-			if (tempcall$family %in% c("gaussian") & !("numerator" %in% names(tempcall))) stop("Numerator necessary for family = \"gaussian\"")
-			if (!("data" %in% names(tempcall))) stop("No data specified")
-			if (!is.null(tempcall$trunc)) {if(tempcall$trunc < 0 | tempcall$trunc > 0.5) stop("Invalid truncation percentage specified (0-0.5)")}
-		#record original order of dataframe so that the output can be returned in the same order
-			order.orig <- 1:nrow(data)
-			order.orig <- order.orig[order(
-				eval(parse(text = paste("data$", deparse(tempcall$id, width.cutoff = 500), sep = ""))),
-				eval(parse(text = paste("data$", deparse(tempcall$timevar, width.cutoff = 500), sep = "")))
-				)] #sort as below
-		#sort dataframe on follow-up time within each individual, necessary for cumulative products below
-			data <- data[order(
-				eval(parse(text = paste("data$", deparse(tempcall$id, width.cutoff = 500), sep = ""))),
-				eval(parse(text = paste("data$", deparse(tempcall$timevar, width.cutoff = 500), sep = "")))
-				),]
-		#make new dataframe for newly computed variables, to prevent variable name conflicts
-			tempdat <- data.frame(
-				id = data[,as.character(tempcall$id)],
-				timevar = data[,as.character(tempcall$timevar)],
-				exposure = data[,as.character(tempcall$exposure)]
-			)
-		#make selection variable, time points up to first switch from lowest value, or all time points
-			if (type %in% c("first", "cens") & (family == "binomial" | family == "survival"))
-				{tempdat$selvar <- do.call("c", lapply(split(tempdat$exposure, tempdat$id),function(x)if (!is.na(match(1, x))) return(c(rep(1,match(1, x)),rep(0,length(x)-match(1, x)))) else return(rep(1,length(x)))))}
-			if (type %in% c("first", "cens") & (family == "multinomial" | family == "ordinal")){
-					z <- unique(tempdat$exposure)[unique(tempdat$exposure) != sort(unique(tempdat$exposure))[1]]
-					min2 <- function(x)ifelse(min(is.na(unique(x))) == 1, NA, min(x, na.rm = TRUE))
-					tempdat$selvar <- do.call("c", lapply(split(tempdat$exposure, tempdat$id),function(x)if (!is.na(min2(match(z, x)))) return(c(rep(1,min2(match(z, x))),rep(0,length(x)-min2(match(z, x))))) else return(rep(1,length(x)))))
-			}
-			if (type == "all")
-				{tempdat$selvar <- rep(1, nrow(tempdat))}
+    exposure, family, link, denominator, id, tstart, timevar, type, data, 
+    numerator = NULL, corstr = "ar1", trunc = NULL, ...
+) {
+  # Save input
+  tempcall <- match.call()
+  
+  # Helper function to check for required arguments
+  check_required <- function(arg, msg) {
+    if (!(arg %in% names(tempcall))) stop(msg)
+  }
+  
+  # Validate required arguments
+  check_required("exposure", "No exposure variable specified")
+  check_required("denominator", "No denominator model specified")
+  check_required("id", "No patient id specified")
+  check_required("timevar", "No timevar specified")
+  check_required("type", "No type specified (\"first\" or \"all\")")
+  check_required("data", "No data specified")
+  
+  # Family-specific checks
+  valid_families <- c("binomial", "survival", "multinomial", "ordinal", "gaussian")
+  valid_links <- list(
+    binomial = c("logit", "probit", "cauchit", "log", "cloglog"),
+    ordinal = c("logit", "probit", "cauchit", "cloglog")
+  )
+  
+  if (!family %in% valid_families) stop("Invalid family specified")
+  if (family %in% names(valid_links) && !(link %in% valid_links[[family]])) {
+    stop(paste("No valid link function specified for family =", family))
+  }
+  
+  # Special handling for survival
+  if (family == "survival") check_required("tstart", "No tstart specified, necessary for family = \"survival\"")
+  
+  # Formula validation
+  if (!is.null(numerator) && !is(eval(numerator), "formula")) stop("Invalid numerator formula specified")
+  if (!is(eval(denominator), "formula")) stop("Invalid denominator formula specified")
+  
+  # Type validation
+  valid_types <- c("first", "all", "cens")
+  if (!(type %in% valid_types)) stop("Invalid type specified (\"first\", \"all\" or \"cens\")")
+  
+  unsupported_combinations <- list(
+    survival = "all",
+    multinomial = c("all", "cens"),
+    ordinal = c("all", "cens"),
+    gaussian = c("first", "cens")
+  )
+  
+  if (family %in% names(unsupported_combinations) && type %in% unsupported_combinations[[family]]) {
+    stop(paste("Type", type, "not yet implemented for family =", family))
+  }
+  
+  if (family == "gaussian" && is.null(numerator)) stop("Numerator necessary for family = \"gaussian\"")
+  
+  # Truncation check
+  if (!is.null(trunc) && (trunc < 0 || trunc > 0.5)) stop("Invalid truncation percentage specified (0-0.5)")
+  
+  # Record original order of dataframe
+  order.orig <- order(data[[deparse(tempcall$id)]], data[[deparse(tempcall$timevar)]])
+  
+  # Sort dataframe by id and timevar
+  data <- data[order.orig, ]
+  # Create a new dataframe to prevent variable name conflicts
+  tempdat <- data.frame(
+    id = data[[deparse(tempcall$id)]],
+    timevar = data[[deparse(tempcall$timevar)]],
+    exposure = data[[deparse(tempcall$exposure)]]
+  )
+  
+  # Define selection variable based on type and family
+  if (type %in% c("first", "cens")) {
+    if (family %in% c("binomial", "survival")) {
+      tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {
+        idx <- match(1, x)
+        if (!is.na(idx)) c(rep(1, idx), rep(0, length(x) - idx)) else rep(1, length(x))
+      }))
+    } else if (family %in% c("multinomial", "ordinal")) {
+      z <- setdiff(unique(tempdat$exposure), sort(unique(tempdat$exposure))[1])
+      min2 <- function(x) ifelse(all(is.na(x)), NA, min(x, na.rm = TRUE))
+      tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {
+        idx <- min2(match(z, x))
+        if (!is.na(idx)) c(rep(1, idx), rep(0, length(x) - idx)) else rep(1, length(x))
+      }))
+    }
+  } else if (type == "all") {
+    tempdat$selvar <- rep(1, nrow(tempdat))
+  }
+  
 		#weights binomial, type "first"
 			if (tempcall$family == "binomial" & tempcall$type %in% c("first", "cens")) {
 				if(tempcall$link == "logit") lf <- binomial(link = logit)

From a6baecd4acb8e7620d2f9c4ca91f1b9d0cd53177 Mon Sep 17 00:00:00 2001
From: Hung Tran Thai <hungtt@oucru.org>
Date: Wed, 16 Oct 2024 15:18:39 +0700
Subject: [PATCH 02/14] concise the input of family

---
 R/ipwpoint.R |  6 +-----
 R/ipwtm.R    | 12 ++----------
 2 files changed, 3 insertions(+), 15 deletions(-)

diff --git a/R/ipwpoint.R b/R/ipwpoint.R
index 34930d2..cbf1f5d 100644
--- a/R/ipwpoint.R
+++ b/R/ipwpoint.R
@@ -27,11 +27,7 @@ ipwpoint <- function(
 			)
 		#weights binomial
 			if (tempcall$family == "binomial") {
-				if(tempcall$link == "logit") lf <- binomial(link = logit)
-				if(tempcall$link == "probit") lf  <- binomial(link = probit)
-				if(tempcall$link == "cauchit") lf  <- binomial(link = cauchit)
-				if(tempcall$link == "log") lf  <- binomial(link = log)
-				if(tempcall$link == "cloglog") lf  <- binomial(link = cloglog)
+			  lf <- binomial(link = tempcall$link)
 				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
 				else {
 					mod1 <- glm(
diff --git a/R/ipwtm.R b/R/ipwtm.R
index 0c08f05..d595a88 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -90,11 +90,7 @@ ipwtm <- function(
   
 		#weights binomial, type "first"
 			if (tempcall$family == "binomial" & tempcall$type %in% c("first", "cens")) {
-				if(tempcall$link == "logit") lf <- binomial(link = logit)
-				if(tempcall$link == "probit") lf  <- binomial(link = probit)
-				if(tempcall$link == "cauchit") lf  <- binomial(link = cauchit)
-				if(tempcall$link == "log") lf  <- binomial(link = log)
-				if(tempcall$link == "cloglog") lf  <- binomial(link = cloglog)
+			  lf <- binomial(link = tempcall$link)
 				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
 				else {
 					mod1 <- glm(
@@ -136,11 +132,7 @@ ipwtm <- function(
 			}
 		#weights binomial, type "all"
 			if (tempcall$family == "binomial" & tempcall$type == "all") {
-				if(tempcall$link == "logit") lf <- binomial(link = logit)
-				if(tempcall$link == "probit") lf  <- binomial(link = probit)
-				if(tempcall$link == "cauchit") lf  <- binomial(link = cauchit)
-				if(tempcall$link == "log") lf  <- binomial(link = log)
-				if(tempcall$link == "cloglog") lf  <- binomial(link = cloglog)
+			  lf <- binomial(link = tempcall$link)
 				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
 				else {
 					mod1 <- glm(

From 6e29859478f97bc498472c6248ffce289fa8c821 Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Fri, 1 Nov 2024 12:10:06 +0700
Subject: [PATCH 03/14] Update exposure argument to ipwtm

Allow user to input character factor variable for exposure.
So, now the function
Add flow control for wrong format input.
---
 R/ipwtm.R | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)

diff --git a/R/ipwtm.R b/R/ipwtm.R
index d595a88..9b922ad 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -62,6 +62,57 @@ ipwtm <- function(
   
   # Sort dataframe by id and timevar
   data <- data[order.orig, ]
+  
+  #check if exposure is a binary variable of (0 and 1) or factor char variable 
+  check_exposure <- function(exposure) {
+    # Check if exposure has exactly 2 unique values
+    unique_vals <- unique(exposure)
+    num_unique_vals <- length(unique_vals)
+    
+    if (num_unique_vals > 2) {
+      stop(paste(
+        "Error: The 'exposure' variable must be binary, with only two unique values.\n",
+        "Currently, it has", num_unique_vals, "unique values:",
+        paste(unique_vals, collapse = ", "), "."
+      ))
+    } else if (num_unique_vals == 1) {
+      warning("Warning: The 'exposure' variable contains only one unique value.\nEnsure this is intentional.")
+    }
+    
+    # Check for numeric type with values 0 and 1 only
+    if (is.numeric(exposure)) {
+      non_binary_vals <- unique_vals[!unique_vals %in% c(0, 1)]
+      if (length(non_binary_vals) > 0) {
+        stop(paste(
+          "Error: Numeric 'exposure' values must be either 0 or 1.\n",
+          "Found invalid values:",
+          paste(non_binary_vals, collapse = ", "), "."
+        ))
+      }
+    }
+    
+    # Check for character type and factor requirement
+    else if (is.character(exposure)) {
+      if (!is.factor(exposure)) {
+        stop("Error: Character 'exposure' variable must be converted to a factor before use.")
+      } else {
+        warning("Notice: The factor levels of 'exposure' will be recoded.\nThe lower level will become 0, and the higher level will become 1.")
+        # Convert factor to 0 and 1 based on levels
+        exposure <- as.numeric(exposure) - 1
+      }
+    }
+    
+    # Error if exposure is not numeric or character
+    else {
+      stop("Error: The 'exposure' variable must be either a numeric variable with values 0 and 1, or a character variable that has been converted to a factor with exactly two levels.\nIf your data type is incorrect, please update it to meet these requirements for the function to work correctly.")
+    }
+    
+    # Return the modified exposure variable
+    return(exposure)
+  }
+  
+  data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]])
+  
   # Create a new dataframe to prevent variable name conflicts
   tempdat <- data.frame(
     id = data[[deparse(tempcall$id)]],

From 1e623e175edab71b8eb931e264143274d1674b1d Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Fri, 1 Nov 2024 13:28:26 +0700
Subject: [PATCH 04/14] debug check_exposure

debug check_exposure
---
 R/ipwtm.R | 34 +++++++++++++++++++++-------------
 1 file changed, 21 insertions(+), 13 deletions(-)

diff --git a/R/ipwtm.R b/R/ipwtm.R
index 9b922ad..f69c59a 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -91,26 +91,34 @@ ipwtm <- function(
       }
     }
     
-    # Check for character type and factor requirement
-    else if (is.character(exposure)) {
-      if (!is.factor(exposure)) {
-        stop("Error: Character 'exposure' variable must be converted to a factor before use.")
-      } else {
-        warning("Notice: The factor levels of 'exposure' will be recoded.\nThe lower level will become 0, and the higher level will become 1.")
-        # Convert factor to 0 and 1 based on levels
-        exposure <- as.numeric(exposure) - 1
+    # If non-numeric variable, check factor requirement
+    else if (is.factor(exposure)) {
+      # Extract factor levels
+      factor_levels <- levels(exposure)
+      
+      if (length(factor_levels) != 2) {
+        stop("Error: The factor 'exposure' must have exactly two levels.")
       }
-    }
-    
-    # Error if exposure is not numeric or character
-    else {
-      stop("Error: The 'exposure' variable must be either a numeric variable with values 0 and 1, or a character variable that has been converted to a factor with exactly two levels.\nIf your data type is incorrect, please update it to meet these requirements for the function to work correctly.")
+      
+      lower_level <- factor_levels[1]
+      higher_level <- factor_levels[2]
+      
+      warning(paste0(
+        "Notice: The factor levels of 'exposure' will be recoded.\n",
+        "The lower level (", lower_level, ") will be assigned as unexposed, and the higher level (", higher_level, ") will be assigned as exposed."
+      ))
+      
+      # Convert factor to 0 and 1 based on levels
+      exposure <- as.numeric(exposure) - 1
+    } else {
+      stop("Error: non-numeric 'exposure' variable must be converted to a factor before use.")
     }
     
     # Return the modified exposure variable
     return(exposure)
   }
   
+  
   data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]])
   
   # Create a new dataframe to prevent variable name conflicts

From 2dd5926135c4f75bf52a13fb9c3c1f87c90ee8ed Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Fri, 1 Nov 2024 14:09:38 +0700
Subject: [PATCH 05/14] move check function to separate file

move check function to separate file
---
 R/check_input.R | 54 ++++++++++++++++++++++++++++++++++++++++++++++
 R/ipwtm.R       | 57 +------------------------------------------------
 2 files changed, 55 insertions(+), 56 deletions(-)
 create mode 100644 R/check_input.R

diff --git a/R/check_input.R b/R/check_input.R
new file mode 100644
index 0000000..0a82b38
--- /dev/null
+++ b/R/check_input.R
@@ -0,0 +1,54 @@
+#check if exposure is a binary variable of (0 and 1) or factor char variable 
+check_exposure <- function(exposure) {
+  # Check if exposure has exactly 2 unique values
+  unique_vals <- unique(exposure)
+  num_unique_vals <- length(unique_vals)
+  
+  if (num_unique_vals > 2) {
+    stop(paste(
+      "Error: The 'exposure' variable must be binary, with only two unique values.\n",
+      "Currently, it has", num_unique_vals, "unique values:",
+      paste(unique_vals, collapse = ", "), "."
+    ))
+  } else if (num_unique_vals == 1) {
+    warning("Warning: The 'exposure' variable contains only one unique value.\nEnsure this is intentional.")
+  }
+  
+  # Check for numeric type with values 0 and 1 only
+  if (is.numeric(exposure)) {
+    non_binary_vals <- unique_vals[!unique_vals %in% c(0, 1)]
+    if (length(non_binary_vals) > 0) {
+      stop(paste(
+        "Error: Numeric 'exposure' values must be either 0 or 1.\n",
+        "Found invalid values:",
+        paste(non_binary_vals, collapse = ", "), "."
+      ))
+    }
+  }
+  
+  # If non-numeric variable, check factor requirement
+  else if (is.factor(exposure)) {
+    # Extract factor levels
+    factor_levels <- levels(exposure)
+    
+    if (length(factor_levels) != 2) {
+      stop("Error: The factor 'exposure' must have exactly two levels.")
+    }
+    
+    lower_level <- factor_levels[1]
+    higher_level <- factor_levels[2]
+    
+    warning(paste0(
+      "Notice: The factor levels of 'exposure' will be recoded.\n",
+      "The lower level (", lower_level, ") will be assigned as unexposed, and the higher level (", higher_level, ") will be assigned as exposed."
+    ))
+    
+    # Convert factor to 0 and 1 based on levels
+    exposure <- as.numeric(exposure) - 1
+  } else {
+    stop("Error: non-numeric 'exposure' variable must be converted to a factor before use.")
+  }
+  
+  # Return the modified exposure variable
+  return(exposure)
+}
\ No newline at end of file
diff --git a/R/ipwtm.R b/R/ipwtm.R
index f69c59a..7dc6dfd 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -63,62 +63,7 @@ ipwtm <- function(
   # Sort dataframe by id and timevar
   data <- data[order.orig, ]
   
-  #check if exposure is a binary variable of (0 and 1) or factor char variable 
-  check_exposure <- function(exposure) {
-    # Check if exposure has exactly 2 unique values
-    unique_vals <- unique(exposure)
-    num_unique_vals <- length(unique_vals)
-    
-    if (num_unique_vals > 2) {
-      stop(paste(
-        "Error: The 'exposure' variable must be binary, with only two unique values.\n",
-        "Currently, it has", num_unique_vals, "unique values:",
-        paste(unique_vals, collapse = ", "), "."
-      ))
-    } else if (num_unique_vals == 1) {
-      warning("Warning: The 'exposure' variable contains only one unique value.\nEnsure this is intentional.")
-    }
-    
-    # Check for numeric type with values 0 and 1 only
-    if (is.numeric(exposure)) {
-      non_binary_vals <- unique_vals[!unique_vals %in% c(0, 1)]
-      if (length(non_binary_vals) > 0) {
-        stop(paste(
-          "Error: Numeric 'exposure' values must be either 0 or 1.\n",
-          "Found invalid values:",
-          paste(non_binary_vals, collapse = ", "), "."
-        ))
-      }
-    }
-    
-    # If non-numeric variable, check factor requirement
-    else if (is.factor(exposure)) {
-      # Extract factor levels
-      factor_levels <- levels(exposure)
-      
-      if (length(factor_levels) != 2) {
-        stop("Error: The factor 'exposure' must have exactly two levels.")
-      }
-      
-      lower_level <- factor_levels[1]
-      higher_level <- factor_levels[2]
-      
-      warning(paste0(
-        "Notice: The factor levels of 'exposure' will be recoded.\n",
-        "The lower level (", lower_level, ") will be assigned as unexposed, and the higher level (", higher_level, ") will be assigned as exposed."
-      ))
-      
-      # Convert factor to 0 and 1 based on levels
-      exposure <- as.numeric(exposure) - 1
-    } else {
-      stop("Error: non-numeric 'exposure' variable must be converted to a factor before use.")
-    }
-    
-    # Return the modified exposure variable
-    return(exposure)
-  }
-  
-  
+  # Check if exposure is in correct format
   data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]])
   
   # Create a new dataframe to prevent variable name conflicts

From 1fb1186f73542a9efc80ab0ad79a23f7533c2d13 Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Fri, 1 Nov 2024 14:27:27 +0700
Subject: [PATCH 06/14] update ipwpoint

update check_exposure and concise check_requires
---
 R/ipwpoint.R | 56 ++++++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 46 insertions(+), 10 deletions(-)

diff --git a/R/ipwpoint.R b/R/ipwpoint.R
index cbf1f5d..e22d4ea 100644
--- a/R/ipwpoint.R
+++ b/R/ipwpoint.R
@@ -11,16 +11,52 @@ ipwpoint <- function(
 		#save input
 			tempcall <- match.call()
 		#some basic input checks
-			if (!("exposure" %in% names(tempcall))) stop("No exposure variable specified")
-			if (!("family" %in% names(tempcall)) | ("family" %in% names(tempcall) & !(tempcall$family %in% c("binomial", "multinomial", "ordinal", "gaussian")))) stop("No valid family specified (\"binomial\", \"multinomial\", \"ordinal\", \"gaussian\")")
-			if (tempcall$family == "binomial") {if(!(tempcall$link %in% c("logit", "probit", "cauchit", "log", "cloglog"))) stop("No valid link function specified for family = binomial (\"logit\", \"probit\", \"cauchit\", \"log\", \"cloglog\")")}
-			if (tempcall$family == "ordinal" ) {if(!(tempcall$link %in% c("logit", "probit", "cauchit", "cloglog"))) stop("No valid link function specified for family = binomial (\"logit\", \"probit\", \"cauchit\", \"cloglog\")")}
-			if (!("denominator" %in% names(tempcall))) stop("No denominator model specified")
-			if (!is.null(tempcall$numerator) & !is(eval(tempcall$numerator), "formula")) stop("Invalid numerator formula specified")
-			if (!is.null(tempcall$denominator) & !is(eval(tempcall$denominator), "formula")) stop("Invalid denominator formula specified")
-			if (tempcall$family %in% c("gaussian") & !("numerator" %in% names(tempcall))) stop("Numerator necessary for family = \"gaussian\"")
-			if (!("data" %in% names(tempcall))) stop("No data specified")
-			if (!is.null(tempcall$trunc)) {if(tempcall$trunc < 0 | tempcall$trunc > 0.5) stop("Invalid truncation percentage specified (0-0.5)")}
+			# Helper function to check if a variable exists in tempcall
+			check_var <- function(var, msg) {
+			  if (!(var %in% names(tempcall))) stop(msg)
+			}
+			
+			# Helper function to check if a value is in a set
+			check_in_set <- function(value, valid_set, msg) {
+			  if (!(value %in% valid_set)) stop(msg)
+			}
+			
+			# Check required variables
+			check_var("exposure", "No exposure variable specified")
+			check_var("denominator", "No denominator model specified")
+			check_var("data", "No data specified")
+			
+			# Check family validity
+			check_var("family", "No valid family specified (\"binomial\", \"multinomial\", \"ordinal\", \"gaussian\")")
+			check_in_set(tempcall$family, c("binomial", "multinomial", "ordinal", "gaussian"), "Invalid family specified")
+			
+			# Check link function for specific families
+			valid_links <- c("logit", "probit", "cauchit", "log", "cloglog")
+			if (tempcall$family %in% c("binomial", "ordinal")) {
+			  check_var("link", paste("No valid link function specified for family =", tempcall$family, "(", paste(valid_links, collapse = ", "), ")"))
+			  check_in_set(tempcall$link, valid_links, paste("No valid link function specified for family =", tempcall$family, "(", paste(valid_links, collapse = ", "), ")"))
+			}
+			
+			# Validate numerator and denominator formulas
+			if (!is.null(tempcall$numerator)) {
+			  if (!is(eval(tempcall$numerator), "formula")) stop("Invalid numerator formula specified")
+			}
+			if (!is.null(tempcall$denominator)) {
+			  if (!is(eval(tempcall$denominator), "formula")) stop("Invalid denominator formula specified")
+			}
+			
+			# Check for numerator in Gaussian family
+			if (tempcall$family == "gaussian") {
+			  check_var("numerator", "Numerator necessary for family = \"gaussian\"")
+			}
+			
+			# Validate truncation percentage
+			if (!is.null(tempcall$trunc)) {
+			  if (tempcall$trunc < 0 | tempcall$trunc > 0.5) stop("Invalid truncation percentage specified (0-0.5)")
+			}
+			
+		# Check if exposure is in correct format
+			data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]])
 		#make new dataframe for newly computed variables, to prevent variable name conflicts
 			tempdat <- data.frame(
 				exposure = data[,as.character(tempcall$exposure)]

From ea183ee9011340601ca40ab5b3b6e8c7cb96a8d8 Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Fri, 1 Nov 2024 14:49:35 +0700
Subject: [PATCH 07/14] make code more coherent

---
 R/ipwpoint.R | 18 +++++++++---------
 R/ipwtm.R    | 14 +++++++++++---
 2 files changed, 20 insertions(+), 12 deletions(-)

diff --git a/R/ipwpoint.R b/R/ipwpoint.R
index e22d4ea..00738e9 100644
--- a/R/ipwpoint.R
+++ b/R/ipwpoint.R
@@ -11,9 +11,9 @@ ipwpoint <- function(
 		#save input
 			tempcall <- match.call()
 		#some basic input checks
-			# Helper function to check if a variable exists in tempcall
-			check_var <- function(var, msg) {
-			  if (!(var %in% names(tempcall))) stop(msg)
+			# Helper function to check for required arguments
+			check_required <- function(arg, msg) {
+			  if (!(arg %in% names(tempcall))) stop(msg)
 			}
 			
 			# Helper function to check if a value is in a set
@@ -22,18 +22,18 @@ ipwpoint <- function(
 			}
 			
 			# Check required variables
-			check_var("exposure", "No exposure variable specified")
-			check_var("denominator", "No denominator model specified")
-			check_var("data", "No data specified")
+			check_required("exposure", "No exposure variable specified")
+			check_required("denominator", "No denominator model specified")
+			check_required("data", "No data specified")
 			
 			# Check family validity
-			check_var("family", "No valid family specified (\"binomial\", \"multinomial\", \"ordinal\", \"gaussian\")")
+			check_required("family", "No valid family specified (\"binomial\", \"multinomial\", \"ordinal\", \"gaussian\")")
 			check_in_set(tempcall$family, c("binomial", "multinomial", "ordinal", "gaussian"), "Invalid family specified")
 			
 			# Check link function for specific families
 			valid_links <- c("logit", "probit", "cauchit", "log", "cloglog")
 			if (tempcall$family %in% c("binomial", "ordinal")) {
-			  check_var("link", paste("No valid link function specified for family =", tempcall$family, "(", paste(valid_links, collapse = ", "), ")"))
+			  check_required("link", paste("No valid link function specified for family =", tempcall$family, "(", paste(valid_links, collapse = ", "), ")"))
 			  check_in_set(tempcall$link, valid_links, paste("No valid link function specified for family =", tempcall$family, "(", paste(valid_links, collapse = ", "), ")"))
 			}
 			
@@ -47,7 +47,7 @@ ipwpoint <- function(
 			
 			# Check for numerator in Gaussian family
 			if (tempcall$family == "gaussian") {
-			  check_var("numerator", "Numerator necessary for family = \"gaussian\"")
+			  check_required("numerator", "Numerator necessary for family = \"gaussian\"")
 			}
 			
 			# Validate truncation percentage
diff --git a/R/ipwtm.R b/R/ipwtm.R
index 7dc6dfd..4d8f2c0 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -10,6 +10,11 @@ ipwtm <- function(
     if (!(arg %in% names(tempcall))) stop(msg)
   }
   
+  # Helper function to check if a value is in a set
+  check_in_set <- function(value, valid_set, msg) {
+    if (!(value %in% valid_set)) stop(msg)
+  }
+
   # Validate required arguments
   check_required("exposure", "No exposure variable specified")
   check_required("denominator", "No denominator model specified")
@@ -25,13 +30,15 @@ ipwtm <- function(
     ordinal = c("logit", "probit", "cauchit", "cloglog")
   )
   
-  if (!family %in% valid_families) stop("Invalid family specified")
+  check_in_set(family, valid_families, "Invalid family specified")
   if (family %in% names(valid_links) && !(link %in% valid_links[[family]])) {
     stop(paste("No valid link function specified for family =", family))
   }
   
   # Special handling for survival
-  if (family == "survival") check_required("tstart", "No tstart specified, necessary for family = \"survival\"")
+  if (family == "survival") {
+    check_required("tstart", "No tstart specified, necessary for family = \"survival\"")
+  }
   
   # Formula validation
   if (!is.null(numerator) && !is(eval(numerator), "formula")) stop("Invalid numerator formula specified")
@@ -39,8 +46,9 @@ ipwtm <- function(
   
   # Type validation
   valid_types <- c("first", "all", "cens")
-  if (!(type %in% valid_types)) stop("Invalid type specified (\"first\", \"all\" or \"cens\")")
+  check_in_set(type, valid_types, "Invalid type specified (\"first\", \"all\" or \"cens\")")
   
+  # Unsupported combinations
   unsupported_combinations <- list(
     survival = "all",
     multinomial = c("all", "cens"),

From 9e3046ccffcf5c72d2e675dd934ae7d4b6b3ebe8 Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Tue, 5 Nov 2024 10:10:15 +0700
Subject: [PATCH 08/14] debug the ipwpoint

make binomial and ordinal family and link function check distinct
---
 R/ipwpoint.R | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/R/ipwpoint.R b/R/ipwpoint.R
index 00738e9..d0d738c 100644
--- a/R/ipwpoint.R
+++ b/R/ipwpoint.R
@@ -26,15 +26,16 @@ ipwpoint <- function(
 			check_required("denominator", "No denominator model specified")
 			check_required("data", "No data specified")
 			
-			# Check family validity
-			check_required("family", "No valid family specified (\"binomial\", \"multinomial\", \"ordinal\", \"gaussian\")")
-			check_in_set(tempcall$family, c("binomial", "multinomial", "ordinal", "gaussian"), "Invalid family specified")
+			# Family-specific checks
+			valid_families <- c("binomial", "multinomial", "ordinal", "gaussian")
+			valid_links <- list(
+			  binomial = c("logit", "probit", "cauchit", "log", "cloglog"),
+			  ordinal = c("logit", "probit", "cauchit", "cloglog")
+			)
 			
-			# Check link function for specific families
-			valid_links <- c("logit", "probit", "cauchit", "log", "cloglog")
-			if (tempcall$family %in% c("binomial", "ordinal")) {
-			  check_required("link", paste("No valid link function specified for family =", tempcall$family, "(", paste(valid_links, collapse = ", "), ")"))
-			  check_in_set(tempcall$link, valid_links, paste("No valid link function specified for family =", tempcall$family, "(", paste(valid_links, collapse = ", "), ")"))
+			check_in_set(family, valid_families, "Invalid family specified")
+			if (family %in% names(valid_links) && !(link %in% valid_links[[family]])) {
+			  stop(paste("No valid link function specified for family =", family))
 			}
 			
 			# Validate numerator and denominator formulas

From 5e1f8ba0b08cea16cd9da029dc728bbe9a18f4a1 Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Tue, 5 Nov 2024 14:53:40 +0700
Subject: [PATCH 09/14] change check_exposure

change check_exposure to make it match with each family
---
 R/check_input.R | 128 +++++++++++++++++++++++++++++++++---------------
 R/ipwpoint.R    |   2 +-
 R/ipwtm.R       |   2 +-
 3 files changed, 90 insertions(+), 42 deletions(-)

diff --git a/R/check_input.R b/R/check_input.R
index 0a82b38..c7ffbed 100644
--- a/R/check_input.R
+++ b/R/check_input.R
@@ -1,54 +1,102 @@
-#check if exposure is a binary variable of (0 and 1) or factor char variable 
-check_exposure <- function(exposure) {
-  # Check if exposure has exactly 2 unique values
+check_exposure <- function(exposure, family) {
   unique_vals <- unique(exposure)
   num_unique_vals <- length(unique_vals)
   
-  if (num_unique_vals > 2) {
-    stop(paste(
-      "Error: The 'exposure' variable must be binary, with only two unique values.\n",
-      "Currently, it has", num_unique_vals, "unique values:",
-      paste(unique_vals, collapse = ", "), "."
-    ))
-  } else if (num_unique_vals == 1) {
-    warning("Warning: The 'exposure' variable contains only one unique value.\nEnsure this is intentional.")
-  }
-  
-  # Check for numeric type with values 0 and 1 only
-  if (is.numeric(exposure)) {
-    non_binary_vals <- unique_vals[!unique_vals %in% c(0, 1)]
-    if (length(non_binary_vals) > 0) {
+  if (family == "binomial") {
+    # Binomial check: binary values (0, 1) or two levels if factor
+    if (num_unique_vals > 2) {
       stop(paste(
-        "Error: Numeric 'exposure' values must be either 0 or 1.\n",
-        "Found invalid values:",
-        paste(non_binary_vals, collapse = ", "), "."
+        "Error: For family 'binomial', the 'exposure' variable must be binary with only two unique values.\n",
+        "Currently, it has", num_unique_vals, "unique values:", 
+        paste(unique_vals, collapse = ", "), "."
       ))
+    } else if (num_unique_vals == 1) {
+      warning("Warning: The 'exposure' variable contains only one unique value.\nEnsure this is intentional.")
+    }
+    if (is.numeric(exposure)) {
+      non_binary_vals <- unique_vals[!unique_vals %in% c(0, 1)]
+      if (length(non_binary_vals) > 0) {
+        stop(paste(
+          "Error: Numeric 'exposure' values must be either 0 or 1.\n",
+          "Found invalid values:", paste(non_binary_vals, collapse = ", "), "."
+        ))
+      }
+    } else if (is.factor(exposure)) {
+      if (length(levels(exposure)) != 2) {
+        stop("Error: The factor 'exposure' must have exactly two levels.")
+      }
+      exposure <- as.numeric(exposure) - 1
+    } else {
+      stop("Error: non-numeric 'exposure' variable must be converted to a factor before use.")
     }
   }
   
-  # If non-numeric variable, check factor requirement
-  else if (is.factor(exposure)) {
-    # Extract factor levels
-    factor_levels <- levels(exposure)
+  else if (family == "survival") {
+    # Survival check: status indicator (0 for alive, 1 for dead, etc.)
+    valid_status_values <- c(0, 1, TRUE, FALSE, 2, 3)
     
-    if (length(factor_levels) != 2) {
-      stop("Error: The factor 'exposure' must have exactly two levels.")
+    if (!is.numeric(exposure) && !is.logical(exposure) && !is.factor(exposure)) {
+      stop("Error: For family 'survival', 'exposure' should be a numeric, logical, or factor variable representing the status indicator.")
     }
     
-    lower_level <- factor_levels[1]
-    higher_level <- factor_levels[2]
-    
-    warning(paste0(
-      "Notice: The factor levels of 'exposure' will be recoded.\n",
-      "The lower level (", lower_level, ") will be assigned as unexposed, and the higher level (", higher_level, ") will be assigned as exposed."
-    ))
-    
-    # Convert factor to 0 and 1 based on levels
-    exposure <- as.numeric(exposure) - 1
-  } else {
-    stop("Error: non-numeric 'exposure' variable must be converted to a factor before use.")
+    if (is.numeric(exposure)) {
+      # Check for valid numeric status indicators
+      non_valid_vals <- unique_vals[!unique_vals %in% valid_status_values]
+      if (length(non_valid_vals) > 0) {
+        stop(paste(
+          "Error: For family 'survival', 'exposure' values must be 0, 1, TRUE, FALSE, or 2/3 for interval censoring.\n",
+          "Found invalid values:", paste(non_valid_vals, collapse = ", "), "."
+        ))
+      }
+    } else if (is.logical(exposure)) {
+      # Convert logical TRUE/FALSE to numeric
+      exposure <- as.numeric(exposure)
+      valid_status_values <- as.numeric(valid_status_values)
+      non_valid_vals <- unique(exposure[!exposure %in% valid_status_values])
+      if (length(non_valid_vals) > 0) {
+        stop("Error: For family 'survival', logical exposure must represent 0 (alive) or 1 (dead).")
+      }
+    } else if (is.factor(exposure)) {
+      # Convert factor to numeric and check levels
+      factor_levels <- levels(exposure)
+      if (!all(factor_levels %in% as.character(valid_status_values))) {
+        stop("Error: The factor 'exposure' must have levels corresponding to valid status indicators (0, 1, TRUE, FALSE, 2, 3).")
+      }
+      exposure <- as.numeric(exposure)
+    }
+  }
+  
+  else if (family == "multinomial") {
+    # Multinomial check: categorical variable with more than 2 levels
+    if (!is.factor(exposure)) {
+      stop("Error: For family 'multinomial', 'exposure' should be a factor with more than two levels.")
+    }
+    if (num_unique_vals <= 2) {
+      stop("Error: For family 'multinomial', 'exposure' must have more than two unique values.")
+    }
+  }
+  
+  else if (family == "ordinal") {
+    # Ordinal check: ordered factor variable with more than 2 levels
+    if (!is.ordered(exposure)) {
+      stop("Error: For family 'ordinal', 'exposure' should be an ordered factor.")
+    }
+    if (num_unique_vals <= 2) {
+      stop("Error: For family 'ordinal', 'exposure' must have more than two unique values.")
+    }
+  }
+  
+  else if (family == "gaussian") {
+    # Gaussian check: continuous numeric variable
+    if (!is.numeric(exposure)) {
+      stop("Error: For family 'gaussian', 'exposure' should be a continuous numeric variable.")
+    }
+  }
+  
+  else {
+    stop("Error: Unknown family. Choose one of 'binomial', 'survival', 'multinomial', 'ordinal', or 'gaussian'.")
   }
   
-  # Return the modified exposure variable
+  # Return the modified or validated exposure variable
   return(exposure)
-}
\ No newline at end of file
+}
diff --git a/R/ipwpoint.R b/R/ipwpoint.R
index d0d738c..ce4839f 100644
--- a/R/ipwpoint.R
+++ b/R/ipwpoint.R
@@ -57,7 +57,7 @@ ipwpoint <- function(
 			}
 			
 		# Check if exposure is in correct format
-			data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]])
+			data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]], family)
 		#make new dataframe for newly computed variables, to prevent variable name conflicts
 			tempdat <- data.frame(
 				exposure = data[,as.character(tempcall$exposure)]
diff --git a/R/ipwtm.R b/R/ipwtm.R
index 4d8f2c0..39425ac 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -72,7 +72,7 @@ ipwtm <- function(
   data <- data[order.orig, ]
   
   # Check if exposure is in correct format
-  data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]])
+  data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]], family)
   
   # Create a new dataframe to prevent variable name conflicts
   tempdat <- data.frame(

From 5e9593507579fd441a1f387e1f1b475f94cdc09a Mon Sep 17 00:00:00 2001
From: TranHung93 <hungtt@oucru.org>
Date: Thu, 7 Nov 2024 12:13:54 +0700
Subject: [PATCH 10/14] Fix check_exposure

Adapt exposure requirement match to current response format in each model
---
 R/check_input.R | 37 ++++++++++++++++---------------------
 R/ipwpoint.R    |  4 ++--
 R/ipwtm.R       |  4 ++--
 3 files changed, 20 insertions(+), 25 deletions(-)

diff --git a/R/check_input.R b/R/check_input.R
index c7ffbed..cde56b6 100644
--- a/R/check_input.R
+++ b/R/check_input.R
@@ -11,7 +11,7 @@ check_exposure <- function(exposure, family) {
         paste(unique_vals, collapse = ", "), "."
       ))
     } else if (num_unique_vals == 1) {
-      warning("Warning: The 'exposure' variable contains only one unique value.\nEnsure this is intentional.")
+      warning("Warning: The 'exposure' variable contains only one unique value.\n Hence there is no contrast. Ensure this is intentional.")
     }
     if (is.numeric(exposure)) {
       non_binary_vals <- unique_vals[!unique_vals %in% c(0, 1)]
@@ -33,7 +33,15 @@ check_exposure <- function(exposure, family) {
   
   else if (family == "survival") {
     # Survival check: status indicator (0 for alive, 1 for dead, etc.)
-    valid_status_values <- c(0, 1, TRUE, FALSE, 2, 3)
+    if (num_unique_vals > 2) {
+      stop(paste(
+        "Error: For family 'survival', the 'exposure' variable must be binary with only two unique values.\n",
+        "Currently, it has", num_unique_vals, "unique values:", 
+        paste(unique_vals, collapse = ", "), "."
+      ))
+    } else if (num_unique_vals == 1) {
+      warning("Warning: The 'exposure' variable contains only one unique value.\n Hence there is no contrast. Ensure this is intentional.")
+    }
     
     if (!is.numeric(exposure) && !is.logical(exposure) && !is.factor(exposure)) {
       stop("Error: For family 'survival', 'exposure' should be a numeric, logical, or factor variable representing the status indicator.")
@@ -41,36 +49,23 @@ check_exposure <- function(exposure, family) {
     
     if (is.numeric(exposure)) {
       # Check for valid numeric status indicators
-      non_valid_vals <- unique_vals[!unique_vals %in% valid_status_values]
-      if (length(non_valid_vals) > 0) {
+      if (!all(unique_vals %in% c(0,1)) & !all(unique_vals  %in% c(1,2))) {
         stop(paste(
-          "Error: For family 'survival', 'exposure' values must be 0, 1, TRUE, FALSE, or 2/3 for interval censoring.\n",
-          "Found invalid values:", paste(non_valid_vals, collapse = ", "), "."
+          "Error: For family 'survival', numeric 'exposure' values must be (0 and 1) or (1 and 2).\n",
+          "Found invalid values:", paste(unique_vals, collapse = ", "), "."
         ))
       }
     } else if (is.logical(exposure)) {
       # Convert logical TRUE/FALSE to numeric
       exposure <- as.numeric(exposure)
-      valid_status_values <- as.numeric(valid_status_values)
-      non_valid_vals <- unique(exposure[!exposure %in% valid_status_values])
-      if (length(non_valid_vals) > 0) {
-        stop("Error: For family 'survival', logical exposure must represent 0 (alive) or 1 (dead).")
-      }
     } else if (is.factor(exposure)) {
       # Convert factor to numeric and check levels
-      factor_levels <- levels(exposure)
-      if (!all(factor_levels %in% as.character(valid_status_values))) {
-        stop("Error: The factor 'exposure' must have levels corresponding to valid status indicators (0, 1, TRUE, FALSE, 2, 3).")
-      }
       exposure <- as.numeric(exposure)
     }
   }
   
   else if (family == "multinomial") {
-    # Multinomial check: categorical variable with more than 2 levels
-    if (!is.factor(exposure)) {
-      stop("Error: For family 'multinomial', 'exposure' should be a factor with more than two levels.")
-    }
+    # Multinomial check: categorical variable with more than 2 unique values
     if (num_unique_vals <= 2) {
       stop("Error: For family 'multinomial', 'exposure' must have more than two unique values.")
     }
@@ -78,8 +73,8 @@ check_exposure <- function(exposure, family) {
   
   else if (family == "ordinal") {
     # Ordinal check: ordered factor variable with more than 2 levels
-    if (!is.ordered(exposure)) {
-      stop("Error: For family 'ordinal', 'exposure' should be an ordered factor.")
+    if (!is.factor(exposure)) {
+      stop("Error: For family 'ordinal', 'exposure' should be a factor.")
     }
     if (num_unique_vals <= 2) {
       stop("Error: For family 'ordinal', 'exposure' must have more than two unique values.")
diff --git a/R/ipwpoint.R b/R/ipwpoint.R
index ce4839f..9c11a09 100644
--- a/R/ipwpoint.R
+++ b/R/ipwpoint.R
@@ -12,8 +12,8 @@ ipwpoint <- function(
 			tempcall <- match.call()
 		#some basic input checks
 			# Helper function to check for required arguments
-			check_required <- function(arg, msg) {
-			  if (!(arg %in% names(tempcall))) stop(msg)
+			check_required <- function(arg, msg, arg_list = tempcall) {
+			  if (!(arg %in% names(arg_list))) stop(msg)
 			}
 			
 			# Helper function to check if a value is in a set
diff --git a/R/ipwtm.R b/R/ipwtm.R
index 39425ac..3e2c84d 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -6,8 +6,8 @@ ipwtm <- function(
   tempcall <- match.call()
   
   # Helper function to check for required arguments
-  check_required <- function(arg, msg) {
-    if (!(arg %in% names(tempcall))) stop(msg)
+  check_required <- function(arg, msg, arg_list = tempcall) {
+    if (!(arg %in% names(arg_list))) stop(msg)
   }
   
   # Helper function to check if a value is in a set

From 4038bdd5efc732710387bf22a47c7f1dc7723fec Mon Sep 17 00:00:00 2001
From: Hung Tran Thai <hungtt@oucru.org>
Date: Tue, 29 Apr 2025 16:17:39 +0700
Subject: [PATCH 11/14] shorten the tempcall link

---
 R/ipwtm.R | 19 +++++--------------
 1 file changed, 5 insertions(+), 14 deletions(-)

diff --git a/R/ipwtm.R b/R/ipwtm.R
index f1e3dc6..b2ad4d0 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -63,11 +63,7 @@ ipwtm <- function(
 				{tempdat$selvar <- rep(1, nrow(tempdat))}
 		#weights binomial, type "first"
 			if (tempcall$family == "binomial" & tempcall$type %in% c("first", "cens")) {
-				if(tempcall$link == "logit") lf <- binomial(link = logit)
-				if(tempcall$link == "probit") lf  <- binomial(link = probit)
-				if(tempcall$link == "cauchit") lf  <- binomial(link = cauchit)
-				if(tempcall$link == "log") lf  <- binomial(link = log)
-				if(tempcall$link == "cloglog") lf  <- binomial(link = cloglog)
+			  lf <- binomial(link = tempcall$link)
 				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
 				else {
 					mod1 <- glm(
@@ -109,11 +105,7 @@ ipwtm <- function(
 			}
 		#weights binomial, type "all"
 			if (tempcall$family == "binomial" & tempcall$type == "all") {
-				if(tempcall$link == "logit") lf <- binomial(link = logit)
-				if(tempcall$link == "probit") lf  <- binomial(link = probit)
-				if(tempcall$link == "cauchit") lf  <- binomial(link = cauchit)
-				if(tempcall$link == "log") lf  <- binomial(link = log)
-				if(tempcall$link == "cloglog") lf  <- binomial(link = cloglog)
+			  lf <- binomial(link = tempcall$link)
 				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
 				else {
 					mod1 <- glm(
@@ -238,10 +230,7 @@ ipwtm <- function(
 			}
 		#weights ordinal
 			if (tempcall$family == "ordinal") {
-				if(tempcall$link == "logit") m <- "logistic"
-				if(tempcall$link == "probit") m  <- "probit"
-				if(tempcall$link == "cloglog") m  <- "cloglog"
-				if(tempcall$link == "cauchit") m  <- "cauchit"
+				if(tempcall$link == "logit") m <- "logistic" else m <- tempcall$link
 				if (is.null(tempcall$numerator)) tempdat$p.numerator <- 1
 				else {
 					mod1 <- polr(
@@ -307,12 +296,14 @@ ipwtm <- function(
 			}
 		#check for NA's in weights
 			if (sum(is.na(tempdat$ipw.weights)) > 0) stop ("NA's in weights!")
+			
 		#truncate weights, when trunc value is specified (0-0.5)
 			if (!(is.null(tempcall$trunc))){
 				tempdat$weights.trunc <- tempdat$ipw.weights
 				tempdat$weights.trunc[tempdat$ipw.weights <= quantile(tempdat$ipw.weights, 0+trunc)] <- quantile(tempdat$ipw.weights, 0+trunc)
 				tempdat$weights.trunc[tempdat$ipw.weights >  quantile(tempdat$ipw.weights, 1-trunc)] <- quantile(tempdat$ipw.weights, 1-trunc)
 			}
+			
 		#return results in the same order as the original input dataframe
 			if (is.null(tempcall$trunc)){
 				if (is.null(tempcall$numerator)) return(list(ipw.weights = tempdat$ipw.weights[order(order.orig)], call = tempcall, selvar = tempdat$selvar[order(order.orig)], den.mod = mod2))

From 915b2ed370ca91746d9a40efa653375e40b0f3d4 Mon Sep 17 00:00:00 2001
From: Hung Tran Thai <hungtt@oucru.org>
Date: Mon, 5 May 2025 11:54:10 +0700
Subject: [PATCH 12/14] family to tempcall$family

---
 R/ipwtm.R | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/R/ipwtm.R b/R/ipwtm.R
index 3e2c84d..7480453 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -30,13 +30,13 @@ ipwtm <- function(
     ordinal = c("logit", "probit", "cauchit", "cloglog")
   )
   
-  check_in_set(family, valid_families, "Invalid family specified")
-  if (family %in% names(valid_links) && !(link %in% valid_links[[family]])) {
+  check_in_set(tempcall$family, valid_families, "Invalid family specified")
+  if (tempcall$family %in% names(valid_links) && !(tempcall$link %in% valid_links[[tempcall$family]])) {
     stop(paste("No valid link function specified for family =", family))
   }
   
   # Special handling for survival
-  if (family == "survival") {
+  if (tempcall$family == "survival") {
     check_required("tstart", "No tstart specified, necessary for family = \"survival\"")
   }
   
@@ -56,11 +56,11 @@ ipwtm <- function(
     gaussian = c("first", "cens")
   )
   
-  if (family %in% names(unsupported_combinations) && type %in% unsupported_combinations[[family]]) {
-    stop(paste("Type", type, "not yet implemented for family =", family))
+  if (tempcall$family %in% names(unsupported_combinations) && type %in% unsupported_combinations[[tempcall$family]]) {
+    stop(paste("Type", type, "not yet implemented for family =", tempcall$family))
   }
   
-  if (family == "gaussian" && is.null(numerator)) stop("Numerator necessary for family = \"gaussian\"")
+  if (tempcall$family == "gaussian" && is.null(numerator)) stop("Numerator necessary for family = \"gaussian\"")
   
   # Truncation check
   if (!is.null(trunc) && (trunc < 0 || trunc > 0.5)) stop("Invalid truncation percentage specified (0-0.5)")
@@ -72,7 +72,7 @@ ipwtm <- function(
   data <- data[order.orig, ]
   
   # Check if exposure is in correct format
-  data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]], family)
+  data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]], tempcall$family)
   
   # Create a new dataframe to prevent variable name conflicts
   tempdat <- data.frame(
@@ -83,12 +83,12 @@ ipwtm <- function(
   
   # Define selection variable based on type and family
   if (type %in% c("first", "cens")) {
-    if (family %in% c("binomial", "survival")) {
+    if (tempcall$family %in% c("binomial", "survival")) {
       tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {
         idx <- match(1, x)
         if (!is.na(idx)) c(rep(1, idx), rep(0, length(x) - idx)) else rep(1, length(x))
       }))
-    } else if (family %in% c("multinomial", "ordinal")) {
+    } else if (tempcall$family %in% c("multinomial", "ordinal")) {
       z <- setdiff(unique(tempdat$exposure), sort(unique(tempdat$exposure))[1])
       min2 <- function(x) ifelse(all(is.na(x)), NA, min(x, na.rm = TRUE))
       tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {

From 5eccca8cab84577bf3756410f994d1e5b979c9f1 Mon Sep 17 00:00:00 2001
From: Hung Tran Thai <hungtt@oucru.org>
Date: Mon, 5 May 2025 12:10:37 +0700
Subject: [PATCH 13/14] edit wrong push

---
 R/ipwtm.R | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/R/ipwtm.R b/R/ipwtm.R
index 7480453..3e2c84d 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -30,13 +30,13 @@ ipwtm <- function(
     ordinal = c("logit", "probit", "cauchit", "cloglog")
   )
   
-  check_in_set(tempcall$family, valid_families, "Invalid family specified")
-  if (tempcall$family %in% names(valid_links) && !(tempcall$link %in% valid_links[[tempcall$family]])) {
+  check_in_set(family, valid_families, "Invalid family specified")
+  if (family %in% names(valid_links) && !(link %in% valid_links[[family]])) {
     stop(paste("No valid link function specified for family =", family))
   }
   
   # Special handling for survival
-  if (tempcall$family == "survival") {
+  if (family == "survival") {
     check_required("tstart", "No tstart specified, necessary for family = \"survival\"")
   }
   
@@ -56,11 +56,11 @@ ipwtm <- function(
     gaussian = c("first", "cens")
   )
   
-  if (tempcall$family %in% names(unsupported_combinations) && type %in% unsupported_combinations[[tempcall$family]]) {
-    stop(paste("Type", type, "not yet implemented for family =", tempcall$family))
+  if (family %in% names(unsupported_combinations) && type %in% unsupported_combinations[[family]]) {
+    stop(paste("Type", type, "not yet implemented for family =", family))
   }
   
-  if (tempcall$family == "gaussian" && is.null(numerator)) stop("Numerator necessary for family = \"gaussian\"")
+  if (family == "gaussian" && is.null(numerator)) stop("Numerator necessary for family = \"gaussian\"")
   
   # Truncation check
   if (!is.null(trunc) && (trunc < 0 || trunc > 0.5)) stop("Invalid truncation percentage specified (0-0.5)")
@@ -72,7 +72,7 @@ ipwtm <- function(
   data <- data[order.orig, ]
   
   # Check if exposure is in correct format
-  data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]], tempcall$family)
+  data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]], family)
   
   # Create a new dataframe to prevent variable name conflicts
   tempdat <- data.frame(
@@ -83,12 +83,12 @@ ipwtm <- function(
   
   # Define selection variable based on type and family
   if (type %in% c("first", "cens")) {
-    if (tempcall$family %in% c("binomial", "survival")) {
+    if (family %in% c("binomial", "survival")) {
       tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {
         idx <- match(1, x)
         if (!is.na(idx)) c(rep(1, idx), rep(0, length(x) - idx)) else rep(1, length(x))
       }))
-    } else if (tempcall$family %in% c("multinomial", "ordinal")) {
+    } else if (family %in% c("multinomial", "ordinal")) {
       z <- setdiff(unique(tempdat$exposure), sort(unique(tempdat$exposure))[1])
       min2 <- function(x) ifelse(all(is.na(x)), NA, min(x, na.rm = TRUE))
       tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {

From b24c9a43179768d4535a03b8a08104d065e45308 Mon Sep 17 00:00:00 2001
From: Hung Tran Thai <hungtt@oucru.org>
Date: Mon, 5 May 2025 12:21:03 +0700
Subject: [PATCH 14/14] reformat code

---
 R/ipwtm.R | 909 +++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 634 insertions(+), 275 deletions(-)

diff --git a/R/ipwtm.R b/R/ipwtm.R
index 6606fd9..b40f61d 100644
--- a/R/ipwtm.R
+++ b/R/ipwtm.R
@@ -1,20 +1,31 @@
-ipwtm <- function(
-    exposure, family, link, denominator, id, tstart, timevar, type, data, 
-    numerator = NULL, corstr = "ar1", trunc = NULL, ...
-) {
+ipwtm <- function(exposure,
+                  family,
+                  link,
+                  denominator,
+                  id,
+                  tstart,
+                  timevar,
+                  type,
+                  data,
+                  numerator = NULL,
+                  corstr = "ar1",
+                  trunc = NULL,
+                  ...) {
   # Save input
   tempcall <- match.call()
   
   # Helper function to check for required arguments
   check_required <- function(arg, msg, arg_list = tempcall) {
-    if (!(arg %in% names(arg_list))) stop(msg)
+    if (!(arg %in% names(arg_list)))
+      stop(msg)
   }
   
   # Helper function to check if a value is in a set
   check_in_set <- function(value, valid_set, msg) {
-    if (!(value %in% valid_set)) stop(msg)
+    if (!(value %in% valid_set))
+      stop(msg)
   }
-
+  
   # Validate required arguments
   check_required("exposure", "No exposure variable specified")
   check_required("denominator", "No denominator model specified")
@@ -31,22 +42,29 @@ ipwtm <- function(
   )
   
   check_in_set(family, valid_families, "Invalid family specified")
-  if (family %in% names(valid_links) && !(link %in% valid_links[[family]])) {
+  if (family %in% names(valid_links) &&
+      !(link %in% valid_links[[family]])) {
     stop(paste("No valid link function specified for family =", family))
   }
   
   # Special handling for survival
   if (family == "survival") {
-    check_required("tstart", "No tstart specified, necessary for family = \"survival\"")
+    check_required("tstart",
+                   "No tstart specified, necessary for family = \"survival\"")
   }
   
   # Formula validation
-  if (!is.null(numerator) && !is(eval(numerator), "formula")) stop("Invalid numerator formula specified")
-  if (!is(eval(denominator), "formula")) stop("Invalid denominator formula specified")
+  if (!is.null(numerator) &&
+      !is(eval(numerator), "formula"))
+    stop("Invalid numerator formula specified")
+  if (!is(eval(denominator), "formula"))
+    stop("Invalid denominator formula specified")
   
   # Type validation
   valid_types <- c("first", "all", "cens")
-  check_in_set(type, valid_types, "Invalid type specified (\"first\", \"all\" or \"cens\")")
+  check_in_set(type,
+               valid_types,
+               "Invalid type specified (\"first\", \"all\" or \"cens\")")
   
   # Unsupported combinations
   unsupported_combinations <- list(
@@ -56,14 +74,20 @@ ipwtm <- function(
     gaussian = c("first", "cens")
   )
   
-  if (family %in% names(unsupported_combinations) && type %in% unsupported_combinations[[family]]) {
+  if (family %in% names(unsupported_combinations) &&
+      type %in% unsupported_combinations[[family]]) {
     stop(paste("Type", type, "not yet implemented for family =", family))
   }
   
-  if (family == "gaussian" && is.null(numerator)) stop("Numerator necessary for family = \"gaussian\"")
+  if (family == "gaussian" &&
+      is.null(numerator))
+    stop("Numerator necessary for family = \"gaussian\"")
   
   # Truncation check
-  if (!is.null(trunc) && (trunc < 0 || trunc > 0.5)) stop("Invalid truncation percentage specified (0-0.5)")
+  if (!is.null(trunc) &&
+      (trunc < 0 ||
+       trunc > 0.5))
+    stop("Invalid truncation percentage specified (0-0.5)")
   
   # Record original order of dataframe
   order.orig <- order(data[[deparse(tempcall$id)]], data[[deparse(tempcall$timevar)]])
@@ -75,281 +99,616 @@ ipwtm <- function(
   data[[deparse(tempcall$exposure)]] <- check_exposure(data[[deparse(tempcall$exposure)]], family)
   
   # Create a new dataframe to prevent variable name conflicts
-  tempdat <- data.frame(
-    id = data[[deparse(tempcall$id)]],
-    timevar = data[[deparse(tempcall$timevar)]],
-    exposure = data[[deparse(tempcall$exposure)]]
-  )
+  tempdat <- data.frame(id = data[[deparse(tempcall$id)]],
+                        timevar = data[[deparse(tempcall$timevar)]],
+                        exposure = data[[deparse(tempcall$exposure)]])
   
   # Define selection variable based on type and family
   if (type %in% c("first", "cens")) {
     if (family %in% c("binomial", "survival")) {
       tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {
         idx <- match(1, x)
-        if (!is.na(idx)) c(rep(1, idx), rep(0, length(x) - idx)) else rep(1, length(x))
+        if (!is.na(idx))
+          c(rep(1, idx), rep(0, length(x) - idx))
+        else
+          rep(1, length(x))
       }))
     } else if (family %in% c("multinomial", "ordinal")) {
       z <- setdiff(unique(tempdat$exposure), sort(unique(tempdat$exposure))[1])
-      min2 <- function(x) ifelse(all(is.na(x)), NA, min(x, na.rm = TRUE))
+      min2 <- function(x)
+        ifelse(all(is.na(x)), NA, min(x, na.rm = TRUE))
       tempdat$selvar <- unlist(lapply(split(tempdat$exposure, tempdat$id), function(x) {
         idx <- min2(match(z, x))
-        if (!is.na(idx)) c(rep(1, idx), rep(0, length(x) - idx)) else rep(1, length(x))
+        if (!is.na(idx))
+          c(rep(1, idx), rep(0, length(x) - idx))
+        else
+          rep(1, length(x))
       }))
     }
   } else if (type == "all") {
     tempdat$selvar <- rep(1, nrow(tempdat))
   }
   
-		#weights binomial, type "first"
-			if (tempcall$family == "binomial" & tempcall$type %in% c("first", "cens")) {
-			  lf <- binomial(link = tempcall$link)
-				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
-				else {
-					mod1 <- glm(
-						formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = ""))),
-						family = lf,
-						data = data,
-						subset = tempdat$selvar == 1,
-						na.action = na.fail,
-						...)
-					tempdat$p.numerator <- vector("numeric", nrow(tempdat))
-					tempdat$p.numerator[tempdat$exposure == 0 & tempdat$selvar == 1] <- 1 - predict.glm(mod1, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 0]
-					if(type == "first"){tempdat$p.numerator[tempdat$exposure == 1 & tempdat$selvar == 1] <- predict.glm(mod1, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]}
-					if(type == "cens"){tempdat$p.numerator[tempdat$exposure == 1 & tempdat$selvar == 1] <- 1 - predict.glm(mod1, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]}
-					tempdat$p.numerator[tempdat$selvar == 0] <- 1
-					tempdat$w.numerator <- unlist(lapply(split(tempdat$p.numerator, tempdat$id), function(x)cumprod(x)))
-						mod1$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = "")))
-						mod1$call$family <- tempcall$link
-						mod1$call$data <- tempcall$data
-						mod1$call$subset <- paste("up to first instance of ", deparse(tempcall$exposure, width.cutoff = 500), " = 1 (selvar == 1)", sep = "")
-				}
-				mod2 <- glm(
-					formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = ""))),
-					family = lf,
-					data = data,
-					subset = tempdat$selvar == 1,
-					na.action = na.fail,
-					...)
-				tempdat$p.denominator <- vector("numeric", nrow(tempdat))
-				tempdat$p.denominator[tempdat$exposure == 0 & tempdat$selvar == 1] <- 1 - predict.glm(mod2, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 0]
-				if(type == "first"){tempdat$p.denominator[tempdat$exposure == 1 & tempdat$selvar == 1] <- predict.glm(mod2, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]}
-				if(type == "cens"){tempdat$p.denominator[tempdat$exposure == 1 & tempdat$selvar == 1] <- 1 - predict.glm(mod2, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]}
-				tempdat$p.denominator[tempdat$selvar == 0] <- 1
-				tempdat$w.denominator <- unlist(lapply(split(tempdat$p.denominator, tempdat$id), function(x)cumprod(x)))
-				mod2$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = "")))
-				mod2$call$family <- tempcall$link
-				mod2$call$data <- tempcall$data
-				mod2$call$subset <- paste("up to first instance of ", deparse(tempcall$exposure, width.cutoff = 500), " = 1 (selvar == 1)", sep = "")
-				tempdat$ipw.weights <- tempdat$w.numerator/tempdat$w.denominator
-			}
-		#weights binomial, type "all"
-			if (tempcall$family == "binomial" & tempcall$type == "all") {
-			  lf <- binomial(link = tempcall$link)
-				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
-				else {
-					mod1 <- glm(
-						formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = ""))),
-						family = lf,
-						data = data,
-						na.action = na.fail,
-						...)
-					tempdat$p.numerator <- vector("numeric", nrow(tempdat))
-					tempdat$p.numerator[tempdat$exposure == 0] <- 1 - predict.glm(mod1, type = "response")[tempdat$exposure == 0]
-					tempdat$p.numerator[tempdat$exposure == 1] <- predict.glm(mod1, type = "response")[tempdat$exposure == 1]
-					tempdat$w.numerator <- unlist(lapply(split(tempdat$p.numerator, tempdat$id), function(x)cumprod(x)))
-						mod1$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = "")))
-						mod1$call$family <- tempcall$link
-						mod1$call$data <- tempcall$data
-				}
-				mod2 <- glm(
-					formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = ""))),
-					family = lf,
-					data = data,
-					na.action = na.fail,
-					...)
-				tempdat$p.denominator <- vector("numeric", nrow(tempdat))
-				tempdat$p.denominator[tempdat$exposure == 0] <- 1 - predict.glm(mod2, type = "response")[tempdat$exposure == 0]
-				tempdat$p.denominator[tempdat$exposure == 1] <- predict.glm(mod2, type = "response")[tempdat$exposure == 1]
-				tempdat$w.denominator <- unlist(lapply(split(tempdat$p.denominator, tempdat$id), function(x)cumprod(x)))
-				mod2$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = "")))
-				mod2$call$family <- tempcall$link
-				mod2$call$data <- tempcall$data
-				tempdat$ipw.weights <- tempdat$w.numerator/tempdat$w.denominator
-			}
-		#weights Cox
-			if (tempcall$family == "survival") {
-				if (is.null(tempcall$numerator)) tempdat$w.numerator <- 1
-					else {
-					mod1 <- coxph(
-						formula = eval(parse(text = paste("Surv(", deparse(tempcall$tstart), ", ", deparse(tempcall$timevar, width.cutoff = 500), ", ", deparse(tempcall$exposure, width.cutoff = 500), ") ", deparse(tempcall$numerator, width.cutoff = 500), sep = ""))),
-						data = data,
-						subset = tempdat$selvar == 1,
-						na.action = na.fail,
-						method = "efron",
-						...)
-					bh <- basehaz(mod1, centered = TRUE)
-					temp <- data.frame(timevar = sort(unique(tempdat$timevar)))
-					temp <- merge(temp, data.frame(timevar = bh$time, bashaz.cum.numerator = bh$hazard), by = "timevar", all.x = TRUE);rm(bh)
-					if (is.na(temp$bashaz.cum.numerator[temp$timevar == min(unique(tempdat$timevar))])) temp$bashaz.cum.numerator[temp$timevar == min(unique(tempdat$timevar))] <- 0
-					temp$bashaz.cum.numerator <- approx(x = temp$timevar, y = temp$bashaz.cum.numerator, xout = temp$timevar, method = "constant", rule = 2)$y
-					temp$bashaz.numerator[1] <- temp$bashaz.cum.numerator[1]
-					temp$bashaz.numerator[2:nrow(temp)] <- diff(temp$bashaz.cum.numerator, 1)
-					temp$bashaz.cum.numerator <- NULL
-					tempdat <- merge(tempdat, temp, by = "timevar", all.x = TRUE);rm(temp)
-					tempdat <- tempdat[order(tempdat$id, tempdat$timevar),]
-					tempdat$risk.numerator[tempdat$selvar == 1] <-predict(mod1, type="risk", centered = TRUE)
-					tempdat$hazard.numerator[tempdat$selvar == 1] <- with(tempdat[tempdat$selvar == 1,], bashaz.numerator*risk.numerator)
-					tempdat$p.numerator[with(tempdat, selvar == 1 & exposure == 0)] <- with(tempdat[with(tempdat, selvar == 1 & exposure == 0),], exp(-1*bashaz.numerator*risk.numerator))
-					if(type == "first"){tempdat$p.numerator[with(tempdat, selvar == 1 & exposure == 1)] <- 1 - with(tempdat[with(tempdat, selvar == 1 & exposure == 1),], exp(-1*bashaz.numerator*risk.numerator))}
-					if(type == "cens"){tempdat$p.numerator[with(tempdat, selvar == 1 & exposure == 1)] <- with(tempdat[with(tempdat, selvar == 1 & exposure == 1),], exp(-1*bashaz.numerator*risk.numerator))}
-					tempdat$p.numerator[tempdat$selvar == 0] <- 1
-					tempdat$w.numerator <- unsplit(lapply(split(tempdat$p.numerator, tempdat$id), function(x)cumprod(x)), tempdat$id)
-					mod1$call$formula <- eval(parse(text = paste("Surv(", deparse(tempcall$tstart), ", ", deparse(tempcall$timevar, width.cutoff = 500), ", ", deparse(tempcall$exposure, width.cutoff = 500), ") ", deparse(tempcall$numerator, width.cutoff = 500), sep = "")))
-					mod1$call$data <- tempcall$data
-				}
-				mod2 <- coxph(
-					formula = eval(parse(text = paste("Surv(", deparse(tempcall$tstart), ", ", deparse(tempcall$timevar, width.cutoff = 500), ", ", deparse(tempcall$exposure, width.cutoff = 500), ") ", deparse(tempcall$denominator, width.cutoff = 500), sep = ""))),
-					data = data,
-					subset = tempdat$selvar == 1,
-					na.action = na.fail,
-					method = "efron",
-					...)
-				bh <- basehaz(mod2, centered = TRUE)
-				temp <- data.frame(timevar = sort(unique(tempdat$timevar)))
-				temp <- merge(temp, data.frame(timevar = bh$time, bashaz.cum.denominator = bh$hazard), by = "timevar", all.x = TRUE);rm(bh)
-				if (is.na(temp$bashaz.cum.denominator[temp$timevar == min(unique(tempdat$timevar))])) temp$bashaz.cum.denominator[temp$timevar == min(unique(tempdat$timevar))] <- 0
-				temp$bashaz.cum.denominator <- approx(x = temp$timevar, y = temp$bashaz.cum.denominator, xout = temp$timevar, method = "constant", rule = 2)$y
-				temp$bashaz.denominator[1] <- temp$bashaz.cum.denominator[1]
-				temp$bashaz.denominator[2:nrow(temp)] <- diff(temp$bashaz.cum.denominator, 1)
-				temp$bashaz.cum.denominator <- NULL
-				tempdat <- merge(tempdat, temp, by = "timevar", all.x = TRUE);rm(temp)
-				tempdat <- tempdat[order(tempdat$id, tempdat$timevar),]
-				tempdat$risk.denominator[tempdat$selvar == 1] <-predict(mod2, type="risk", centered = TRUE)
-				tempdat$hazard.denominator[tempdat$selvar == 1] <- with(tempdat[tempdat$selvar == 1,], bashaz.denominator*risk.denominator)
-				tempdat$p.denominator[with(tempdat, selvar == 1 & exposure == 0)] <- with(tempdat[with(tempdat, selvar == 1 & exposure == 0),], exp(-1*bashaz.denominator*risk.denominator))
-				if(type == "first"){tempdat$p.denominator[with(tempdat, selvar == 1 & exposure == 1)] <- 1 - with(tempdat[with(tempdat, selvar == 1 & exposure == 1),], exp(-1*bashaz.denominator*risk.denominator))}
-				if(type == "cens"){tempdat$p.denominator[with(tempdat, selvar == 1 & exposure == 1)] <- with(tempdat[with(tempdat, selvar == 1 & exposure == 1),], exp(-1*bashaz.denominator*risk.denominator))}
-				tempdat$p.denominator[tempdat$selvar == 0] <- 1
-				tempdat$w.denominator <- unsplit(lapply(split(tempdat$p.denominator, tempdat$id), function(x)cumprod(x)), tempdat$id)
-				mod2$call$formula <- eval(parse(text = paste("Surv(", deparse(tempcall$tstart), ", ", deparse(tempcall$timevar, width.cutoff = 500), ", ", deparse(tempcall$exposure, width.cutoff = 500), ") ", deparse(tempcall$denominator, width.cutoff = 500), sep = "")))
-				mod2$call$data <- tempcall$data
-				mod2$call$subset <- paste("up to first instance of ", deparse(tempcall$exposure, width.cutoff = 500), " = 1 (selvar == 1)", sep = "")
-				tempdat$ipw.weights <- tempdat$w.numerator/tempdat$w.denominator
-			}
-		#weights multinomial
-			if (tempcall$family == "multinomial") {
-				if (is.null(tempcall$numerator)) tempdat$p.numerator <- 1
-				else {
-					mod1 <- multinom(
-						formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = ""))),
-						data = data,
-						subset = tempdat$selvar == 1,
-						na.action = na.fail,
-						...)
-					pred1 <- as.data.frame(predict(mod1, type = "probs"))
-					tempdat$p.numerator[tempdat$selvar == 0] <- 1
-					for (i in 1:length(unique(tempdat$exposure)))tempdat$p.numerator[with(tempdat, tempdat$selvar == 1 & exposure == sort(unique(tempdat$exposure))[i])] <- pred1[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i],i]
-					mod1$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = "")))
-					mod1$call$data <- tempcall$data
-					mod1$call$subset <- paste("up to first instance of ", deparse(tempcall$exposure, width.cutoff = 500), " = 1 (selvar == 1)", sep = "")
-				}
-				mod2 <- multinom(
-					formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = ""))),
-					data = data,
-					subset = tempdat$selvar == 1,
-					na.action = na.fail,
-					...)
-				pred2 <- as.data.frame(predict(mod2, type = "probs"))
-				tempdat$p.denominator[tempdat$selvar == 0] <- 1
-				for (i in 1:length(unique(tempdat$exposure)))tempdat$p.denominator[with(tempdat, tempdat$selvar == 1 & exposure == sort(unique(tempdat$exposure))[i])] <- pred2[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i],i]
-				mod2$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = "")))
-				mod2$call$data <- tempcall$data
-				mod2$call$subset <- paste("up to first instance of ", deparse(tempcall$exposure, width.cutoff = 500), " = 1 (selvar == 1)", sep = "")
-				tempdat$ipw.weights <- unsplit(lapply(split(with(tempdat, p.numerator/p.denominator), tempdat$id), function(x)cumprod(x)), tempdat$id)
-			}
-		#weights ordinal
-			if (tempcall$family == "ordinal") {
-				if(tempcall$link == "logit") m <- "logistic" else m <- tempcall$link
-				if (is.null(tempcall$numerator)) tempdat$p.numerator <- 1
-				else {
-					mod1 <- polr(
-						formula = eval(parse(text = paste("as.factor(", deparse(tempcall$exposure, width.cutoff = 500), ")", deparse(tempcall$numerator, width.cutoff = 500), sep = ""))),
-						data = data,
-						method = m,
-						subset = tempdat$selvar == 1,
-						na.action = na.fail,
-						...)
-					pred1 <- as.data.frame(predict(mod1, type = "probs"))
-					tempdat$p.numerator[tempdat$selvar == 0] <- 1
-					for (i in 1:length(unique(tempdat$exposure)))tempdat$p.numerator[with(tempdat, tempdat$selvar == 1 & exposure == sort(unique(tempdat$exposure))[i])] <- pred1[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i],i]
-					mod1$call$formula <- eval(parse(text = paste("as.factor(", deparse(tempcall$exposure, width.cutoff = 500), ")", deparse(tempcall$numerator, width.cutoff = 500), sep = "")))
-					mod1$call$data <- tempcall$data
-					mod1$call$method <- m
-					mod1$call$subset <- paste("up to first instance of ", deparse(tempcall$exposure, width.cutoff = 500), " = 1 (selvar == 1)", sep = "")
-				}
-				mod2 <- polr(
-					formula = eval(parse(text = paste("as.factor(", deparse(tempcall$exposure, width.cutoff = 500), ")", deparse(tempcall$denominator, width.cutoff = 500), sep = ""))),
-					data = data,
-					method = m,
-					subset = tempdat$selvar == 1,
-					na.action = na.fail,
-					...)
-				pred2 <- as.data.frame(predict(mod2, type = "probs"))
-				tempdat$p.denominator[tempdat$selvar == 0] <- 1
-				for (i in 1:length(unique(tempdat$exposure)))tempdat$p.denominator[with(tempdat, tempdat$selvar == 1 & exposure == sort(unique(tempdat$exposure))[i])] <- pred2[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i],i]
-				mod2$call$formula <- eval(parse(text = paste("as.factor(", deparse(tempcall$exposure, width.cutoff = 500), ")", deparse(tempcall$denominator, width.cutoff = 500), sep = "")))
-				mod2$call$data <- tempcall$data
-				mod2$call$method <- m
-				mod2$call$subset <- paste("up to first instance of ", deparse(tempcall$exposure, width.cutoff = 500), " = 1 (selvar == 1)", sep = "")
-				tempdat$ipw.weights <- unsplit(lapply(split(with(tempdat, p.numerator/p.denominator), tempdat$id), function(x)cumprod(x)), tempdat$id)
-			}
-		#weights gaussian
-			if (tempcall$family == "gaussian") {
-				mod1 <- geeglm(
-					formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = ""))),
-					data = data,
-					id = tempdat$id,
-					corstr = tempcall$corstr,
-					waves = tempdat$timevar,
-					...)
-				mod1$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$numerator, width.cutoff = 500), sep = "")))
-				mod1$call$data <- tempcall$data
-				mod1$call$id <- tempcall$id
-				mod1$call$corstr <- tempcall$corstr
-				mod1$call$waves <- tempcall$waves
-				mod2 <- geeglm(
-					formula = eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = ""))),
-					data = data,
-					id = tempdat$id,
-					corstr = tempcall$corstr,
-					waves = tempdat$timevar,
-					...)
-				mod2$call$formula <- eval(parse(text = paste(deparse(tempcall$exposure, width.cutoff = 500), deparse(tempcall$denominator, width.cutoff = 500), sep = "")))
-				mod2$call$data <- tempcall$data
-				mod2$call$id <- tempcall$id
-				mod2$call$corstr <- tempcall$corstr
-				mod2$call$waves <- tempcall$waves
-				tempdat$kdens1 <- dnorm(tempdat$exposure, predict(mod1), as.numeric(sqrt(summary(mod1)$dispersion)[1]))
-				tempdat$kdens2 <- dnorm(tempdat$exposure, predict(mod2), as.numeric(sqrt(summary(mod2)$dispersion)[1]))
-				tempdat$ipw.weights <- unsplit(lapply(split(with(tempdat, kdens1/kdens2), tempdat$id), function(x)cumprod(x)), tempdat$id)
-			}
-		#check for NA's in weights
-			if (sum(is.na(tempdat$ipw.weights)) > 0) stop ("NA's in weights!")
-			
-		#truncate weights, when trunc value is specified (0-0.5)
-			if (!(is.null(tempcall$trunc))){
-				tempdat$weights.trunc <- tempdat$ipw.weights
-				tempdat$weights.trunc[tempdat$ipw.weights <= quantile(tempdat$ipw.weights, 0+trunc)] <- quantile(tempdat$ipw.weights, 0+trunc)
-				tempdat$weights.trunc[tempdat$ipw.weights >  quantile(tempdat$ipw.weights, 1-trunc)] <- quantile(tempdat$ipw.weights, 1-trunc)
-			}
-			
-		#return results in the same order as the original input dataframe
-			if (is.null(tempcall$trunc)){
-				if (is.null(tempcall$numerator)) return(list(ipw.weights = tempdat$ipw.weights[order(order.orig)], call = tempcall, selvar = tempdat$selvar[order(order.orig)], den.mod = mod2))
-				else return(list(ipw.weights = tempdat$ipw.weights[order(order.orig)], call = tempcall, selvar = tempdat$selvar[order(order.orig)], num.mod = mod1, den.mod = mod2))
-			}
-			else{
-				if (is.null(tempcall$numerator)) return(list(ipw.weights = tempdat$ipw.weights[order(order.orig)], weights.trunc = tempdat$weights.trunc[order(order.orig)], call = tempcall, selvar = tempdat$selvar[order(order.orig)], den.mod = mod2))
-				else return(list(ipw.weights = tempdat$ipw.weights[order(order.orig)], weights.trunc = tempdat$weights.trunc[order(order.orig)], call = tempcall, selvar = tempdat$selvar[order(order.orig)], num.mod = mod1, den.mod = mod2))
-			}
+  #weights binomial, type "first"
+  if (tempcall$family == "binomial" &
+      tempcall$type %in% c("first", "cens")) {
+    lf <- binomial(link = tempcall$link)
+    if (is.null(tempcall$numerator))
+      tempdat$w.numerator <- 1
+    else {
+      mod1 <- glm(
+        formula = eval(parse(text = paste(
+          deparse(tempcall$exposure, width.cutoff = 500),
+          deparse(tempcall$numerator, width.cutoff = 500),
+          sep = ""
+        ))),
+        family = lf,
+        data = data,
+        subset = tempdat$selvar == 1,
+        na.action = na.fail,
+        ...
+      )
+      tempdat$p.numerator <- vector("numeric", nrow(tempdat))
+      tempdat$p.numerator[tempdat$exposure == 0 &
+                            tempdat$selvar == 1] <- 1 - predict.glm(mod1, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 0]
+      if (type == "first") {
+        tempdat$p.numerator[tempdat$exposure == 1 &
+                              tempdat$selvar == 1] <- predict.glm(mod1, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]
+      }
+      if (type == "cens") {
+        tempdat$p.numerator[tempdat$exposure == 1 &
+                              tempdat$selvar == 1] <- 1 - predict.glm(mod1, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]
+      }
+      tempdat$p.numerator[tempdat$selvar == 0] <- 1
+      tempdat$w.numerator <- unlist(lapply(split(tempdat$p.numerator, tempdat$id), function(x)
+        cumprod(x)))
+      mod1$call$formula <- eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$numerator, width.cutoff = 500),
+        sep = ""
+      )))
+      mod1$call$family <- tempcall$link
+      mod1$call$data <- tempcall$data
+      mod1$call$subset <- paste(
+        "up to first instance of ",
+        deparse(tempcall$exposure, width.cutoff = 500),
+        " = 1 (selvar == 1)",
+        sep = ""
+      )
+    }
+    mod2 <- glm(
+      formula = eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$denominator, width.cutoff = 500),
+        sep = ""
+      ))),
+      family = lf,
+      data = data,
+      subset = tempdat$selvar == 1,
+      na.action = na.fail,
+      ...
+    )
+    tempdat$p.denominator <- vector("numeric", nrow(tempdat))
+    tempdat$p.denominator[tempdat$exposure == 0 &
+                            tempdat$selvar == 1] <- 1 - predict.glm(mod2, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 0]
+    if (type == "first") {
+      tempdat$p.denominator[tempdat$exposure == 1 &
+                              tempdat$selvar == 1] <- predict.glm(mod2, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]
+    }
+    if (type == "cens") {
+      tempdat$p.denominator[tempdat$exposure == 1 &
+                              tempdat$selvar == 1] <- 1 - predict.glm(mod2, type = "response")[tempdat$exposure[tempdat$selvar == 1] == 1]
+    }
+    tempdat$p.denominator[tempdat$selvar == 0] <- 1
+    tempdat$w.denominator <- unlist(lapply(split(tempdat$p.denominator, tempdat$id), function(x)
+      cumprod(x)))
+    mod2$call$formula <- eval(parse(text = paste(
+      deparse(tempcall$exposure, width.cutoff = 500),
+      deparse(tempcall$denominator, width.cutoff = 500),
+      sep = ""
+    )))
+    mod2$call$family <- tempcall$link
+    mod2$call$data <- tempcall$data
+    mod2$call$subset <- paste(
+      "up to first instance of ",
+      deparse(tempcall$exposure, width.cutoff = 500),
+      " = 1 (selvar == 1)",
+      sep = ""
+    )
+    tempdat$ipw.weights <- tempdat$w.numerator / tempdat$w.denominator
+  }
+  #weights binomial, type "all"
+  if (tempcall$family == "binomial" & tempcall$type == "all") {
+    lf <- binomial(link = tempcall$link)
+    if (is.null(tempcall$numerator))
+      tempdat$w.numerator <- 1
+    else {
+      mod1 <- glm(
+        formula = eval(parse(text = paste(
+          deparse(tempcall$exposure, width.cutoff = 500),
+          deparse(tempcall$numerator, width.cutoff = 500),
+          sep = ""
+        ))),
+        family = lf,
+        data = data,
+        na.action = na.fail,
+        ...
+      )
+      tempdat$p.numerator <- vector("numeric", nrow(tempdat))
+      tempdat$p.numerator[tempdat$exposure == 0] <- 1 - predict.glm(mod1, type = "response")[tempdat$exposure == 0]
+      tempdat$p.numerator[tempdat$exposure == 1] <- predict.glm(mod1, type = "response")[tempdat$exposure == 1]
+      tempdat$w.numerator <- unlist(lapply(split(tempdat$p.numerator, tempdat$id), function(x)
+        cumprod(x)))
+      mod1$call$formula <- eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$numerator, width.cutoff = 500),
+        sep = ""
+      )))
+      mod1$call$family <- tempcall$link
+      mod1$call$data <- tempcall$data
+    }
+    mod2 <- glm(
+      formula = eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$denominator, width.cutoff = 500),
+        sep = ""
+      ))),
+      family = lf,
+      data = data,
+      na.action = na.fail,
+      ...
+    )
+    tempdat$p.denominator <- vector("numeric", nrow(tempdat))
+    tempdat$p.denominator[tempdat$exposure == 0] <- 1 - predict.glm(mod2, type = "response")[tempdat$exposure == 0]
+    tempdat$p.denominator[tempdat$exposure == 1] <- predict.glm(mod2, type = "response")[tempdat$exposure == 1]
+    tempdat$w.denominator <- unlist(lapply(split(tempdat$p.denominator, tempdat$id), function(x)
+      cumprod(x)))
+    mod2$call$formula <- eval(parse(text = paste(
+      deparse(tempcall$exposure, width.cutoff = 500),
+      deparse(tempcall$denominator, width.cutoff = 500),
+      sep = ""
+    )))
+    mod2$call$family <- tempcall$link
+    mod2$call$data <- tempcall$data
+    tempdat$ipw.weights <- tempdat$w.numerator / tempdat$w.denominator
+  }
+  #weights Cox
+  if (tempcall$family == "survival") {
+    if (is.null(tempcall$numerator))
+      tempdat$w.numerator <- 1
+    else {
+      mod1 <- coxph(
+        formula = eval(parse(
+          text = paste(
+            "Surv(",
+            deparse(tempcall$tstart),
+            ", ",
+            deparse(tempcall$timevar, width.cutoff = 500),
+            ", ",
+            deparse(tempcall$exposure, width.cutoff = 500),
+            ") ",
+            deparse(tempcall$numerator, width.cutoff = 500),
+            sep = ""
+          )
+        )),
+        data = data,
+        subset = tempdat$selvar == 1,
+        na.action = na.fail,
+        method = "efron",
+        ...
+      )
+      bh <- basehaz(mod1, centered = TRUE)
+      temp <- data.frame(timevar = sort(unique(tempdat$timevar)))
+      temp <- merge(
+        temp,
+        data.frame(
+          timevar = bh$time,
+          bashaz.cum.numerator = bh$hazard
+        ),
+        by = "timevar",
+        all.x = TRUE
+      )
+      rm(bh)
+      if (is.na(temp$bashaz.cum.numerator[temp$timevar == min(unique(tempdat$timevar))]))
+        temp$bashaz.cum.numerator[temp$timevar == min(unique(tempdat$timevar))] <- 0
+      temp$bashaz.cum.numerator <- approx(
+        x = temp$timevar,
+        y = temp$bashaz.cum.numerator,
+        xout = temp$timevar,
+        method = "constant",
+        rule = 2
+      )$y
+      temp$bashaz.numerator[1] <- temp$bashaz.cum.numerator[1]
+      temp$bashaz.numerator[2:nrow(temp)] <- diff(temp$bashaz.cum.numerator, 1)
+      temp$bashaz.cum.numerator <- NULL
+      tempdat <- merge(tempdat, temp, by = "timevar", all.x = TRUE)
+      rm(temp)
+      tempdat <- tempdat[order(tempdat$id, tempdat$timevar), ]
+      tempdat$risk.numerator[tempdat$selvar == 1] <- predict(mod1, type =
+                                                               "risk", centered = TRUE)
+      tempdat$hazard.numerator[tempdat$selvar == 1] <- with(tempdat[tempdat$selvar == 1, ], bashaz.numerator *
+                                                              risk.numerator)
+      tempdat$p.numerator[with(tempdat, selvar == 1 &
+                                 exposure == 0)] <- with(tempdat[with(tempdat, selvar == 1 &
+                                                                        exposure == 0), ], exp(-1 * bashaz.numerator * risk.numerator))
+      if (type == "first") {
+        tempdat$p.numerator[with(tempdat, selvar == 1 &
+                                   exposure == 1)] <- 1 - with(tempdat[with(tempdat, selvar == 1 &
+                                                                              exposure == 1), ], exp(-1 * bashaz.numerator * risk.numerator))
+      }
+      if (type == "cens") {
+        tempdat$p.numerator[with(tempdat, selvar == 1 &
+                                   exposure == 1)] <- with(tempdat[with(tempdat, selvar == 1 &
+                                                                          exposure == 1), ], exp(-1 * bashaz.numerator * risk.numerator))
+      }
+      tempdat$p.numerator[tempdat$selvar == 0] <- 1
+      tempdat$w.numerator <- unsplit(lapply(split(tempdat$p.numerator, tempdat$id), function(x)
+        cumprod(x)), tempdat$id)
+      mod1$call$formula <- eval(parse(
+        text = paste(
+          "Surv(",
+          deparse(tempcall$tstart),
+          ", ",
+          deparse(tempcall$timevar, width.cutoff = 500),
+          ", ",
+          deparse(tempcall$exposure, width.cutoff = 500),
+          ") ",
+          deparse(tempcall$numerator, width.cutoff = 500),
+          sep = ""
+        )
+      ))
+      mod1$call$data <- tempcall$data
+    }
+    mod2 <- coxph(
+      formula = eval(parse(
+        text = paste(
+          "Surv(",
+          deparse(tempcall$tstart),
+          ", ",
+          deparse(tempcall$timevar, width.cutoff = 500),
+          ", ",
+          deparse(tempcall$exposure, width.cutoff = 500),
+          ") ",
+          deparse(tempcall$denominator, width.cutoff = 500),
+          sep = ""
+        )
+      )),
+      data = data,
+      subset = tempdat$selvar == 1,
+      na.action = na.fail,
+      method = "efron",
+      ...
+    )
+    bh <- basehaz(mod2, centered = TRUE)
+    temp <- data.frame(timevar = sort(unique(tempdat$timevar)))
+    temp <- merge(
+      temp,
+      data.frame(
+        timevar = bh$time,
+        bashaz.cum.denominator = bh$hazard
+      ),
+      by = "timevar",
+      all.x = TRUE
+    )
+    rm(bh)
+    if (is.na(temp$bashaz.cum.denominator[temp$timevar == min(unique(tempdat$timevar))]))
+      temp$bashaz.cum.denominator[temp$timevar == min(unique(tempdat$timevar))] <- 0
+    temp$bashaz.cum.denominator <- approx(
+      x = temp$timevar,
+      y = temp$bashaz.cum.denominator,
+      xout = temp$timevar,
+      method = "constant",
+      rule = 2
+    )$y
+    temp$bashaz.denominator[1] <- temp$bashaz.cum.denominator[1]
+    temp$bashaz.denominator[2:nrow(temp)] <- diff(temp$bashaz.cum.denominator, 1)
+    temp$bashaz.cum.denominator <- NULL
+    tempdat <- merge(tempdat, temp, by = "timevar", all.x = TRUE)
+    rm(temp)
+    tempdat <- tempdat[order(tempdat$id, tempdat$timevar), ]
+    tempdat$risk.denominator[tempdat$selvar == 1] <- predict(mod2, type =
+                                                               "risk", centered = TRUE)
+    tempdat$hazard.denominator[tempdat$selvar == 1] <- with(tempdat[tempdat$selvar == 1, ], bashaz.denominator *
+                                                              risk.denominator)
+    tempdat$p.denominator[with(tempdat, selvar == 1 &
+                                 exposure == 0)] <- with(tempdat[with(tempdat, selvar == 1 &
+                                                                        exposure == 0), ], exp(-1 * bashaz.denominator * risk.denominator))
+    if (type == "first") {
+      tempdat$p.denominator[with(tempdat, selvar == 1 &
+                                   exposure == 1)] <- 1 - with(tempdat[with(tempdat, selvar == 1 &
+                                                                              exposure == 1), ], exp(-1 * bashaz.denominator * risk.denominator))
+    }
+    if (type == "cens") {
+      tempdat$p.denominator[with(tempdat, selvar == 1 &
+                                   exposure == 1)] <- with(tempdat[with(tempdat, selvar == 1 &
+                                                                          exposure == 1), ], exp(-1 * bashaz.denominator * risk.denominator))
+    }
+    tempdat$p.denominator[tempdat$selvar == 0] <- 1
+    tempdat$w.denominator <- unsplit(lapply(split(tempdat$p.denominator, tempdat$id), function(x)
+      cumprod(x)), tempdat$id)
+    mod2$call$formula <- eval(parse(
+      text = paste(
+        "Surv(",
+        deparse(tempcall$tstart),
+        ", ",
+        deparse(tempcall$timevar, width.cutoff = 500),
+        ", ",
+        deparse(tempcall$exposure, width.cutoff = 500),
+        ") ",
+        deparse(tempcall$denominator, width.cutoff = 500),
+        sep = ""
+      )
+    ))
+    mod2$call$data <- tempcall$data
+    mod2$call$subset <- paste(
+      "up to first instance of ",
+      deparse(tempcall$exposure, width.cutoff = 500),
+      " = 1 (selvar == 1)",
+      sep = ""
+    )
+    tempdat$ipw.weights <- tempdat$w.numerator / tempdat$w.denominator
+  }
+  #weights multinomial
+  if (tempcall$family == "multinomial") {
+    if (is.null(tempcall$numerator))
+      tempdat$p.numerator <- 1
+    else {
+      mod1 <- multinom(
+        formula = eval(parse(text = paste(
+          deparse(tempcall$exposure, width.cutoff = 500),
+          deparse(tempcall$numerator, width.cutoff = 500),
+          sep = ""
+        ))),
+        data = data,
+        subset = tempdat$selvar == 1,
+        na.action = na.fail,
+        ...
+      )
+      pred1 <- as.data.frame(predict(mod1, type = "probs"))
+      tempdat$p.numerator[tempdat$selvar == 0] <- 1
+      for (i in 1:length(unique(tempdat$exposure)))
+        tempdat$p.numerator[with(tempdat, tempdat$selvar == 1 &
+                                   exposure == sort(unique(tempdat$exposure))[i])] <- pred1[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i], i]
+      mod1$call$formula <- eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$numerator, width.cutoff = 500),
+        sep = ""
+      )))
+      mod1$call$data <- tempcall$data
+      mod1$call$subset <- paste(
+        "up to first instance of ",
+        deparse(tempcall$exposure, width.cutoff = 500),
+        " = 1 (selvar == 1)",
+        sep = ""
+      )
+    }
+    mod2 <- multinom(
+      formula = eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$denominator, width.cutoff = 500),
+        sep = ""
+      ))),
+      data = data,
+      subset = tempdat$selvar == 1,
+      na.action = na.fail,
+      ...
+    )
+    pred2 <- as.data.frame(predict(mod2, type = "probs"))
+    tempdat$p.denominator[tempdat$selvar == 0] <- 1
+    for (i in 1:length(unique(tempdat$exposure)))
+      tempdat$p.denominator[with(tempdat, tempdat$selvar == 1 &
+                                   exposure == sort(unique(tempdat$exposure))[i])] <- pred2[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i], i]
+    mod2$call$formula <- eval(parse(text = paste(
+      deparse(tempcall$exposure, width.cutoff = 500),
+      deparse(tempcall$denominator, width.cutoff = 500),
+      sep = ""
+    )))
+    mod2$call$data <- tempcall$data
+    mod2$call$subset <- paste(
+      "up to first instance of ",
+      deparse(tempcall$exposure, width.cutoff = 500),
+      " = 1 (selvar == 1)",
+      sep = ""
+    )
+    tempdat$ipw.weights <- unsplit(lapply(split(
+      with(tempdat, p.numerator / p.denominator), tempdat$id
+    ), function(x)
+      cumprod(x)), tempdat$id)
+  }
+  #weights ordinal
+  if (tempcall$family == "ordinal") {
+    if (tempcall$link == "logit")
+      m <- "logistic"
+    else
+      m <- tempcall$link
+    if (is.null(tempcall$numerator))
+      tempdat$p.numerator <- 1
+    else {
+      mod1 <- polr(
+        formula = eval(parse(
+          text = paste(
+            "as.factor(",
+            deparse(tempcall$exposure, width.cutoff = 500),
+            ")",
+            deparse(tempcall$numerator, width.cutoff = 500),
+            sep = ""
+          )
+        )),
+        data = data,
+        method = m,
+        subset = tempdat$selvar == 1,
+        na.action = na.fail,
+        ...
+      )
+      pred1 <- as.data.frame(predict(mod1, type = "probs"))
+      tempdat$p.numerator[tempdat$selvar == 0] <- 1
+      for (i in 1:length(unique(tempdat$exposure)))
+        tempdat$p.numerator[with(tempdat, tempdat$selvar == 1 &
+                                   exposure == sort(unique(tempdat$exposure))[i])] <- pred1[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i], i]
+      mod1$call$formula <- eval(parse(text = paste(
+        "as.factor(",
+        deparse(tempcall$exposure, width.cutoff = 500),
+        ")",
+        deparse(tempcall$numerator, width.cutoff = 500),
+        sep = ""
+      )))
+      mod1$call$data <- tempcall$data
+      mod1$call$method <- m
+      mod1$call$subset <- paste(
+        "up to first instance of ",
+        deparse(tempcall$exposure, width.cutoff = 500),
+        " = 1 (selvar == 1)",
+        sep = ""
+      )
+    }
+    mod2 <- polr(
+      formula = eval(parse(
+        text = paste(
+          "as.factor(",
+          deparse(tempcall$exposure, width.cutoff = 500),
+          ")",
+          deparse(tempcall$denominator, width.cutoff = 500),
+          sep = ""
+        )
+      )),
+      data = data,
+      method = m,
+      subset = tempdat$selvar == 1,
+      na.action = na.fail,
+      ...
+    )
+    pred2 <- as.data.frame(predict(mod2, type = "probs"))
+    tempdat$p.denominator[tempdat$selvar == 0] <- 1
+    for (i in 1:length(unique(tempdat$exposure)))
+      tempdat$p.denominator[with(tempdat, tempdat$selvar == 1 &
+                                   exposure == sort(unique(tempdat$exposure))[i])] <- pred2[tempdat$exposure[tempdat$selvar == 1] == sort(unique(tempdat$exposure))[i], i]
+    mod2$call$formula <- eval(parse(text = paste(
+      "as.factor(",
+      deparse(tempcall$exposure, width.cutoff = 500),
+      ")",
+      deparse(tempcall$denominator, width.cutoff = 500),
+      sep = ""
+    )))
+    mod2$call$data <- tempcall$data
+    mod2$call$method <- m
+    mod2$call$subset <- paste(
+      "up to first instance of ",
+      deparse(tempcall$exposure, width.cutoff = 500),
+      " = 1 (selvar == 1)",
+      sep = ""
+    )
+    tempdat$ipw.weights <- unsplit(lapply(split(
+      with(tempdat, p.numerator / p.denominator), tempdat$id
+    ), function(x)
+      cumprod(x)), tempdat$id)
+  }
+  #weights gaussian
+  if (tempcall$family == "gaussian") {
+    mod1 <- geeglm(
+      formula = eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$numerator, width.cutoff = 500),
+        sep = ""
+      ))),
+      data = data,
+      id = tempdat$id,
+      corstr = tempcall$corstr,
+      waves = tempdat$timevar,
+      ...
+    )
+    mod1$call$formula <- eval(parse(text = paste(
+      deparse(tempcall$exposure, width.cutoff = 500),
+      deparse(tempcall$numerator, width.cutoff = 500),
+      sep = ""
+    )))
+    mod1$call$data <- tempcall$data
+    mod1$call$id <- tempcall$id
+    mod1$call$corstr <- tempcall$corstr
+    mod1$call$waves <- tempcall$waves
+    mod2 <- geeglm(
+      formula = eval(parse(text = paste(
+        deparse(tempcall$exposure, width.cutoff = 500),
+        deparse(tempcall$denominator, width.cutoff = 500),
+        sep = ""
+      ))),
+      data = data,
+      id = tempdat$id,
+      corstr = tempcall$corstr,
+      waves = tempdat$timevar,
+      ...
+    )
+    mod2$call$formula <- eval(parse(text = paste(
+      deparse(tempcall$exposure, width.cutoff = 500),
+      deparse(tempcall$denominator, width.cutoff = 500),
+      sep = ""
+    )))
+    mod2$call$data <- tempcall$data
+    mod2$call$id <- tempcall$id
+    mod2$call$corstr <- tempcall$corstr
+    mod2$call$waves <- tempcall$waves
+    tempdat$kdens1 <- dnorm(tempdat$exposure, predict(mod1), as.numeric(sqrt(summary(mod1)$dispersion)[1]))
+    tempdat$kdens2 <- dnorm(tempdat$exposure, predict(mod2), as.numeric(sqrt(summary(mod2)$dispersion)[1]))
+    tempdat$ipw.weights <- unsplit(lapply(split(
+      with(tempdat, kdens1 / kdens2), tempdat$id
+    ), function(x)
+      cumprod(x)), tempdat$id)
+  }
+  #check for NA's in weights
+  if (sum(is.na(tempdat$ipw.weights)) > 0)
+    stop ("NA's in weights!")
+  
+  #truncate weights, when trunc value is specified (0-0.5)
+  if (!(is.null(tempcall$trunc))) {
+    tempdat$weights.trunc <- tempdat$ipw.weights
+    tempdat$weights.trunc[tempdat$ipw.weights <= quantile(tempdat$ipw.weights, 0 +
+                                                            trunc)] <- quantile(tempdat$ipw.weights, 0 + trunc)
+    tempdat$weights.trunc[tempdat$ipw.weights >  quantile(tempdat$ipw.weights, 1 -
+                                                            trunc)] <- quantile(tempdat$ipw.weights, 1 - trunc)
+  }
+  
+  #return results in the same order as the original input dataframe
+  if (is.null(tempcall$trunc)) {
+    if (is.null(tempcall$numerator))
+      return(
+        list(
+          ipw.weights = tempdat$ipw.weights[order(order.orig)],
+          call = tempcall,
+          selvar = tempdat$selvar[order(order.orig)],
+          den.mod = mod2
+        )
+      )
+    else
+      return(
+        list(
+          ipw.weights = tempdat$ipw.weights[order(order.orig)],
+          call = tempcall,
+          selvar = tempdat$selvar[order(order.orig)],
+          num.mod = mod1,
+          den.mod = mod2
+        )
+      )
+  }
+  else{
+    if (is.null(tempcall$numerator))
+      return(
+        list(
+          ipw.weights = tempdat$ipw.weights[order(order.orig)],
+          weights.trunc = tempdat$weights.trunc[order(order.orig)],
+          call = tempcall,
+          selvar = tempdat$selvar[order(order.orig)],
+          den.mod = mod2
+        )
+      )
+    else
+      return(
+        list(
+          ipw.weights = tempdat$ipw.weights[order(order.orig)],
+          weights.trunc = tempdat$weights.trunc[order(order.orig)],
+          call = tempcall,
+          selvar = tempdat$selvar[order(order.orig)],
+          num.mod = mod1,
+          den.mod = mod2
+        )
+      )
+  }
 }