R/DGP_func.R

In saudiwin/idealstan: Robust Measurement with 'Stan'

# DGP Functions -----------------------------------------------------------

# Functions used to generate all the models
#' @noRd
.binary <- function(pr_absence=NULL,
                    pr_vote=NULL,
                    y=NULL,
                    N=NULL,
                    inflate=NULL,
                    person_points=NULL,
                    item_points=NULL,
                    time_points=NULL,
                    type='simulate',
                    outcome=NULL,
                    latent_space=NULL,
                    cov_effect=NULL,
                    person_x=NULL,
                    ...) {

  
  # need a factor to multiply the probability for latent-space models
  if(latent_space && inflate) {
    mul_fac <- 2
  } else {
    mul_fac <- 1
  }

  #standard IRT 2-PL model
  if(type=='simulate') {
    votes <- as.numeric((mul_fac*pr_vote)>runif(N))
  } else if(type=='predict') {
    votes <- t(apply(pr_vote,1,function(c) as.numeric((c*mul_fac)>runif(N))))
  } else if(type=="epred") {
    
    votes <- pr_vote
    
  } else if(type=='log_lik') {
    if(inflate) {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- ifelse(is.na(outcome), 
                          dbinom(1,size = 1,prob=pr_absence[,c]*mul_fac,log=T),
                          dbinom(outcome-1,size=1,prob=pr_vote[,c]*mul_fac,log=T))
      })
      
    } else {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- dbinom(outcome-1,size=1,prob=pr_vote[,c]*mul_fac,log=T)})
    }
    return(t(over_iters))
  }
  
  
  # remove pr of inflate if model is not inflated
  
  if(!inflate) {
    pr_boost <- 1
  } else {
    pr_boost <- 0
  }
  
  if(type=='simulate') {
    combined <- ifelse((pr_absence*mul_fac)<(runif(N)+pr_boost),votes,NA)
    
    combined <- combined - min(combined,na.rm=T)
    
    # Create a score dataset
    
    out_data <- data_frame(outcome_disc=combined,
                           person_id=person_points,
                           time_id=time_points,
                           item_id=item_points,
                           group_id='G')
    
        if(!is.null(cov_effect) && cov_effect!=0) {
      
      cov_data <- distinct(out_data,person_id,time_id) %>% 
        arrange(person_id, time_id) %>% 
        mutate(person_x=person_x)
      
      out_data <- left_join(out_data,
                            cov_data,by=c("person_id","time_id"))
      
      out_data <- id_make(score_data=out_data,person_cov = ~ person_x)
      
    } else {
      
      out_data <- id_make(score_data=out_data)
      
    }
    
    return(out_data) 
  } else if(type=='predict') {
    combined <- sapply(1:nrow(pr_absence), function(c) ifelse(pr_absence[c,]<(runif(N)+pr_boost),votes[c,],2))
    # add one to have minimum = 1
    combined <- combined + 1
    attr(combined,'output') <- 'all'
    # transpose to make S x N matrix
    return(t(combined))
    
  } else if(type=="epred") {
    
    # need to calculate pr(Yes)pr(Present) + Pr(Yes)pr(Absent)
    
    if(inflate) {
      
      combined <- (1 - pr_absence) * pr_vote

    } else {
      
      combined <- pr_vote
      
    }

    attr(combined,'output') <- 'all'
    # transpose to make S x N matrix
    return(combined)
    
  }
                   
}

