时间戳

首页 > 解决方案 > 逆向工程多项逻辑回归数据

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

问题描述

我正在研究多项逻辑回归问题（即，我想对标称结果变量的一些无序、独立水平进行分类）。我的问题是我知道结果变量的水平（在这个例子中，y=c('a','b','c')）并且我知道预测变量、它们的水平和它们的类（这里是数字/整数或名义）。我知道每个预测变量和结果之间的近似分布应该是什么（例如，较高的值x更频繁地出现在 中y='a'，否则较低的值x会随机分布在 的其他级别上y）。

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

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

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 将为不完整案例估算数据，而我想指定有限数量的完整案例并从那里推断。类似地，假设样本分布接近总体分布，自举将重新采样数据。同样，我不知道如何指定有限数量的情况下有效地做到这一点（也许是引导加上排列/改组？）。

无论如何，这里非常感谢任何帮助/建议。感谢任何真正阅读这篇冗长文章的人！

标签： rrandomprobabilitymultinomialnnet

解决方案

因此，我的解决方案是修改随机生成的数据，然后使用修改后的数据（更接近所需分布）来运行 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

推荐阅读