R语言使用二元回归将序数数据建模为多元GLM

用于分析序数数据的最常见模型是 逻辑模型 。本质上，您将结果视为连续潜在变量的分类表现。此结果的预测变量仅以一种方式对其产生影响，因此 为每个预测变量获得一个回归系数。但是该模型有几个截距，它们代表将变量切分以创建观察到的分类表现的点。

就像在普通回归模型中一样，每个预测变量都会以一种方式影响结果，这就是比例赔率假设或约束。或者，可以让每个预测变量在每个切入点对结果产生不同的影响。

如何使用单变量GLM软件对此建模？UCLA idre页面上有关于多元随机系数模型的文章。在这里很重要，因为他们使用nlme（单变量线性混合模型软件）对多元结果进行建模。基本思想是将数据堆叠起来，使其成为一种重复测量，但是找到一种向软件发出信号的信号，即结果是不同的，从而对预测变量要求不同的截距和斜率。

因此，我们要做的是将数据从宽转换为长，将其建模为常规二项式，但是我们需要告诉模型为每个级别估计不同的截距。为此，我使用具有unstructured工作相关性结构的通用估计方程（GEE）。3

示范

library(ordinal) # For ordinal regression to check our results
library(geepack) # For GEE with binary data

数据集。

soup <- ordinal::soup
soup$ID <- 1:nrow(soup) # Create a person ID variable
str(soup)

‘data.frame‘:	1847 obs. of  13 variables:
 $ RESP    : Factor w/ 185 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ PROD    : Factor w/ 2 levels "Ref","Test": 1 2 1 2 1 2 2 2 2 1 ...
 $ PRODID  : Factor w/ 6 levels "1","2","3","4",..: 1 2 1 3 1 6 2 4 5 1 ...
 $ SURENESS: Ord.factor w/ 6 levels "1"<"2"<"3"<"4"<..: 6 5 5 6 5 5 2 5 5 2 ...
 $ DAY     : Factor w/ 2 levels "1","2": 1 1 1 1 2 2 2 2 2 2 ...
 $ SOUPTYPE: Factor w/ 3 levels "Self-made","Canned",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ SOUPFREQ: Factor w/ 3 levels ">1/week","1-4/month",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ COLD    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
 $ EASY    : Factor w/ 10 levels "1","2","3","4",..: 7 7 7 7 7 7 7 7 7 7 ...
 $ GENDER  : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 2 2 2 2 2 ...
 $ AGEGROUP: Factor w/ 4 levels "18-30","31-40",..: 4 4 4 4 4 4 4 4 4 4 ...
 $ LOCATION: Factor w/ 3 levels "Region 1","Region 2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ ID      : int  1 2 3 4 5 6 7 8 9 10 ...

我使用SURENESS变量。它有6个级别。使用DAY和GENDER变量对其进行建模。

# Select variables to work with
soup <- dplyr::select(soup, ID, SURENESS, DAY, GENDER)
# I like dummy variables with recognizable names
soup$girl <- ifelse(soup$GENDER == "Female", 1, 0) # Make male reference group
soup$day2 <- ifelse(soup$DAY == "2", 1, 0) # Make day 1 reference group

下一步是将顺序结果转换为代表每个阈值的5个结果

完成此操作后，我们准备对这5个新的结果变量进行转换。


head(soup.long) # Let‘s look at the data

     ID SURENESS DAY GENDER girl day2 SURE VAL SURE.f
1     1        6   1 Female    1    0    2   1      2
1848  1        6   1 Female    1    0    3   1      3
3695  1        6   1 Female    1    0    4   1      4
5542  1        6   1 Female    1    0    5   1      5
7389  1        6   1 Female    1    0    6   1      6
2     2        5   1 Female    1    0    2   1      2

让我们看看没有选择最高响应类别的人：



     ID SURENESS DAY GENDER girl day2 SURE VAL SURE.f
22   22        4   1 Female    1    0    2   1      2
1869 22        4   1 Female    1    0    3   1      3
3716 22        4   1 Female    1    0    4   1      4
5563 22        4   1 Female    1    0    5   0      5
7410 22        4   1 Female    1    0    6   0      6

这个人选择了SURENESSVAL中的a 。她的前三个分数是1，她的最后两个分数是0，因为4小于4-5阈值和5-6阈值。

下一步是为阈值创建虚拟变量。这些变量将用于表示模型中的截距。