#' @importFrom ordbetareg rordbeta dordbeta
.ordbeta <- function(pr_absence=NULL,
                                 pr_vote=NULL,
                                 N=NULL,
                                 inflate=NULL,
                                 person_points=NULL,
                                 item_points=NULL,
                                 time_points=NULL,
                                 cutpoints=NULL,
                                 ordinal_outcomes=NULL,
                                 type='simulate',
                                output=NULL,
                                 y=NULL,
                                 phi=NULL,
                                 outcome=NULL,
                                 cov_effect=NULL,
                                 person_x=NULL,
                                 miss_val=NULL,
                                 ...)
{
  
  if(!inflate) {
    pr_boost <- 1
  } else {
    pr_boost <- 0
  }
  
  if(type=='simulate') {
    
    cutpoints <- quantile(pr_vote,probs=seq(0,1,length.out = 3))
    cutpoints <- cutpoints[c(1,3)]
    
  } else if(type=='predict') {
    # over posterior draws
    
    cutpoints <- filter(cutpoints, item_id==item_points[1]) %>% 
      ungroup %>% 
      spread(key="cut",value=".value") %>% 
      select(`1`,`2`) %>% 
      as.matrix
    # 
    # cuts_iters <- sapply(1:nrow(cutpoints), function(i) {
    #   cuts <- sapply(1:ncol(cutpoints),function(y) {
    #     qlogis(pr_vote[i,]) - cutpoints[i,y]
    #   })
    # },simplify='array')
    
  } else if(type=='log_lik') {
    if(inflate) {
      n_outcomes <- length(unique(outcome)) - 1
    } else {
      n_outcomes <- length(unique(outcome))
    }
    # over posterior draws
    cuts_iters <- sapply(1:nrow(cutpoints), function(i) {
      cuts <- sapply(1:ncol(cutpoints),function(y) {
        qlogis(pr_vote[i,]) - c(cutpoints[i,y])
      })
    },simplify='array')
    
  }
  
  # Now we pick votes as a function of the number of categories
  # This code should work for any number of categories
  
  if(type=='simulate') {
    
    votes <- rordbeta(n=length(pr_vote),
                      mu=pr_vote,
                      phi=phi,
                      cutpoints=qlogis(cutpoints))
    
    combined <- ifelse(pr_absence<(runif(N)+pr_boost),votes,NA)
    
    # Create a score dataset
    
    out_data <- data_frame(outcome_cont=combined,
                           person_id=person_points,
                           time_id=time_points,
                           item_id=item_points,
                           group_id='G')
    
    if(!is.null(cov_effect) && cov_effect!=0) {
      
      cov_data <- distinct(out_data,person_id,time_id) %>% 
        arrange(person_id, time_id) %>% 
        mutate(person_x=person_x)
      
      out_data <- left_join(out_data,
                            cov_data,by=c("person_id","time_id"))
      
      out_data <- id_make(score_data=out_data,person_cov = ~ person_x)
      
    } else {
      
      out_data <- id_make(score_data=out_data)
      
    }
    
    return(out_data) 
    
  } else if(type=='predict') {
    
    if(output=='observed') {
      
      combined <- sapply(1:nrow(pr_vote), function(i) {
        
        rordbeta(n=ncol(pr_vote),mu=pr_vote[i,],phi=phi[i],
                 cutpoints=cutpoints[i,])
        
              })
      
      combined <- t(combined)
      
      attr(combined,'output') <- 'observed'
      attr(combined,'output_type') <- 'continuous'
    } else if(output=='missing') {
      
      combined <- apply(pr_absence, 2,function(c) as.numeric(c>runif(N)))
      
      attr(combined,'output') <- 'missing'
      attr(combined,'output_type') <- 'discrete'
    }
    
    return(combined)
    
  } else if(type=='log_lik') {
    
    out_num <- as.numeric(outcome)
    
    over_iters <- sapply(1:ncol(pr_vote), function(d) {
        
        dordbeta(x=out_num,
                 mu=pr_vote[,d],
                 phi=phi[d],
                 cutpoints=cuts_iters[,d])

    },simplify='array')
    
    if(inflate) {
      # remove top category for vote prediction
      out_num[out_num==max(out_num)] <- max(out_num) - 1
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- ifelse(is.na(outcome), 
                          dbinom(1,size = 1,prob=pr_absence[,c],log=T),
                          log(over_iters[out_num,,c]))
      })
    } else {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- log(over_iters[out_num,,c])
      })
    }
    
    attr(outdens,'output') <- 'all'
    return(t(outdens))
  }
  
}

