逆向工程多项逻辑回归数据

Question

我正在研究多项逻辑回归问题（即，我想 class 确定名义结果变量的一些无序、独立水平）。我的问题是我知道结果变量的水平（在这个例子中，y=c('a','b','c')）并且我知道预测变量、它们的水平和它们的 class（这里，要么是 numeric/integer 要么标称）。我知道每个预测变量和结果之间的近似分布应该是什么（例如，x 的较高值与 y='a' 一起出现的频率更高，否则 x 的低值随机分布在另一个y 的水平）。

基本上，我想做 4 件事：1) 生成这些变量的数据集，这些变量近似于我指定的分布； 2) 运行数据的多项逻辑回归，nnet::multinom(y~.,df)； 3) 使用生成的模型predict()每个y级别的概率使用新数据；和 4) 检索进一步处理的概率。我对 MLR 模型的准确性或 p 值不感兴趣，所以我不需要将我的数据拆分为 train/test 个样本并进行 k 折交叉验证或任何事情。

我最初的想法是，这种基于某些用户指定分布的逆向工程数据集不能太不常见，并且可能 R package/function 可以做到这一点。到目前为止我还没有找到。到目前为止，我的方法如下：针对每个结果级别手动指定每个预测变量的每个级别的分布，如下所示：

rm(list=ls())
set.seed(123)

# specify vars and levels -- y=outcome var
y <- c('a','b','c')
x <- c(1:5)
p <- c(1:4)
r <- c(1:8)
q <- c('foo','bar','hello','world') # nominal var

# sample data based on user-specified distributions/probs
df1 <- data.frame(x1=sample(x,100,T,prob=c(0.1,0.1,0.2,0.25,0.35)),
                  y='b')
df2 <- data.frame(x1=sample(x,200,T,prob=c(0.35,0.25,0.2,0.1,0.1)),
                  y=sample(c('a','c'),200,T))
df <- rbind(df1,df2)

# check distribution of x1 levels v. y levels
table(df$x1,df$y)
     b  a  c
  1  7 38 30
  2 11 29 26
  3 22 17 22
  4 26 14  7
  5 34 12  5

问题在于，随着预测变量的数量越来越多且级别越来越多，这很乏味。我的下一个方法是生成随机数据样本，运行 MLR 模型，并调整模型权重。

# create random sample
df <- ldply(mget(ls()),
            function(x) sample(x,1000,T)) %>% 
  gather(k,v,-`.id`) %>%
  spread(`.id`,v) %>% select(-k)
str(df)
# change back vars to numeric
df[,c('p','r','x')] <- 
  apply(df[,c('p','r','x')],2,function(x) as.numeric(x))

glimpse(df)

Observations: 1,000
Variables: 5
$ p <dbl> 2, 2, 3, 1, 3, 2, 2, 4, 2, 4, 4, 3, 2, 4, 1, 4, 2, 1, 4, 3, 1, 3, 4, 3, 2, 2, 3...
$ q <chr> "bar", "bar", "foo", "bar", "world", "hello", "foo", "hello", "world", "hello",...
$ r <dbl> 2, 2, 1, 6, 6, 3, 4, 8, 6, 6, 2, 2, 8, 7, 7, 6, 3, 2, 4, 5, 2, 7, 1, 6, 3, 7, 8...
$ x <dbl> 2, 5, 1, 3, 3, 5, 2, 4, 1, 3, 5, 1, 5, 5, 2, 1, 1, 4, 4, 1, 5, 1, 5, 4, 4, 3, 2...
$ y <chr> "a", "c", "b", "a", "b", "a", "b", "c", "c", "b", "c", "c", "b", "a", "c", "b",...

# graph distribution of each predictor against each outcome -- not run here
# df %>% gather(k,v,-y) %>% group_by(y,k,v) %>%
#   summarise(n=n()) %>%
#   mutate(prop=n/sum(n)) %>%
#   ggplot(aes(y,prop,fill=v)) + geom_bar(stat='identity',position='dodge') +
#   facet_wrap(~k,scales='free') + theme(legend.position = 'none')

# run MLR model
m <- multinom(y~.,df)
summary(m)$coefficients
m$wts # coefficients from model

# adjust weight 16, which is x against y=b
m$wts[16] <- 1

同样，当预测变量和级别的数量很大时，这会很乏味。另外当我继续改变模型权重和预测新数据时，我得到了一些意想不到的概率（显然，我对 MLR 的了解不够自信，无法自信地使用这种方法）。