请注意，我将虚拟变量乘以-1。在序数回归中，这样做使解释更容易。总之，它确保正系数增加了从较低类别（例如3）移至较高类别（4）或对较高响应类别做出响应的几率。

现在，我们准备运行模型。我们使用GEE。相关结构为unstructured。

接下来，我使用标准序数回归估算模型：

让我们比较系数和标准误差：

        Estimate Estimate.1 Std.err Std. Error     Wald  z value Pr(>|W|) Pr(>|z|)
SURE.f2 -2.13244   -2.13155 0.10454    0.10450 416.0946 -20.3971   0.0000   0.0000
SURE.f3 -1.19345   -1.19259 0.09142    0.09232 170.4284 -12.9179   0.0000   0.0000
SURE.f4 -0.89164   -0.89079 0.08979    0.09011  98.5995  -9.8857   0.0000   0.0000
SURE.f5 -0.65782   -0.65697 0.08945    0.08898  54.0791  -7.3833   0.0000   0.0000
SURE.f6 -0.04558   -0.04477 0.08801    0.08789   0.2682  -0.5093   0.6046   0.6105
girl    -0.04932   -0.04917 0.09036    0.09074   0.2980  -0.5419   0.5851   0.5879
day2    -0.26172   -0.26037 0.08584    0.08579   9.2954  -3.0351   0.0023   0.0024

可以看到结果非常接近。

但是，使用估计glm()不能建立一个人的结果之间的依存关系的估计会产生不同的结果。


        Estimate Std. Error z value Pr(>|z|)
SURE.f2 -2.15144    0.08255 -26.062   0.0000
SURE.f3 -1.21271    0.06736 -18.004   0.0000
SURE.f4 -0.91149    0.06472 -14.084   0.0000
SURE.f5 -0.67782    0.06327 -10.713   0.0000
SURE.f6 -0.06523    0.06178  -1.056   0.2911
girl    -0.07326    0.04961  -1.477   0.1398
day2    -0.26898    0.04653  -5.780   0.0000

估计值和标准误均不足。

我们可以轻松地放宽pom.bin模型中的比例赔率约束。让我们通过放宽对预测变量的约束来运行某些人所说的偏比例赔率模型day2。我们通过估计阈值虚拟变量和day2预测变量之间的相互作用来做到这一点。

我还使用名义参数运行了相同的模型进行比较day2。



             Estimate Estimate.1 Std.err Std. Error     Wald  z value Pr(>|W|) Pr(>|z|)
SURE.f2      -2.02982   -2.03106 0.11800    0.11834 295.8986 -17.1630  0.00000  0.00000
SURE.f3      -1.22087   -1.22213 0.09829    0.09857 154.2801 -12.3980  0.00000  0.00000
SURE.f4      -0.92773   -0.92899 0.09458    0.09443  96.2112  -9.8375  0.00000  0.00000
SURE.f5      -0.65744   -0.65870 0.09246    0.09188  50.5554  -7.1693  0.00000  0.00000
SURE.f6      -0.04733   -0.04859 0.08955    0.08965   0.2793  -0.5420  0.59714  0.58784
SURE.f2:day2  0.07359    0.07360 0.14148    0.14155   0.2705   0.5199  0.60298  0.60312
SURE.f3:day2  0.31691    0.31697 0.10607    0.10613   8.9270   2.9867  0.00281  0.00282
SURE.f4:day2  0.33301    0.33308 0.09970    0.09973  11.1551   3.3398  0.00084  0.00084
SURE.f5:day2  0.26330    0.26339 0.09618    0.09616   7.4938   2.7391  0.00619  0.00616
SURE.f6:day2  0.26741    0.26748 0.09347    0.09345   8.1842   2.8622  0.00423  0.00421
girl         -0.04809   -0.04994 0.09048    0.09077   0.2825  -0.5502  0.59507  0.58221

结果是可比较的。

现在，我们可以将比例比例赔率二进制模型与比例赔率二进制模型进行比较，以测试day2变量的约束条件。geepack允许anova()对两种模型进行Wald测试：


Analysis of ‘Wald statistic‘ Table

Model 1 VAL ~ 0 + SURE.f2 + SURE.f3 + SURE.f4 + SURE.f5 + SURE.f6 + girl + SURE.f2:day2 + SURE.f3:day2 + SURE.f4:day2 + SURE.f5:day2 + SURE.f6:day2
Model 2 VAL ~ 0 + SURE.f2 + SURE.f3 + SURE.f4 + SURE.f5 + SURE.f6 + girl + day2
  Df   X2 P(>|Chi|)