.ordinal_ratingscale <- function(pr_absence=NULL,
                    pr_vote=NULL,
                    N=NULL,
                    inflate=NULL,
                    person_points=NULL,
                    item_points=NULL,
                    time_points=NULL,
                    cutpoints=NULL,
                    ordinal_outcomes=NULL,
                    type='simulate',
                    y=NULL,
                    outcome=NULL,
                    cov_effect=NULL,
                    person_x=NULL,
                    miss_val=NULL,
                    ...)
{
  
  # if(inflate && type!='simulate') {
  #   ordinal_outcomes <- ordinal_outcomes -1
  # }
  
  if(!inflate) {
    pr_boost <- 1
  } else {
    pr_boost <- 0
  }

  if(type=='simulate') {
    cutpoints <- quantile(pr_vote,probs=seq(0,1,length.out = ordinal_outcomes+1))
    cutpoints <- cutpoints[2:(length(cutpoints)-1)]
    
    #Generate outcomes by personlator
    
    cuts <- sapply(cutpoints,function(y) {
      qlogis(pr_vote) - qlogis(y)
    },simplify='array')
  } else if(type=='predict') {
    # over posterior draws
    cuts_iters <- sapply(1:nrow(cutpoints), function(i) {
      cuts <- sapply(1:ncol(cutpoints),function(y) {
        qlogis(pr_vote[i,]) - c(cutpoints[i,y])
      })
    },simplify='array')

  } else if(type=='log_lik') {
    if(inflate) {
      n_outcomes <- length(unique(outcome)) - 1
    } else {
      n_outcomes <- length(unique(outcome))
    }
    # over posterior draws
    cuts_iters <- sapply(1:nrow(cutpoints), function(i) {
      cuts <- sapply(1:ncol(cutpoints),function(y) {
        qlogis(pr_vote[i,]) - c(cutpoints[i,y])
      })
    },simplify='array')
    
  }
  
  # Now we pick votes as a function of the number of categories
  # This code should work for any number of categories
  
  if(type=='simulate') {
    votes <- sapply(1:nrow(cuts), function(i) {
      this_cut <- cuts[i,]
      
      pr_bottom <- 1 - plogis(this_cut[1])
      
      mid_prs <- sapply(1:(length(this_cut)-1), function(c) {
        plogis(this_cut[c]) - plogis(this_cut[c+1])
      })
      
      pr_top <- plogis(this_cut[length(this_cut)])
      
      return(sample(1:(length(this_cut)+1),size=1,prob=c(pr_bottom,mid_prs,pr_top)))
    })
    
    combined <- ifelse(pr_absence<(runif(N)+pr_boost),votes,NA)
    
    # Create a score dataset
    
    out_data <- data_frame(outcome_disc=combined,
                           person_id=person_points,
                           time_id=time_points,
                           item_id=item_points,
                           ordered_id=ordinal_outcomes,
                           group_id='G')
    
    if(!is.null(cov_effect) && cov_effect!=0) {
      
      cov_data <- distinct(out_data,person_id,time_id) %>% 
        arrange(person_id, time_id) %>% 
        mutate(person_x=person_x)
      
      out_data <- left_join(out_data,
                            cov_data,by=c("person_id","time_id"))
      
      out_data <- id_make(score_data=out_data,person_cov = ~ person_x)
      
    } else {
      
      out_data <- id_make(score_data=out_data)
      
    }
    
    return(out_data) 
  } else if(type=='predict') {

    over_iters <- sapply(1:nrow(pr_vote), function(d) {
      votes <- sapply(1:dim(cuts_iters)[1], function(i) {
        
        this_cut <- cuts_iters[i,,d]
        
        pr_bottom <- 1 - plogis(this_cut[1])
        
        mid_prs <- sapply(1:(length(this_cut)-1), function(c) {
          plogis(this_cut[c]) - plogis(this_cut[c+1])
        })
        
        pr_top <- plogis(this_cut[length(this_cut)])
        
        return(sample(as.integer(ordinal_outcomes),
                      size=1,prob=c(pr_bottom,mid_prs,pr_top)))
      })
    })

    combined <- sapply(1:ncol(pr_absence), function(c) ifelse(pr_absence[,c]<(runif(N)+pr_boost),over_iters[c,],miss_val))
    
    attr(combined,'output') <- 'all'
    return(combined)
  } else if(type=='log_lik') {
    over_iters <- sapply(1:ncol(pr_vote), function(d) {
      votes <- sapply(1:dim(cuts_iters)[1], function(i) {
        
        this_cut <- cuts_iters[i,,d]
        
        pr_bottom <- 1 - plogis(this_cut[1])
        
        mid_prs <- sapply(1:(length(this_cut)-1), function(c) {
          plogis(this_cut[c]) - plogis(this_cut[c+1])
        })
        
        pr_top <- plogis(this_cut[length(this_cut)])
        
        return(c(pr_bottom,mid_prs,pr_top))
      })
    },simplify='array')
    out_num <- as.numeric(outcome)

    if(inflate) {
      # remove top category for vote prediction
      out_num[out_num==max(out_num)] <- max(out_num) - 1
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- ifelse(is.na(outcome), 
                          dbinom(1,size = 1,prob=pr_absence[,c],log=T),
                          log(over_iters[out_num,,c]))
      })
    } else {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- log(over_iters[out_num,,c])
      })
    }
    
    attr(outdens,'output') <- 'all'
    return(t(outdens))
  }
                   
}

