如何绘制由极坐标估计的模型的潜变量？

Question

我想知道是否有可能在 R （即 RStudio）中做一个类似于这个的情节：

我估计的模型是：

library(MASS)

# with logit
mod1 <- polr(lifesatisfaction) ~ gender + age + income + education + health + work less + work much), data = surveywave5, method = "logistic", Hess = TRUE) 

# with probit
mod1 <- polr(lifesatisfaction) ~ gender + age + income + education + health + work less + work much), data = surveywave5, method = "probit", Hess = TRUE) 


# my data:

 $ lifesatisfaction  : Ord.factor w/ 11 levels "0"<"1"<"2"<"3"<..: 9 9 10 10 10 9 11 10 11 7 ...        
 $ gender            : Factor w/ 2 levels "1","2": 2 1 1 1 1 1 2 1 2 1 ...
 $ income            : Factor w/ 10 levels "1","2","3","4",..: NA 2 4 5 5 10 7 7 6 3 ...
 $ age               : int  44 40 36 25 39 80 48 32 74 30 ...
 $ education         : Factor w/ 7 levels "1","2","3","4",..: 3 2 3 7 1 7 3 3 3 5 ...
 $ health            : Ord.factor w/ 5 levels "1","2","3","4",..: 3 4 1 3 4 5 5 4 4 3 ...
 $ work less         : Factor w/ 2 levels "0","1": 1 2 1 1 NA 1 1 1 2 1 ...
 $ work much         : Factor w/ 2 levels "0","1": 2 1 2 2 NA 1 2 2 1 2 ...

编辑*
我找到了这种方式..但是它似乎类似于 str().. 但不知道您是否可以将其用作可重现的：/

dput(head(surveywave5))
structure(list(gender = c(2, 1, 1, 1, 2, 2), maritalstatus = c(4, 6, NA, NA, 6, 6), age = c(62, 30, 44, 34, 58, 26), education = c(2, 7, 7, 7, 6, 4), lifesatisfaction = c(7, 8, 10, 7, 7, 8), health = c(4, 5, 5, 4, 5, 5), work.much = c(0, 1, 0, 0, 0, 0), work.less = c(1, 0, 1, 1, 1, 1), income = c(6, 1, 10, 6, 4, 1)), row.names = c(NA, -6L), class = c("tbl_df", "tbl", "data.frame"))  

＃＃＃编辑＃＃＃
每条曲线代表模型中使用的每个 x 变量，像这样

所以，一条曲线代表年龄，一条曲线代表性别、健康、收入等。

Answer 1

这是一个例子。首先，我们可以制作一些数据并运行模型

x <- rnorm(250)
z <- sample(0:1, 250, replace=TRUE)
b <- c(2, .75)

xb <- cbind(x,z) %*% b

tau <- c(-Inf, -3, 1, Inf)

q <- sapply(tau, (t)plogis(t-xb))
p <- q[,2:4]-q[,1:3]
y <- apply(p, 1, (x)which.max(rmultinom(1, 1, x)))

dat <- data.frame(x=x, z=z, y=as.factor(y))
mod <- MASS::polr(y ~ x + z, data=dat)

接下来，我们将创建两个数据集 - 一个将绘制线 (fake2)，另一个将包含我们将绘制分布的点 (fake)。在这些数据集中，您需要模型中的所有变量。您在x 上绘制的每一个都应该是恒定的（大概是某个中心值）。

fake <- data.frame(y = factor(1, levels=1:3), z=0, x=seq(-3,3,length=4))
fake2 <- data.frame(y = factor(1, levels=1:3), z=0, x=seq(-3,4.1,length=4))

接下来，我们可以使用两个假数据集制作线性预测器：

fxb <- model.matrix(formula(mod), data=fake)[,-1] %*% coef(mod)
fxb2 <- model.matrix(formula(mod), data=fake2)[,-1] %*% coef(mod)

现在，我们必须制作绘制分布的 y 值。最初这有点反复试验 - 试图获得正确的下限和上限。我做了这个，所以在最小值和第一个阈值 (mod$zeta[1]) 之间有 25 个值，在阈值之间有 25 个值：mod$zeta[1] 和 mod$zeta[2]，最后是从上限阈值到最大值的 25 个值。