1  4 6.94      0.14

两种模型之间的差异在统计上均不显着，表明day2变量的比例约束已足够。

我们可以使用或使用函数ordinal进行比较pom.ord和npom.ord建模anova()，从而进行相同的测试nomimal_test()。两者都是似然比检验，比上述GEE的Wald检验更充分。



Likelihood ratio tests of cumulative link models:

         formula:               nominal: link: threshold:
pom.ord  SURENESS ~ girl + day2 ~1       logit flexible
npom.ord SURENESS ~ girl        ~day2    logit flexible  

         no.par  AIC logLik LR.stat df Pr(>Chisq)
pom.ord       7 5554  -2770
npom.ord     11 5555  -2766    6.91  4       0.14

nominal_test(pom.ord)

Tests of nominal effects

formula: SURENESS ~ girl + day2
       Df logLik  AIC  LRT Pr(>Chi)
<none>     -2770 5554
girl    4  -2766 5554 8.02    0.091 .
day2    4  -2766 5555 6.91    0.141
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这两个测试收敛到相同的结果，并且在比较GEE模型的Wald测试中也给出了相同的p值。然而，Wald- χ 2χ2 测试统计数据略高。


完成此操作后，使用序数数据包当然要容易得多。但是，将模型视为二进制可能会有一些好处，但是所有这些都是出于好奇而非必要。由于某种原因，我仍未弄清楚，当一个人尝试使用fitted()函数从模型中获得预测的概率时，它仅返回一组拟合的概率。理想情况下，它应该为每个阈值返回拟合概率。使用geepack，可以直接获得每个级别的预测概率。但是，这种优势是微不足道的。


而且，如果熟悉最大似然估计，则可以简单地对似然函数进行编程。

上面的例子在比例赔率情况下的语法为：


coef(summary(res))
      Estimate Std. Error
a1 -2.13155603 0.10450286
a2 -1.19259266 0.09232077
a3 -0.89079068 0.09010891
a4 -0.65697671 0.08898063
a5 -0.04477565 0.08788869
bg -0.04917604 0.09073602
bd -0.26037369 0.08578617

coef(summary(pom.ord))
        Estimate Std. Error     z value     Pr(>|z|)
1|2  -2.13155281 0.10450291 -20.3970663 1.775532e-92
2|3  -1.19259171 0.09232091 -12.9178937 3.567748e-38
3|4  -0.89078590 0.09010896  -9.8856524 4.804418e-23
4|5  -0.65697465 0.08898068  -7.3833401 1.543671e-13
5|6  -0.04476553 0.08788871  -0.5093434 6.105115e-01
girl -0.04917245 0.09073601  -0.5419287 5.878676e-01
day2 -0.26037360 0.08578617  -3.0351465 2.404188e-03

结果非常相似，对于比较模型的更确定的方式，我们总是可以比较对数似然：

logLik(res)
‘log Lik.‘ -2769.784 (df=7)

logLik(pom.ord)
‘log Lik.‘ -2769.784 (df=7)
  1. Agresti，A。（2013）。分类数据分析。Wiley-Interscience。 ?

如果您有任何疑问，请在下面发表评论。

大数据部落 -中国专业的第三方数据服务提供商，提供定制化的一站式数据挖掘和统计分析咨询服务

统计分析和数据挖掘咨询服务：y0.cn/teradat（咨询服务请联系官网客服）

?QQ：3025393450

?QQ交流群：186388004 

【服务场景】  

科研项目; 公司项目外包;线上线下一对一培训;数据爬虫采集;学术研究;报告撰写;市场调查。

【大数据部落】提供定制化的一站式数据挖掘和统计分析咨询

欢迎关注微信公众号，了解更多数据干货资讯！

欢迎选修我们的R语言数据分析挖掘必知必会课程！

原文地址：https://www.cnblogs.com/tecdat/p/12228536.html

时间： 2024-08-29 21:52:35

R语言使用二元回归将序数数据建模为多元GLM的相关文章

用R语言进行分位数回归：基础篇

用R语言进行分位数回归:基础篇 詹鹏 (北京师范大学经济管理学院 北京) http://www.xiaowanxue.com/up_files/2012121819040.html 原文地址:https://www.cnblogs.com/jwg-fendi/p/10069488.html