.ordinal_grm <- function(pr_absence=NULL,
                                 pr_vote=NULL,
                         y=NULL,
                                 N=NULL,
                                 inflate=NULL,
                                 person_points=NULL,
                                 item_points=NULL,
                                 time_points=NULL,
                                 ordinal_outcomes=NULL,
                         cutpoints=NULL,
                         type='simulate',
                         outcome=NULL,
                         cov_effect=NULL,
                         person_x=NULL,
                         miss_val=NULL,
                                 ...)
{
  
  if(inflate && type!='simulate') {
    ordinal_outcomes <- ordinal_outcomes -1
  }
  
  if(!inflate) {
    pr_boost <- 1
  } else {
    pr_boost <- 0
  }

  # need one set of cutpoints for each item
  if(type=='simulate') {
    all_cuts <- sapply(1:max(item_points), function(i) {
      cutpoints <- sort(runif(2))
    })
    all_cuts <- all_cuts[,item_points]
    
    #Generate outcomes by person and item
    
    cuts <- sapply(1:(ordinal_outcomes-1),function(y) {
      qlogis(pr_vote) - qlogis(all_cuts[y,])
    },simplify='array')
  } else {
    # over posterior draws
    
    cuts_iters <- sapply(1:dim(cutpoints)[1], function(i) {
        qlogis(pr_vote[i,]) - cutpoints[i,item_points]
    },simplify='array')
  }

  
  if(type=='simulate') {
    # Now we pick votes as a function of the number of categories
    # This code should work for any number of categories
    votes <- sapply(1:nrow(cuts), function(i) {
      this_cut <- cuts[i,]
      
      pr_bottom <- 1 - plogis(this_cut[1])
      
      mid_prs <- sapply(1:(length(this_cut)-1), function(c) {
        plogis(this_cut[c]) - plogis(this_cut[c+1])
      })
      
      pr_top <- plogis(this_cut[length(this_cut)])
      
      return(sample(1:(length(this_cut)+1),size=1,prob=c(pr_bottom,mid_prs,pr_top)))
    })
    
    # remove pr of absence if model is not inflated
    
    combined <- ifelse(pr_absence<(runif(N)+pr_boost),votes,NA)
    
    # Create a score dataset
    
    out_data <- data_frame(outcome_disc=combined,
                           person_id=person_points,
                           ordered_id=ordinal_outcomes,
                           time_id=time_points,
                           item_id=item_points,
                           group_id='G')
    
    if(!is.null(cov_effect) && cov_effect!=0) {
      
      cov_data <- distinct(out_data,person_id,time_id) %>% 
        arrange(person_id, time_id) %>% 
        mutate(person_x=person_x)
      
      out_data <- left_join(out_data,
                            cov_data,by=c("person_id","time_id"))
      
      out_data <- id_make(score_data=out_data,person_cov = ~ person_x)
      
    } else {
      
      out_data <- id_make(score_data=out_data)
      
    }
    
    return(out_data)                      
  } else if(type=='predict') {
    over_iters <- sapply(1:nrow(pr_vote), function(d) {
      votes <- sapply(1:dim(cuts_iters)[1], function(i) {
        
        this_cut <- cuts_iters[i,d]
        
        pr_bottom <- 1 - plogis(this_cut[1])
        
        mid_prs <- sapply(1:(length(this_cut)-1), function(c) {
          plogis(this_cut[c]) - plogis(this_cut[c+1])
        })
        
        pr_top <- plogis(this_cut[length(this_cut)])
        
        return(sample(1:(length(this_cut)+1),size=1,prob=c(pr_bottom,mid_prs,pr_top)))
      })
    })
    
    combined <- sapply(1:ncol(pr_absence), function(c) ifelse(pr_absence[,c]<(runif(N)+pr_boost),over_iters[c,],miss_val))
    attr(combined,'output') <- 'all'
    return(t(combined))
  } else if(type=='log_lik') {
    
    over_iters <- sapply(1:ncol(pr_vote), function(d) {
      votes <- sapply(1:dim(cuts_iters)[1], function(i) {
        
        this_cut <- cuts_iters[i,,d]
        
        pr_bottom <- 1 - plogis(this_cut[1])
        
        mid_prs <- sapply(1:(length(this_cut)-1), function(c) {
          plogis(this_cut[c]) - plogis(this_cut[c+1])
        })
        
        pr_top <- plogis(this_cut[length(this_cut)])
        
        return(c(pr_bottom,mid_prs,pr_top))
      })
    },simplify='array')
    
    out_num <- as.numeric(outcome)

    if(inflate) {
      # remove top category for vote prediction
      out_num[out_num==max(out_num)] <- max(out_num) - 1
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- ifelse(is.na(outcome), 
                          dbinom(1,size = 1,prob=pr_absence[,c],log=T),
                          log(over_iters[out_num,,c]))
      })
    } else {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- log(over_iters[out_num,,c])
      })
    }
    
    attr(outdens,'output') <- 'all'
    return(t(outdens))
  }
  
}