所以，我或多或少停留在这个阶段。我考虑过使用多重插补或自举方法来生成所需的数据，但我认为这两种方法都不适用。 MI 将估算 不完整案例 的数据，而我想指定有限数量的完整案例并从那里推断。同样，bootstrapping 将假设样本分布近似于总体分布，对数据重新采样。同样，我看不出如何指定有限数量的案例来有效地做到这一点（也许是引导加 permutation/shuffling ？）。

无论如何，任何 help/suggestions 在这里都非常感谢。并感谢任何真正阅读这篇冗长文章的人 post!

Answer 1

因此，我的解决方案是修改随机生成的数据，然后将修改后的数据（更接近所需的分布）用于运行 MLR 模型。

我创建了两个函数，一个重估数值变量，另一个重估名义变量。数值重估函数允许用户指定预测变量的值应该重新分配的方向，以及它们是否应该应用或排除结果变量的指定水平。下面的函数在问题中包含的示例数据上进行了测试。

当我返回运行 MLR 模型并预测新数据时，我得到每个结果的不同概率，更符合我的预期。

# redistribute values for specific predictors -----------------------------
# at specific levels of outcome var
####
# define function for numeric var
revalue.nums <- function(data,yvar.name,yvar.level,xvar.name,
                         direction=1, yvar.level.opposite=FALSE){
  # evaluate dir==-1 & oppo==T first, then dir==-1 & oppo==F,
  # then dir==1 & oppo==T, finally dir==1 & oppo==F
  if (direction==-1 & yvar.level.opposite==TRUE) {
    data[[xvar.name]][data[[yvar.name]] != yvar.level] <- 
      sample(get(xvar.name), 
             length(data[[xvar.name]][data[[yvar.name]] != yvar.level]), T,
             prob = c(seq(from=max(get(xvar.name)), 
                          to=min(get(xvar.name))) / sum(get(xvar.name))))
    return(data)
  } else if (direction==-1 & yvar.level.opposite==FALSE) {
    data[[xvar.name]][data[[yvar.name]]==yvar.level] <- 
      sample(get(xvar.name), 
             length(data[[xvar.name]][data[[yvar.name]]==yvar.level]), T,
             prob = c(seq(from=max(get(xvar.name)), 
                          to=min(get(xvar.name))) / sum(get(xvar.name))))
    return(data)
  } else if (direction==1 & yvar.level.opposite==TRUE) {
    data[[xvar.name]][data[[yvar.name]] != yvar.level] <- 
      sample(get(xvar.name), 
             length(data[[xvar.name]][data[[yvar.name]] != yvar.level]), T,
             prob = c(seq(from=min(get(xvar.name)), 
                          to=max(get(xvar.name))) / sum(get(xvar.name))))
    return(data)
  } else {
    data[[xvar.name]][data[[yvar.name]]==yvar.level] <- 
      sample(get(xvar.name), 
             length(data[[xvar.name]][data[[yvar.name]]==yvar.level]), T,
             prob = c(seq(from=min(get(xvar.name)), 
                          to=max(get(xvar.name))) / sum(get(xvar.name))))
    return(data)
  }
}
####

# define function
revalue.chars <- function(data,yvar.name,yvar.level,xvar.name,xvar.level,probs=0.25){
  data[[xvar.name]][data[[yvar.name]] == yvar.level] <- 
    sample(sort(sub(xvar.level,'1',get(xvar.name))),
           length(data[[xvar.name]][data[[yvar.name]] == yvar.level]), T,
           prob = c(probs, rep(probs / (length(get(xvar.name))-1),
                               rep(length(get(xvar.name))-1))))
  data[[xvar.name]][data[[xvar.name]] == '1'] <- xvar.level
  return(data)
}
###

# test functions on toy data
table(df$y,df$p) # orig
df1 <- revalue.nums(df,'y','a','p')
table(df1$y,df1$p) # changes y=a only, skew p to have higher values
df1 <- revalue.nums(df1,'y','a','p',yvar.level.opposite = T,direction = -1)
table(df1$y,df1$p) # changes y!=a, skew p to have lower values

table(df$y,df$q)
df2 <- revalue.chars(df,'y','b','q','hello',probs=0.5)
table(df2$y,df2$q) # increase num of q=hello and y=b occurrences

逆向工程多项逻辑回归数据

Reverse engineer multinomial logistic regression data

random

r

probability

nnet

multinomial