R语言基础入门之二：数据导入和描述统计

by 写长城的诗 • October 30, 2011 • Comments Off This post was kindly contributed by 数据科学与R语言 - go there to comment and to read  the full post. 一.数据导入 对初学者来讲,面对一片空白的命令行窗口,第一道真正的难关也许就是数据的导入.数据导入有很多途径,例如从网页抓取.公共数据源获得.文本文件导入.为了快速入门,建议初学者采取R语言协同Excel电子表格的方法.也就

回归预测及R语言实现 Part2 回归R语言实现

下面是回归分析的各种变体的简单介绍,解释变量和相应变量就是指自变量和因变量. 常用普通最小二乘(OLS)回归法来拟合实现简单线性.多项式和多元线性等回归模型.最小二乘法的基本原理前面已经说明了,使得预测值和观察值之差最小. R中实现拟合线性模型最基本的函数是lm(),应用格式为: myfit <- lm(Y~X1+X2+-+Xk,data) data为观测数据,应该为一个data.frame,前面是拟合表达式,Y是因变量,X1-Xk是自变量,+用来分隔不同的自变量的,还有可能用到的其他符号的说明

用R语言分析我的fitbit计步数据

目标:把fitbit的每日运动记录导入到R语言中进行分析,画出统计图表来 已有原始数据:fitbit2014年每日的记录电子表格文件,示例如下: 日期 消耗卡路里数 步 距离 攀爬楼层数 久坐不动的分钟数 不太活跃分钟数 中度活跃分钟数 非常活跃分钟数 2014年4月27日 2736 16581 11.84 7 1111 131 117 81 2014年4月28日 2514 12622 9.01 6 910 136 59 76 2014年4月29日 2231 8357 5.97 9 1208 1

R语言与mysql结合处理交通数据及其算法优化

一.序言 交通数据处理是智能交通的一个很关键的要素,更好的分析交通数据,可以为市政管理.交通信号管制.道路规划.交通设施建设提供更好的咨询和建议.全国各地政府都在寄期望于智能交通,以缓解城市拥堵,甚至一定程度上解决大城市病或者说是市政建设滞后的问题.同时,诸如百度地图.谷歌地图.高德地图.微软地图都推出了相应的交通应用,以期找到更大的商机. 用好的存储方法和好的算法进行分析,在批处理方面可以更多的分析历史数据,分析和发现问题,为未来进行预测以及公共查询服务:在实时计算方面可以更多的进行交通监控.

回归预测及R语言实现 Part1 回归基础综述

Part1 回归基础综述 回归方法有很多种,最常见的是线性回归(又有一元和多元之分).多项式回归.非线性回归.另外还将简单说明对预测结果的检验方法.   线性回归 一元线性回归,是最简单最常见的回归模型,类似初中数学中的一元一次方程,它的基本模型如下: 我们常见的一元线性回归方程一般没有最后一项,确切的说,我们在实际的应用中也忽略了最后一项.最后一项ui的现实意义是:它是指除自变量x以外所有对因变量y有影响的其他因素,应用回归预测时,我们假设ui是一个均值为零的随机变量,方差为常值,不同ui间相

R语言之逻辑回归

本文主要将逻辑回归的实现,模型的检验等 参考博文http://blog.csdn.net/tiaaaaa/article/details/58116346;http://blog.csdn.net/ai_vivi/article/details/43836641 1.测试集和训练集(3:7比例)数据来源:http://archive.ics.uci.edu/ml/datasets/statlog+(australian+credit+approval) austra=read.table("au

R语言入门视频笔记--4--R的数据输入

R的数据输入可以大体三种: 1.键盘输出 2.从文本文件导入 3.从Excel中导入数据 一.从键盘输入 首先创建一个数据框,玩玩嘛,瞎建一个 mydata <- data.frame(age =numeric(0),gender= character(0),weight=numeric(0))    #建一个空数据框,但已经声明过元素类型 mydata <- edit(mydata)                        #可以进行编辑 fix(mydata) #跟上面一样可以进行编

R语言做逻辑回归

前面写过一个多分类的逻辑回归,现在要做一个简单的二分类,用glm函数 导入csv格式如下: mydata<-read.csv("D://li.csv",header=T) colnames(mydata)<-c("x1","x2","x3","y") model<-glm(formula = y ~ x1+x2+x3, family = quasibinomial(link = "