.poisson <- function(pr_absence=NULL,
                    pr_vote=NULL,
                    y=NULL,
                    N=NULL,
                    max_val=NULL,
                    inflate=NULL,
                    person_points=NULL,
                    item_points=NULL,
                    time_points=NULL,
                    type='simulate',
                    output=NULL,
                    outcome=NULL,
                    cov_effect=NULL,
                    person_x=NULL,
                    ...)
{

  #standard IRT 2-PL model
  if(type=='simulate') {
    votes <- rpois(n = length(pr_vote),lambda = exp(qlogis(pr_vote)))
  } else if(type=='predict') {
    votes <- apply(pr_vote,2,function(c) rpois(n=length(c),lambda = exp(qlogis(c))))
  } else if(type=='log_lik') {
    if(inflate) {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- ifelse(is.na(outcome), 
                          dbinom(1,size = 1,prob=pr_absence[,c],log=T),
                          dpois(outcome,lambda=exp(qlogis(pr_vote[,c])),log=T))
      })
      
    } else {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- dpois(outcome,lambda=exp(qlogis(pr_vote[,c])),log=T)})
    }
    return(t(over_iters))
  }
  
  
  # remove pr of inflate if model is not inflated
  
  if(!inflate) {
    pr_boost <- 1
  } else {
    pr_boost <- 0
  }
  
  if(type=='simulate') {

    combined <- ifelse(pr_absence<(runif(N)+pr_boost),votes,NA)
    
    # Create a score dataset
    
    out_data <- data_frame(outcome_disc=combined,
                           person_id=person_points,
                           time_id=time_points,
                           item_id=item_points,
                           group_id='G')
    
    if(!is.null(cov_effect) && cov_effect!=0) {
      
      cov_data <- distinct(out_data,person_id,time_id) %>% 
        arrange(person_id, time_id) %>% 
        mutate(person_x=person_x)
      
      out_data <- left_join(out_data,
                            cov_data,by=c("person_id","time_id"))
      
      out_data <- id_make(score_data=out_data,person_cov = ~ person_x,
                          unbounded=T)
      
    } else {
      
      out_data <- id_make(score_data=out_data,
                          unbounded=T)
      
    }
    
    
    
    return(out_data) 
  } else if(type=='predict') {
    if(output=='observed') {
      combined <- votes
      attr(combined,'output') <- 'observed'
      attr(combined,'output_type') <- 'discrete'
    } else if(output=='missing') {
      combined <- apply(pr_absence, 2,function(c) as.numeric(c>runif(N)))
      attr(combined,'output') <- 'missing'
      attr(combined,'output_type') <- 'discrete'
    }
    # transpose to make S x N matrix
    return(combined)
  } 
  
}

