[作图问题求助] 询问stata无法作cubic spline量效曲线图的问题 [推广有奖]

0关注
0粉丝

小学生

35%

还不是VIP/贵宾

威望: 0 级
论坛币: 0 个
通用积分: 0
学术水平: 0 点
热心指数: 0 点
信用等级: 0 点
经验: 93 点
帖子: 4
精华: 0
在线时间: 2 小时
注册时间: 2013-4-7
最后登录: 2013-8-5

楼主

吴山逸话 发表于 2013-4-7 16:50:19 |AI写论文

是否 +2 论坛币

k人参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我在做一个cubic spline的统计，统计饮酒量与肺癌发病RR的量效关系。
方法是是采用glst法，
想重复一下该方法原创者统计学家尼可拉斯提供的一个例子，http://nicolaorsini.altervista.org/stata/tutorial/g/alcohol_lc.txt结果发现最后画图的命令（红字部分）无法实现，说是无效命令或变量缺失。请有经验的朋友帮我看看，给予指点，不知是我的stata问题还是其他问题？谢

// Example 2

use http://nicolaorsini.altervista.org/data/ex_alcohol_lc, clear

// Two-stage fixed-effect dose-response model assuming linearity

glst logrr dose , se(se) cov(peryears cases) pfirst(study type) ts(f)

/* Relative risk for 12 grams/day incremental unit */
lincom dose*12 , eform

// Two-stage random-effect dose-response model assuming linearity

glst logrr dose , se(se) cov(peryears cases) pfirst(study type) ts(r)
/* Relative risk for 12 grams/day incremental unit */
lincom dose*12 , eform

// Non-linearity using fixed-effect

capture drop doses*
_pctile dose , percentile(5 35 65 95)
ret list
mkspline doses  = dose  , knots(`=r(r1)' `=r(r2)' `=r(r3)' `=r(r4)') cubic displayknots
glst logrr doses* , se(se)  cov(peryears cases) pfirst(study type)
testparm doses2 doses3

// Figure 1 B of the paper

glst  logrr dose , se(se) cov(peryears cases) pfirst(study type)  ts(f)
predictnl lrr_lin = _b[dose]*dose
gen rr_lin = exp(lrr_lin)

// using 0 as referent

glst logrr doses*, se(se) cov(peryears cases) pfirst(study type)
predictnl logrrwithref = _b[doses1]*doses1 + _b[doses2]*doses2 + _b[doses3]*doses3, ci(lo hi)
gen rrwithref  = exp(logrrwithref)
gen lbwithref  = exp(lo)
gen ubwithref  = exp(hi)

// Tabulate result using xblc (findit xblc)

levelsof dose, local(level)
xblc doses*, c(dose) at(`r(levels)') ref(0) eform

twoway ///
      (line lbwithref ubwithref rrwithref  dose,  sort lp(longdash longdash l ) lc(black black black) ) ///
(line rr_lin  dose,  sort lp(shortdash) lc(black) )  , ///
scheme(s1mono) ylabel(.9 1 1.2 1.5 1.8, angle(horiz) format(%3.2fc)) ///
xlabel(0(5)45) ///
legend(off) ///
ytitle("Relative Risk", margin(right)) ///
xtitle("Alcohol intake, grams/day" , margin(top_bottom) ) name(figure1B, replace) yscale(log)  ///
  plotregion(style(none))