全部版块 我的主页
论坛 数据科学与人工智能 数据分析与数据科学 R语言论坛
23774 15
2010-06-29
我通过R软件实现有序结果的logistic回归,程序如下:
library(nnet)
library(MASS)
library(epicalc)
C=rep(c(1,2,3,4),4)
M=c(0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3)
F=c(22,81,30,3,57,236,135,26,11,112,105,17,1,4,10,7)
C<- ordered(C)
m<- polr(C~M,weights=F)
m
ordinal.or.display(m)
单单输入m,是可以显示部分结果的,但是运行ordinal.or.display(m)命令时,R软件报错:
错误于model.frame.default(formula = C ~ M, weights = F) :  变数的长度不一样('(weights)')
我感到十分奇怪,敬请各位大侠指教
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

全部回复
2010-6-29 21:14:37
library(nnet)
library(MASS)
library(epicalc)
C=rep(c(1,2,3,4),4)
M=c(0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3)
F=c(22,81,30,3,57,236,135,26,11,112,105,17,1,4,10,7)
mydata=data.frame(cbind(C,M,F))
mydata$C=ordered(mydata$C)
mydata$M=factor(mydata$M)
m<- polr(C~M,weights=F,data=mydata)
m
ordinal.or.display(m)
summary(m)
> summary(m)

Re-fitting to get Hessian

Call:
polr(formula = C ~ M, data = mydata, weights = F)

Coefficients:
       Value Std.      Error         t value
M1 0.4797267  0.1899885 2.525031
M2 1.1012361  0.2071994 5.314862
M3 2.5307873  0.4596561 5.505828

Intercepts:
        Value   Std. Error t value
1|2 -1.5741  0.1801    -8.7411
2|3  1.1020  0.1725     6.3874
3|4  3.4674  0.2196    15.7922

Residual Deviance: 1869.092
AIC: 1881.092
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-6-29 22:01:24
结果和我原来的程序不一样啊
有序结果的logistic回归,回归系数应该只有一个啊
library(nnet)
library(MASS)
library(epicalc)
C=rep(c(1,2,3,4),4)
M=c(0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3)
F=c(22,81,30,3,57,236,135,26,11,112,105,17,1,4,10,7)
C<- ordered(C)
m<- polr(C~M,weights=F)
m
library(nnet)
library(MASS)
library(epicalc)
C=rep(c(1,2,3,4),4)
M=c(0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3)
F=c(22,81,30,3,57,236,135,26,11,112,105,17,1,4,10,7)
C<- ordered(C)
m<- polr(C~M,weights=F)
m
Call:
polr(formula = C ~ M, weights = F)

Coefficients:
        M
0.6373681

Intercepts:
      1|2       2|3       3|4
-1.457851  1.225450  3.563591

Residual Deviance: 1873.039
AIC: 1881.039
敬请指教,谢谢!
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-6-30 20:14:17
这是因为function ordinal.or.display()
有用到function confint()
function confint() return A matrix (or vector)
由于你只用一个predicator "M"
所以有错误产生
若你改为 2 predictors
m<- polr(as.ordered(C)~M+F,data=mydata)
则OK

当然稍修改程序也可运行1 predictor
  Ordinal OR lower95ci upper95ci P value
M 1.891        1.576         2.274         8.38e-12
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-7-1 09:06:38
非常感谢您的解答,我还想请教一下
我这里C代表儿童智力等级,M代表母亲文化程度,F代表频数
通过有序logistic回归,来分析母亲文化程度对儿童智力的影响,因此只有一个predictor
您说稍微修改程序,也可得到如下结果,不知应该如何修改
Ordinal OR lower95ci upper95ci P value
M 1.891        1.576         2.274         8.38e-12
敬请指教,不胜感激!!
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-7-1 11:10:15
library(MASS)
library(epicalc)
C=rep(c(1,2,3,4),4)
M=c(0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3)
F=c(22,81,30,3,57,236,135,26,11,112,105,17,1,4,10,7)
mydata=data.frame(cbind(C,M,F))
m<- polr(as.ordered(C)~M,weights=F,data=mydata)
m
summary(m)
####ordinal.or.display revised
ordinal.or.display1 <- function (ordinal.model, decimal = 3, alpha = 0.05)
{
    model <- ordinal.model
    if (class(model) != "polr")
        stop("The model is not an ordinal logistic regression model")
    s1 <- summary(model)
    t <- s1$coefficients[, 3]
    df <- s1$df.residual
    p.value <- pt(abs(t), df, lower.tail = FALSE)
    coeff <- t(t(model$coefficients))
    coeff.95ci <- cbind(coeff, matrix(confint(model, level = 1 - alpha),dim(coeff)[1],2))
    oor.95ci <- round(exp(coeff.95ci), decimal)
    len.p <- length(p.value)
    k=length(model$zeta)
    oor.95ci <- cbind(oor.95ci, format(p.value[-(len.p:(len.p -(k-1)))], digits = decimal))
    colnames(oor.95ci) <- c("Ordinal OR", paste("lower", 100 -
        100 * alpha, "ci", sep = ""), paste("upper", 100 - 100 *
        alpha, "ci", sep = ""), "P value")
    print.noquote(oor.95ci)
}
ordinal.or.display1(m)
>
  Ordinal OR lower95ci upper95ci P value
M 1.891        1.576         2.274        8.38e-12

##########compare ordinal.or.display() & ordinal.or.display1()
nes96 <- read.table("http://www.stat.washington.edu/quinn/classes/536/data/nes96r.dat", header=TRUE)
polr.out <- polr(as.ordered(ClinLR)~TVnews+selfLR+PID+age+educ+income,data=nes96)
ordinal.or.display(polr.out)
>
Re-fitting to get Hessian
Waiting for profiling to be done...
Re-fitting to get Hessian

       Ordinal OR lower95ci upper95ci P value
TVnews 0.994      0.947     1.043     3.99e-01
selfLR   1.074      0.961     1.201     1.05e-01
PID       0.68        0.631     0.731     3.93e-24
age       0.992      0.984     1            2.25e-02
educ     0.884      0.816     0.959      1.53e-03
income 0.969      0.949     0.99        2.26e-03
Warning message:
In cbind(oor.95ci, format(p.value[-(len.p:(len.p - 1))], digits = decimal)) :
  number of rows of result is not a multiple of vector length (arg 2)

##########
ordinal.or.display1(polr.out)
>
Re-fitting to get Hessian
Waiting for profiling to be done...
Re-fitting to get Hessian

       Ordinal OR lower95ci upper95ci P value
TVnews 0.994      0.947     1.043     3.99e-01
selfLR 1.074         0.961     1.201     1.05e-01
PID      0.68          0.631      0.731     3.93e-24
age      0.992        0.984      1            2.25e-02
educ    0.884        0.816      0.959     1.53e-03
income 0.969       0.949      0.99       2.26e-03
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

点击查看更多内容…
相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群