.normal <- function(pr_absence=NULL,
                     pr_vote=NULL,
                    y=NULL,
                     N=NULL,
                    max_val=NULL,
                     inflate=NULL,
                     person_points=NULL,
                     item_points=NULL,
                     time_points=NULL,
                    sigma_sd=NULL,
                    type='simulate',
                    output='observed',
                    outcome=NULL,
                    cov_effect=NULL,
                    person_x=NULL,
                     ...)
{

  #standard IRT 2-PL model
  if(type=='simulate') {
    votes <- rnorm(n = length(pr_vote),mean = qlogis(pr_vote),sd = sigma_sd)
  } else if(type=='predict') {

    votes <- sapply(1:ncol(pr_vote),function(c) rnorm(n=nrow(pr_vote),mean=qlogis(pr_vote[,c]),sd=sigma_sd[c]))
    
  } else if(type=='log_lik') {
    if(inflate) {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- ifelse(is.na(outcome), 
                          dbinom(1,size = 1,prob=pr_absence[,c],log=T),
                          dnorm(outcome,mean=qlogis(pr_vote[,c]),sd=sigma_sd[c],log=T))
      })
      
    } else {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- dnorm(outcome,mean=qlogis(pr_vote[,c]),sd=sigma_sd[c],log=T)})
    }
    return(t(over_iters))
  }
  
  
  # remove pr of inflate if model is not inflated
  
  if(!inflate) {
    pr_boost <- 1
  } else {
    pr_boost <- 0
  }
  
  if(type=='simulate') {
    
    combined <- ifelse(pr_absence<(runif(N)+pr_boost),votes,NA)
    
    # Create a score dataset
    
    out_data <- data_frame(outcome_cont=combined,
                           person_id=person_points,
                           time_id=time_points,
                           item_id=item_points,
                           group_id='G')
    
    if(!is.null(cov_effect) && cov_effect!=0) {
      
      cov_data <- distinct(out_data,person_id,time_id) %>% 
        arrange(person_id, time_id) %>% 
        mutate(person_x=person_x)
      
      out_data <- left_join(out_data,
                            cov_data,by=c("person_id","time_id"))
      
      out_data <- id_make(score_data=out_data,person_cov = ~ person_x)
      
    } else {
      
      out_data <- id_make(score_data=out_data)
      
    }
    
    return(out_data) 
  } else if(type=='predict') {
    if(output=='observed') {
      combined <- votes
      attr(combined,'output') <- 'observed'
      attr(combined,'output_type') <- 'continuous'
    } else if(output=='missing') {
      combined <- apply(pr_absence, 2,function(c) as.numeric(c>runif(N)))
      attr(combined,'output') <- 'missing'
      attr(combined,'output_type') <- 'discrete'
    }
    # transpose to make S x N matrix
    return(combined)
  } 
  
}

