如何使用 mgcv 包从非线性模型在 R 中创建 ggplot?
How do I create a ggplot in R from a non-linear model using the mgcv package?
我有一个非线性生存模型,我使用 mgcv
包对其进行了编码。我可以生成一个常规图,但我希望能够编写 ggplot2
代码。我该怎么做?
这是我的代码:
df <- structure(list(SurvYear =c(3L, 2L, 3L, 6L, 8L, 3L, 5L, 2L, 9L,
8L, 1L, 7L, 1L, 4L, 6L, 8L, 2L, 5L, 1L, 1L, 7L, 1L, 5L, 3L, 2L,
1L, 9L, 1L, 5L, 2L, 2L, 1L, 2L, 3L, 4L, 8L, 7L, 2L, 2L, 6L, 9L,
7L, 3L, 9L, 6L, 8L, 2L, 8L, 2L, 1L, 1L, 6L, 5L, 3L, 3L, 7L, 2L,
4L, 5L, 2L, 3L, 7L, 4L, 1L, 2L, 2L, 3L, 5L, 1L, 9L, 2L, 2L, 3L,
9L, 6L, 2L, 2L, 4L, 3L, 1L, 9L, 7L, 3L, 1L, 2L, 1L, 6L, 3L, 1L,
5L, 6L, 5L, 6L, 4L, 2L, 1L, 3L, 1L, 1L, 3L, 4L, 3L, 8L, 9L, 7L,
6L, 3L, 5L, 2L, 7L, 9L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 9L, 1L,
4L, 8L, 1L, 8L, 1L, 1L, 8L, 5L, 2L, 9L, 4L, 8L, 4L, 9L, 2L, 2L,
3L, 2L, 9L, 3L, 2L, 1L, 3L, 2L, 1L, 9L, 9L, 2L, 1L, 1L, 1L, 2L,
9L, 1L, 5L, 1L, 6L, 9L, 3L, 2L, 2L, 5L, 7L, 4L, 2L, 7L, 2L, 4L,
5L, 3L, 3L, 9L, 2L, 6L, 1L, 3L, 4L, 5L, 9L, 8L, 1L, 2L, 8L, 2L,
9L, 1L, 7L, 3L, 3L, 1L, 6L, 3L, 4L, 9L, 1L, 3L, 4L, 4L, 2L, 7L,
2L, 3L, 1L, 1L, 7L, 2L, 1L, 1L, 2L, 1L, 9L, 1L, 2L, 9L, 1L, 1L,
2L, 3L, 7L, 3L, 1L, 1L, 2L, 5L, 4L, 6L, 7L, 1L, 9L, 2L, 1L, 8L,
1L, 2L, 1L, 4L, 2L, 3L, 3L, 9L, 9L, 9L, 4L, 1L, 1L, 4L, 9L, 3L,
1L, 1L, 3L, 3L, 4L, 1L, 1L, 1L, 1L, 6L, 9L, 1L, 1L, 8L, 1L, 3L,
3L, 8L, 3L, 5L, 1L, 2L, 1L, 2L, 4L, 3L, 1L, 6L, 1L, 4L, 8L, 1L,
3L, 2L, 2L, 3L, 6L, 2L, 1L, 1L, 1L, 9L, 3L, 1L, 7L, 3L, 9L, 1L,
9L, 5L, 4L), Gender = c(1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L,
1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L,
1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L,
0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L,
1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L,
1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L,
0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L,
0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L,
0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L,
1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L,
1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L,
1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L,
1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L,
1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L,
0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L,
1L, 1L), Age = c(63L, 66L, 34L, 43L, 63L, 21L, 24L, 44L, 52L,
59L, 27L, 32L, 30L, 20L, 56L, 55L, 35L, 26L, 53L, 43L, 39L, 19L,
34L, 28L, 19L, 24L, 50L, 22L, 58L, 24L, 50L, 25L, 37L, 30L, 51L,
69L, 23L, 49L, 22L, 46L, 58L, 31L, 23L, 53L, 59L, 25L, 38L, 44L,
34L, 49L, 19L, 39L, 24L, 51L, 29L, 27L, 48L, 77L, 22L, 43L, 59L,
49L, 60L, 51L, 49L, 47L, 50L, 44L, 41L, 44L, 50L, 42L, 46L, 54L,
35L, 21L, 26L, 26L, 40L, 21L, 48L, 49L, 20L, 20L, 32L, 37L, 22L,
36L, 46L, 28L, 39L, 35L, 51L, 39L, 49L, 57L, 46L, 18L, 52L, 47L,
27L, 32L, 23L, 43L, 42L, 57L, 22L, 40L, 19L, 58L, 71L, 55L, 42L,
20L, 51L, 21L, 20L, 61L, 36L, 54L, 19L, 35L, 38L, 41L, 34L, 22L,
41L, 42L, 56L, 50L, 53L, 53L, 48L, 22L, 59L, 27L, 28L, 32L, 37L,
68L, 24L, 26L, 61L, 21L, 20L, 20L, 50L, 62L, 61L, 29L, 18L, 40L,
67L, 43L, 25L, 43L, 22L, 56L, 47L, 41L, 40L, 43L, 27L, 37L, 61L,
35L, 23L, 54L, 38L, 38L, 39L, 45L, 49L, 63L, 49L, 44L, 44L, 23L,
37L, 58L, 61L, 25L, 18L, 59L, 25L, 51L, 40L, 27L, 42L, 22L, 38L,
22L, 45L, 33L, 32L, 36L, 53L, 52L, 19L, 45L, 53L, 27L, 65L, 25L,
53L, 57L, 29L, 23L, 62L, 36L, 56L, 59L, 41L, 61L, 44L, 24L, 21L,
38L, 29L, 55L, 33L, 18L, 21L, 19L, 65L, 24L, 59L, 34L, 25L, 45L,
48L, 18L, 41L, 61L, 32L, 37L, 21L, 20L, 57L, 25L, 65L, 50L, 61L,
32L, 27L, 19L, 50L, 63L, 19L, 45L, 20L, 36L, 20L, 19L, 53L, 39L,
50L, 20L, 24L, 57L, 28L, 21L, 39L, 49L, 21L, 20L, 39L, 20L, 44L,
19L, 39L, 53L, 29L, 60L, 43L, 21L, 23L, 30L, 42L, 42L, 51L, 35L,
50L, 51L, 56L, 52L, 22L, 36L, 56L, 28L, 57L, 20L, 47L, 48L, 65L,
71L, 21L, 70L, 23L, 63L), Highest_Educationmx = c(4L, 5L, 3L,
2L, 3L, 2L, 3L, 1L, 3L, 1L, 7L, 3L, 2L, 3L, 3L, 2L, 6L, 2L, 3L,
6L, 3L, 2L, 2L, 7L, 2L, 1L, 2L, 3L, 6L, 3L, 5L, 3L, 5L, 6L, 2L,
1L, 5L, 2L, 5L, 1L, 1L, 3L, 2L, 3L, 1L, 7L, 5L, 4L, 7L, 3L, 1L,
1L, 6L, 3L, 3L, 2L, 4L, 6L, 5L, 4L, 2L, 6L, 1L, 3L, 4L, 2L, 1L,
5L, 5L, 3L, 1L, 5L, 3L, 3L, 1L, 4L, 2L, 3L, 5L, 3L, 1L, 4L, 2L,
1L, 2L, 7L, 2L, 5L, 3L, 2L, 6L, 1L, 1L, 3L, 4L, 1L, 5L, 1L, 3L,
4L, 2L, 7L, 2L, 4L, 4L, 7L, 4L, 6L, 3L, 1L, 2L, 1L, 5L, 5L, 1L,
5L, 2L, 7L, 3L, 4L, 2L, 4L, 2L, 4L, 2L, 2L, 4L, 1L, 2L, 1L, 2L,
6L, 1L, 2L, 5L, 2L, 2L, 5L, 1L, 6L, 5L, 2L, 1L, 2L, 1L, 1L, 3L,
2L, 4L, 3L, 2L, 3L, 1L, 5L, 5L, 7L, 1L, 3L, 3L, 2L, 1L, 3L, 4L,
5L, 1L, 1L, 3L, 3L, 3L, 5L, 3L, 6L, 4L, 3L, 1L, 3L, 5L, 7L, 1L,
3L, 4L, 5L, 3L, 3L, 1L, 1L, 1L, 7L, 3L, 1L, 4L, 3L, 3L, 5L, 1L,
4L, 5L, 4L, 2L, 5L, 3L, 1L, 1L, 5L, 4L, 7L, 5L, 2L, 2L, 5L, 3L,
1L, 1L, 2L, 3L, 5L, 3L, 7L, 5L, 1L, 5L, 3L, 1L, 1L, 1L, 1L, 7L,
5L, 7L, 3L, 1L, 5L, 7L, 6L, 3L, 7L, 2L, 2L, 3L, 1L, 2L, 1L, 5L,
5L, 2L, 4L, 1L, 1L, 2L, 1L, 4L, 7L, 3L, 2L, 5L, 3L, 2L, 4L, 2L,
1L, 7L, 5L, 2L, 2L, 2L, 3L, 4L, 1L, 2L, 5L, 2L, 3L, 3L, 1L, 3L,
2L, 3L, 5L, 1L, 3L, 1L, 5L, 4L, 5L, 4L, 5L, 5L, 5L, 1L, 3L, 3L,
1L, 3L, 6L, 3L, 4L, 3L, 3L, 5L, 3L), Censor = c(0L, 1L, 1L, 0L,
0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L,
1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L,
1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L,
1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L,
1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L,
0L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L,
1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L,
0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L,
0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L,
1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L,
1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L,
0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L,
1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L,
1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L)), class = "data.frame",
row.names = c(NA, -300L))
这是脚本:
library(mgcv)
library(ggplot2)
#Run the model
Model1 <- gam(SurvYear~
(Gender)+
s(Age, k=50)+
s(Highest_Educationmx, k=7),
weights=Censor, data=df, gamma=1.5, family=cox.ph())
summary(Model1)
#Build a perspective chart
vis.gam(Model1, view=c("Age","Highest_Educationmx"),
plot.type="persp", color="gray", se=-1, theta=45, phi=25,
xlab="Age", ylab= "Highest Education",
ticktype="detailed", zlim=c(-5.00, 2.00))
#Plot individual predictors using plot command from mgcv
plot(Model1, all.terms=T, rug=T, residuals=F, se=T, shade=T, seWithMean=T)
#Plot individual predictors using ggplot instead of plot command from mgcv
#UNSURE HOW DO TO THIS
我有偏见(我写的)但你可以为此使用 gratia 包。
您可以使用 draw()
函数替代 plot.gam()
,如果您想要完全控制,只需使用 evaluate_smooth()
生成平滑的整洁表示,然后使用 ggplot2.
轻松绘制
下面是根据上面 Gavin Simpson 的建议编写的脚本:
library(gratia)
#Plot individual predictors using ggplot instead of the plot command from mgcv
sm <- gratia::evaluate_smooth(Model1, "Age")
ggplot(sm, aes(x=Age, y=est)) + geom_line(size=1.0) +
geom_ribbon(aes(ymax=est+se, ymin=est-se), alpha=0.20) +
coord_cartesian(xlim=c(20.00, 75.00), ylim=c(-2.00, 1.00)) +
scale_x_continuous(breaks=seq(20.00, 75.00, 5.00)) +
scale_y_continuous(breaks=seq(-2.00, 1.00, 1.00)) +
labs(title="Age") +
xlab("Age") +
ylab("Linear Risk Score") +
theme(plot.title=element_text(size=10)) +
geom_hline(yintercept=0, linetype="dashed", size=0.5) +
geom_vline(xintercept=mean(df$Age), linetype="dashed", size=0.5)
我有一个非线性生存模型,我使用 mgcv
包对其进行了编码。我可以生成一个常规图,但我希望能够编写 ggplot2
代码。我该怎么做?
这是我的代码:
df <- structure(list(SurvYear =c(3L, 2L, 3L, 6L, 8L, 3L, 5L, 2L, 9L,
8L, 1L, 7L, 1L, 4L, 6L, 8L, 2L, 5L, 1L, 1L, 7L, 1L, 5L, 3L, 2L,
1L, 9L, 1L, 5L, 2L, 2L, 1L, 2L, 3L, 4L, 8L, 7L, 2L, 2L, 6L, 9L,
7L, 3L, 9L, 6L, 8L, 2L, 8L, 2L, 1L, 1L, 6L, 5L, 3L, 3L, 7L, 2L,
4L, 5L, 2L, 3L, 7L, 4L, 1L, 2L, 2L, 3L, 5L, 1L, 9L, 2L, 2L, 3L,
9L, 6L, 2L, 2L, 4L, 3L, 1L, 9L, 7L, 3L, 1L, 2L, 1L, 6L, 3L, 1L,
5L, 6L, 5L, 6L, 4L, 2L, 1L, 3L, 1L, 1L, 3L, 4L, 3L, 8L, 9L, 7L,
6L, 3L, 5L, 2L, 7L, 9L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 9L, 1L,
4L, 8L, 1L, 8L, 1L, 1L, 8L, 5L, 2L, 9L, 4L, 8L, 4L, 9L, 2L, 2L,
3L, 2L, 9L, 3L, 2L, 1L, 3L, 2L, 1L, 9L, 9L, 2L, 1L, 1L, 1L, 2L,
9L, 1L, 5L, 1L, 6L, 9L, 3L, 2L, 2L, 5L, 7L, 4L, 2L, 7L, 2L, 4L,
5L, 3L, 3L, 9L, 2L, 6L, 1L, 3L, 4L, 5L, 9L, 8L, 1L, 2L, 8L, 2L,
9L, 1L, 7L, 3L, 3L, 1L, 6L, 3L, 4L, 9L, 1L, 3L, 4L, 4L, 2L, 7L,
2L, 3L, 1L, 1L, 7L, 2L, 1L, 1L, 2L, 1L, 9L, 1L, 2L, 9L, 1L, 1L,
2L, 3L, 7L, 3L, 1L, 1L, 2L, 5L, 4L, 6L, 7L, 1L, 9L, 2L, 1L, 8L,
1L, 2L, 1L, 4L, 2L, 3L, 3L, 9L, 9L, 9L, 4L, 1L, 1L, 4L, 9L, 3L,
1L, 1L, 3L, 3L, 4L, 1L, 1L, 1L, 1L, 6L, 9L, 1L, 1L, 8L, 1L, 3L,
3L, 8L, 3L, 5L, 1L, 2L, 1L, 2L, 4L, 3L, 1L, 6L, 1L, 4L, 8L, 1L,
3L, 2L, 2L, 3L, 6L, 2L, 1L, 1L, 1L, 9L, 3L, 1L, 7L, 3L, 9L, 1L,
9L, 5L, 4L), Gender = c(1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L,
1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L,
1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L,
0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L,
1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L,
1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L,
0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L,
0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L,
0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L,
1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L,
1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L,
1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L,
1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L,
1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L,
0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L,
1L, 1L), Age = c(63L, 66L, 34L, 43L, 63L, 21L, 24L, 44L, 52L,
59L, 27L, 32L, 30L, 20L, 56L, 55L, 35L, 26L, 53L, 43L, 39L, 19L,
34L, 28L, 19L, 24L, 50L, 22L, 58L, 24L, 50L, 25L, 37L, 30L, 51L,
69L, 23L, 49L, 22L, 46L, 58L, 31L, 23L, 53L, 59L, 25L, 38L, 44L,
34L, 49L, 19L, 39L, 24L, 51L, 29L, 27L, 48L, 77L, 22L, 43L, 59L,
49L, 60L, 51L, 49L, 47L, 50L, 44L, 41L, 44L, 50L, 42L, 46L, 54L,
35L, 21L, 26L, 26L, 40L, 21L, 48L, 49L, 20L, 20L, 32L, 37L, 22L,
36L, 46L, 28L, 39L, 35L, 51L, 39L, 49L, 57L, 46L, 18L, 52L, 47L,
27L, 32L, 23L, 43L, 42L, 57L, 22L, 40L, 19L, 58L, 71L, 55L, 42L,
20L, 51L, 21L, 20L, 61L, 36L, 54L, 19L, 35L, 38L, 41L, 34L, 22L,
41L, 42L, 56L, 50L, 53L, 53L, 48L, 22L, 59L, 27L, 28L, 32L, 37L,
68L, 24L, 26L, 61L, 21L, 20L, 20L, 50L, 62L, 61L, 29L, 18L, 40L,
67L, 43L, 25L, 43L, 22L, 56L, 47L, 41L, 40L, 43L, 27L, 37L, 61L,
35L, 23L, 54L, 38L, 38L, 39L, 45L, 49L, 63L, 49L, 44L, 44L, 23L,
37L, 58L, 61L, 25L, 18L, 59L, 25L, 51L, 40L, 27L, 42L, 22L, 38L,
22L, 45L, 33L, 32L, 36L, 53L, 52L, 19L, 45L, 53L, 27L, 65L, 25L,
53L, 57L, 29L, 23L, 62L, 36L, 56L, 59L, 41L, 61L, 44L, 24L, 21L,
38L, 29L, 55L, 33L, 18L, 21L, 19L, 65L, 24L, 59L, 34L, 25L, 45L,
48L, 18L, 41L, 61L, 32L, 37L, 21L, 20L, 57L, 25L, 65L, 50L, 61L,
32L, 27L, 19L, 50L, 63L, 19L, 45L, 20L, 36L, 20L, 19L, 53L, 39L,
50L, 20L, 24L, 57L, 28L, 21L, 39L, 49L, 21L, 20L, 39L, 20L, 44L,
19L, 39L, 53L, 29L, 60L, 43L, 21L, 23L, 30L, 42L, 42L, 51L, 35L,
50L, 51L, 56L, 52L, 22L, 36L, 56L, 28L, 57L, 20L, 47L, 48L, 65L,
71L, 21L, 70L, 23L, 63L), Highest_Educationmx = c(4L, 5L, 3L,
2L, 3L, 2L, 3L, 1L, 3L, 1L, 7L, 3L, 2L, 3L, 3L, 2L, 6L, 2L, 3L,
6L, 3L, 2L, 2L, 7L, 2L, 1L, 2L, 3L, 6L, 3L, 5L, 3L, 5L, 6L, 2L,
1L, 5L, 2L, 5L, 1L, 1L, 3L, 2L, 3L, 1L, 7L, 5L, 4L, 7L, 3L, 1L,
1L, 6L, 3L, 3L, 2L, 4L, 6L, 5L, 4L, 2L, 6L, 1L, 3L, 4L, 2L, 1L,
5L, 5L, 3L, 1L, 5L, 3L, 3L, 1L, 4L, 2L, 3L, 5L, 3L, 1L, 4L, 2L,
1L, 2L, 7L, 2L, 5L, 3L, 2L, 6L, 1L, 1L, 3L, 4L, 1L, 5L, 1L, 3L,
4L, 2L, 7L, 2L, 4L, 4L, 7L, 4L, 6L, 3L, 1L, 2L, 1L, 5L, 5L, 1L,
5L, 2L, 7L, 3L, 4L, 2L, 4L, 2L, 4L, 2L, 2L, 4L, 1L, 2L, 1L, 2L,
6L, 1L, 2L, 5L, 2L, 2L, 5L, 1L, 6L, 5L, 2L, 1L, 2L, 1L, 1L, 3L,
2L, 4L, 3L, 2L, 3L, 1L, 5L, 5L, 7L, 1L, 3L, 3L, 2L, 1L, 3L, 4L,
5L, 1L, 1L, 3L, 3L, 3L, 5L, 3L, 6L, 4L, 3L, 1L, 3L, 5L, 7L, 1L,
3L, 4L, 5L, 3L, 3L, 1L, 1L, 1L, 7L, 3L, 1L, 4L, 3L, 3L, 5L, 1L,
4L, 5L, 4L, 2L, 5L, 3L, 1L, 1L, 5L, 4L, 7L, 5L, 2L, 2L, 5L, 3L,
1L, 1L, 2L, 3L, 5L, 3L, 7L, 5L, 1L, 5L, 3L, 1L, 1L, 1L, 1L, 7L,
5L, 7L, 3L, 1L, 5L, 7L, 6L, 3L, 7L, 2L, 2L, 3L, 1L, 2L, 1L, 5L,
5L, 2L, 4L, 1L, 1L, 2L, 1L, 4L, 7L, 3L, 2L, 5L, 3L, 2L, 4L, 2L,
1L, 7L, 5L, 2L, 2L, 2L, 3L, 4L, 1L, 2L, 5L, 2L, 3L, 3L, 1L, 3L,
2L, 3L, 5L, 1L, 3L, 1L, 5L, 4L, 5L, 4L, 5L, 5L, 5L, 1L, 3L, 3L,
1L, 3L, 6L, 3L, 4L, 3L, 3L, 5L, 3L), Censor = c(0L, 1L, 1L, 0L,
0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L,
1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 0L,
1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L,
1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L,
1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L,
0L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L,
1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L,
0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L,
0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L,
0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L,
1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L,
1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L,
0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L,
1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L,
1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L,
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L)), class = "data.frame",
row.names = c(NA, -300L))
这是脚本:
library(mgcv)
library(ggplot2)
#Run the model
Model1 <- gam(SurvYear~
(Gender)+
s(Age, k=50)+
s(Highest_Educationmx, k=7),
weights=Censor, data=df, gamma=1.5, family=cox.ph())
summary(Model1)
#Build a perspective chart
vis.gam(Model1, view=c("Age","Highest_Educationmx"),
plot.type="persp", color="gray", se=-1, theta=45, phi=25,
xlab="Age", ylab= "Highest Education",
ticktype="detailed", zlim=c(-5.00, 2.00))
#Plot individual predictors using plot command from mgcv
plot(Model1, all.terms=T, rug=T, residuals=F, se=T, shade=T, seWithMean=T)
#Plot individual predictors using ggplot instead of plot command from mgcv
#UNSURE HOW DO TO THIS
我有偏见(我写的)但你可以为此使用 gratia 包。
您可以使用 draw()
函数替代 plot.gam()
,如果您想要完全控制,只需使用 evaluate_smooth()
生成平滑的整洁表示,然后使用 ggplot2.
下面是根据上面 Gavin Simpson 的建议编写的脚本:
library(gratia)
#Plot individual predictors using ggplot instead of the plot command from mgcv
sm <- gratia::evaluate_smooth(Model1, "Age")
ggplot(sm, aes(x=Age, y=est)) + geom_line(size=1.0) +
geom_ribbon(aes(ymax=est+se, ymin=est-se), alpha=0.20) +
coord_cartesian(xlim=c(20.00, 75.00), ylim=c(-2.00, 1.00)) +
scale_x_continuous(breaks=seq(20.00, 75.00, 5.00)) +
scale_y_continuous(breaks=seq(-2.00, 1.00, 1.00)) +
labs(title="Age") +
xlab("Age") +
ylab("Linear Risk Score") +
theme(plot.title=element_text(size=10)) +
geom_hline(yintercept=0, linetype="dashed", size=0.5) +
geom_vline(xintercept=mean(df$Age), linetype="dashed", size=0.5)