s <- c(seq(-12, mod$zeta[1], length=25), 
       seq(mod$zeta[1], mod$zeta[2], length=25), 
       seq(mod$zeta[2], 12, length=25))

然后，我们就可以开始制作情节了。首先，我们绘制这条线：

par(mar = c(3,5.5,1,3))
plot(fake2$x, fxb2, type="l", ylim=c(-12, 12), xlim=c(-3,4.1), 
     axes=FALSE, xlab="", ylab = "")

接下来，我们需要知道要绘制分布的密度。我们可以从以线性预测器的值为中心的逻辑分布 PDF (dlogis()) 中得到这一点。由于这些是密度值，因此需要按比例放大它们以便分布可见。以下是第一个分布的密度值：

d1 <- dlogis(s, fxb[1])*5 + -3

现在，我们可以画出密度线。请注意，密度值在x 上，s 值在y 上：

lines(d1, s)

接下来，我们需要填写分布的各个部分。我们使用polygon() 函数来执行此操作。为此，我们需要指定一组完全连接的 x-y 对。我对x 使用d1 值，对y 使用s 值：

polygon(c(d1[1:25], rep(fake$x[1], 25), d1[1]), 
        y=c(s[1:25], s[25:1], s[1]), 
        col=rgb(1,0,0,.25), lty=0)

接下来，我们对分布的其他两个部分执行相同的操作：

polygon(c(d1[26:50], rep(fake$x[1], 25), d1[26]), 
        y=c(s[26:50], s[50:26], s[26]), 
        col=rgb(0,1,0,.25), lty=0)
polygon(c(d1[51:75], rep(fake$x[1], 25), d1[51]), 
        y=c(s[51:75], s[75:51], s[51]), 
        col=rgb(0,0,1,.25), lty=0)

然后，我们只需对要绘制分布的x 的其他值重复此操作：

d2 <- dlogis(s, fxb[2])*5 + -1
lines(d2, s)
polygon(c(d2[1:25], rep(fake$x[2], 25), d2[1]), 
        y=c(s[1:25], s[25:1], s[1]), 
        col=rgb(1,0,0,.25), lty=0)
polygon(c(d2[26:50], rep(fake$x[2], 25), d2[26]), 
        y=c(s[26:50], s[50:26], s[26]), 
        col=rgb(0,1,0,.25), lty=0)
polygon(c(d2[51:75], rep(fake$x[2], 25), d2[51]), 
        y=c(s[51:75], s[75:51], s[51]), 
        col=rgb(0,0,1,.25), lty=0)

d3 <- dlogis(s, fxb[3])*5 + 1
lines(d3, s)
polygon(c(d3[1:25], rep(fake$x[3], 25), d3[1]), 
        y=c(s[1:25], s[25:1], s[1]), 
        col=rgb(1,0,0,.25), lty=0)
polygon(c(d3[26:50], rep(fake$x[3], 25), d3[26]), 
        y=c(s[26:50], s[50:26], s[26]), 
        col=rgb(0,1,0,.25), lty=0)
polygon(c(d3[51:75], rep(fake$x[3], 25), d3[51]), 
        y=c(s[51:75], s[75:51], s[51]), 
        col=rgb(0,0,1,.25), lty=0)

d4 <- dlogis(s, fxb[4])*5 + 3
lines(d4, s)
polygon(c(d4[1:25], rep(fake$x[4], 25), d4[1]), 
        y=c(s[1:25], s[25:1], s[1]), 
        col=rgb(1,0,0,.25), lty=0)
polygon(c(d4[26:50], rep(fake$x[4], 25), d4[26]), 
        y=c(s[26:50], s[50:26], s[26]), 
        col=rgb(0,1,0,.25), lty=0)
polygon(c(d4[51:75], rep(fake$x[4], 25), d4[51]), 
        y=c(s[51:75], s[75:51], s[51]), 
        col=rgb(0,0,1,.25), lty=0)
abline(h=mod$zeta, lty=2)
axis(1)
axis(2, at=c(-7.5, -.9, 6.55), labels=c("Outcome 1", "Outcome 2", "Outcome 3"), las=1)
mtext("Latent Propensity", 4, line=1)
mtext("Values of X", 1, line=2)
box()

^{由reprex package (v2.0.1) 于 2022 年 10 月 27 日创建}

我怀疑有一种更优雅的方法可以做到这一点，但我不确定它是什么。