.lognormal <- function(pr_absence=NULL,
                    pr_vote=NULL,
                    N=NULL,
                    y=NULL,
                    max_val=NULL,
                    inflate=NULL,
                    person_points=NULL,
                    item_points=NULL,
                    time_points=NULL,
                    sigma_sd=NULL,
                    type='simulate',
                    output='observed',
                    outcome=NULL,
                    cov_effect=NULL,
                    person_x=NULL,
                    ...)
{
  
  #standard IRT 2-PL model
  if(type=='simulate') {
    votes <- rlnorm(n = length(pr_vote),meanlog = qlogis(pr_vote),sdlog = sigma_sd)
  } else if(type=='predict') {
    votes <- sapply(1:ncol(pr_vote),function(c) rlnorm(n=nrow(pr_vote),meanlog=qlogis(pr_vote[,c]),sdlog=sigma_sd[c]))
  } else if(type=='log_lik') {
    if(inflate) {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- ifelse(is.na(outcomel), 
                          dbinom(1,size = 1,prob=pr_absence[,c],log=T),
                          dlnorm(outcome,meanlog=qlogis(pr_vote[,c]),sdlog=sigma_sd[c],log=T))
      })

    } else {
      over_iters <- sapply(1:ncol(pr_vote), function(c) {
        outdens <- dlnorm(outcome,meanlog=qlogis(pr_vote[,c]),sdlog=sigma_sd[c],log=T)})
    }
    return(t(over_iters))
  }
  
  
  # remove pr of inflate if model is not inflated
  
  if(!inflate) {
    pr_boost <- 1
  } else {
    pr_boost <- 0
  }
  
  if(type=='simulate') {
    
    combined <- ifelse(pr_absence<(runif(N)+pr_boost),votes,NA)
    
    # Create a score dataset
    
    out_data <- data_frame(outcome_cont=combined,
                           person_id=person_points,
                           time_id=time_points,
                           item_id=item_points,
                           group_id='G')
    
    if(!is.null(cov_effect) && cov_effect!=0) {
      
      cov_data <- distinct(out_data,person_id,time_id) %>% 
        arrange(person_id, time_id) %>% 
        mutate(person_x=person_x)
      
      out_data <- left_join(out_data,
                            cov_data,by=c("person_id","time_id"))
      
      out_data <- id_make(score_data=out_data,person_cov = ~ person_x,
                          unbounded=T)
      
    } else {
      
      out_data <- id_make(score_data=out_data,
                          unbounded=T)
      
    }
    
    return(out_data) 
  } else if(type=='predict') {
    if(output=='observed') {
      combined <- votes
      attr(combined,'output') <- 'observed'
      attr(combined,'output_type') <- 'continuous'
    } else if(output=='missing') {
      combined <- apply(pr_absence, 2,function(c) as.numeric(c>runif(N)))
      attr(combined,'output') <- 'missing'
      attr(combined,'output_type') <- 'discrete'
    }
    # transpose to make S x N matrix
    return(combined)
  }               
}

#' Function to generate random-walk or AR(1) person parameters
#' Recursively generate data
#' @noRd
.gen_ts_data <- function(t,adj_in,alpha_int,sigma,init_sides) {
  current_val <- new.env()
  current_val$t1 <- 0
  
  out_vec <- lapply(1:t,function(t_1) {
    
    if(t_1==1) {
      t_11 <- alpha_int
      current_val$t1 <- t_11
      return(data_frame(t_11))
    } else {
      if(adj_in==1) {
        t_11 <- adj_in*current_val$t1 + rnorm(n=1,sd=sigma)
      } else {
        t_11 <- alpha_int + adj_in*current_val$t1 + rnorm(n=1,sd=sigma)
      }
      
    }
    current_val$t1 <- t_11
    return(data_frame(t_11))
  })  %>% bind_rows
  return(out_vec)
}