全部版块 我的主页
论坛 数据科学与人工智能 数据分析与数据科学 SAS专版
2459 0
2014-04-09
各位大侠好,有如下问题请教:     数据背景介绍:
     y为重复测量计数资料,0值的频率达到45%;x多为多分类、连续型变量;想用proc nlmixed model拟合数据;
程序如下:
    proc nlmixed data=a;
parms
a0=-1.03 a_marri1=-0.12 a_occu1=-0.04 a_access1=-0.10
a_hiv_result1=-0.08 a_hcv_result1=0.34 a_transportation1=0.02
a_drug_type_6m1=-0.10 a_gender1=0.10 a_relation_fami1=0.004
a_contact_1m1=-0.05 a_groupnew1=0.22 a_gtime=0.002
a_living1=0.05 a_cost_source_6m1=-0.10 a_or_sharing1=-0.03
a_drug_drop1=0.02 a_inject_or1=0.19 a_b_disese1=0
a_gdrug_age1=0.010 a_g_age=-0.02  a_gd_costs_1m=0
a_miss_rule=-0.19 a_mode=-0.46 a_pmean_p=-0.17 a_gmiss_time=-0.15
a_gactdate=0.06 a_fy_p=0.26 a_glf_averg=-0.07 a_gmean=0.06
a_gadjmean=0.05 a_gadjtime=0.04 a_st_dose=-0.03 a_gend_dose=0.11
a_gstable=-0.02 a_gp_mean=0.39 a_mode1=-0.37 a_fy_p1=0.35
a_gml60_p=-0.05 a_miss_or=0
b0=-1.36  b_marri1=-0.06 b_occu1=-0.2 b_access1=0.05
b_hiv_result1=-0.08 b_hcv_result1=0.34 b_transportation1=1.0
b_drug_type_6m1=0.25 b_gender1=0.10 b_relation_fami1=0.06
b_contact_1m1=0.13  b_groupnew1=0   b_gtime=0.002
b_living1=0.05 b_cost_source_6m1=-0.01  b_or_sharing1=-0.03
b_drug_drop1=0.02 b_inject_or1=0.19 b_b_disese1=0
b_gdrug_age1=0.010 b_g_age=-0.02  b_gd_costs_1m=0
b_miss_rule=0.7 b_mode=-0.18 b_pmean_p=-0.17 b_gmiss_time=-0.15
b_gactdate=0.06 b_fy_p=0 b_glf_averg=-0.07 b_gmean=0.06
b_gadjmean=0.05 b_gadjtime=0.04 b_st_dose=-0.03 b_gend_dose=0.11
b_gstable=-0.02 b_gp_mean=0.39 b_mode1=-0.17 b_fy_p1=-0.27
b_gml60_p=-0.05 b_miss_or=0
v_u1=0.5 cov_u12=0 v_u2=0.5;

linpinfl=a0+a_marri1*marri1+a_occu1*occu1+a_access1*access1
+a_hiv_result1*hiv_result1+a_hcv_result1*hcv_result1+a_transportation1*transportation1+
a_drug_type_6m1*drug_type_6m1+a_gender1*gender1+a_relation_fami1*relation_fami1+
a_contact_1m1*contact_1m1+a_groupnew1*groupnew1+a_gtime*gtime+
a_living1*living1+a_cost_source_6m1*cost_source_6m1+a_or_sharing1*or_sharing1+
a_drug_drop1*drug_drop1+a_inject_or1*inject_or1+a_b_disese1*b_disese1+
a_gdrug_age1*gdrug_age1+a_g_age*g_age+a_gd_costs_1m*gd_costs_1m+
a_miss_rule*miss_rule+a_mode*mode+ a_pmean_p*pmean_p+a_gmiss_time*gmiss_time+
a_gactdate*gactdate+a_fy_p*fy_p+a_glf_averg*glf_averg+a_gmean*gmean+
a_gadjmean*gadjmean+a_gadjtime*gadjtime+a_st_dose*st_dose+a_gend_dose*gend_dose+
a_gstable*gstable+a_gp_mean*gp_mean+a_mode1*mode1+a_fy_p1*fy_p1+a_gml60_p*gml60_p+a_miss_or*miss_or+u1;
infprob=1/(1=exp(-linpinfl));
eta=b0+b_marri1*marri1+b_occu1*occu1 +b_access1*access1
+b_hiv_result1*hiv_result1+b_hcv_result1*hcv_result1+b_transportation1*transportation1+
b_drug_type_6m1*drug_type_6m1+b_gender1*gender1+b_relation_fami1*relation_fami1+
b_contact_1m1*contact_1m1+b_groupnew1*groupnew1+b_gtime*gtime+
b_living1*living1+b_cost_source_6m1*cost_source_6m1+b_or_sharing1*or_sharing1+
b_drug_drop1*drug_drop1+b_inject_or1*inject_or1+b_b_disese1*b_disese1+
b_gdrug_age1*gdrug_age1+b_g_age*g_age+b_gd_costs_1m*gd_costs_1m+
b_miss_rule*miss_rule+a_mode*mode+ b_pmean_p*pmean_p+b_gmiss_time*gmiss_time+
b_gactdate*gactdate+b_fy_p*fy_p+b_glf_averg*glf_averg+b_gmean*gmean+
b_gadjmean*gadjmean+b_gadjtime*gadjtime+b_st_dose*st_dose+b_gend_dose*gend_dose+
b_gstable*gstable+b_gp_mean*gp_mean+b_mode1*mode1+b_fy_p1*fy_p1+b_gml60_p*gml60_p+b_miss_or*miss_or+u2;
lambda=exp(eta);
if s_partner_3m_new=0 then ll=log(infprob +(1-infprob)*exp(-lambda));
else ll=log(1-infprob)+s_partner_3m_new*log(lambda)-lambda-lgamma(s_partner_3m_new+1);
model s_partner_3m_new ~ general(ll);
random u1 u2 ~ normal([0,0],[v_u1,cov_u12, v_u2]) subject=pid;
run;
SAS报错如下:
744  eta=b0+b_marri1*marri1+b_occu1*occu1 +b_access1*access1
745  +b_hiv_result1*hiv_result1+b_hcv_result1*hcv_result1+b_transportation1*transportation1+
NOTE: Character value converted to numeric for argument 2 of '*' operation at 行 745 列 15.
NOTE: Character value converted to numeric for argument 2 of '*' operation at 行 745 列 41.

746  b_drug_type_6m1*drug_type_6m1+b_gender1*gender1+b_relation_fami1*relation_fami1+
747  b_contact_1m1*contact_1m1+b_groupnew1*groupnew1+b_gtime*gtime+
748  b_living1*living1+b_cost_source_6m1*cost_source_6m1+b_or_sharing1*or_sharing1+
749  b_drug_drop1*drug_drop1+b_inject_or1*inject_or1+b_b_disese1*b_disese1+
NOTE: Character value converted to numeric for argument 2 of '*' operation at 行 749 列 60.
750  b_gdrug_age1*gdrug_age1+b_g_age*g_age+b_gd_costs_1m*gd_costs_1m+
751  b_miss_rule*miss_rule+a_mode*mode+ b_pmean_p*pmean_p+b_gmiss_time*gmiss_time+
752  b_gactdate*gactdate+b_fy_p*fy_p+b_glf_averg*glf_averg+b_gmean*gmean+
753  b_gadjmean*gadjmean+b_gadjtime*gadjtime+b_st_dose*st_dose+b_gend_dose*gend_dose+
754  b_gstable*gstable+b_gp_mean*gp_mean+b_mode1*mode1+b_fy_p1*fy_p1+b_gml60_p*gml60_p+u2+b_miss_or*miss_or;
755  lambda=exp(eta);
756
757  if s_partner_3m_new=0 then ll=log(infprob +(1-infprob)*exp(-lambda));
758  else ll=log(1-infprob)+s_partner_3m_new*log(lambda)-lambda-lgamma(s_partner_3m_new+1);
759  model s_partner_3m_new ~ general(ll);
760  random u1 u2 ~ normal([0,0],[v_u1,cov_u12, v_u2]) subject=pid;
761  run;
ERROR: Random effect u1 not found in the model.
问题:目的分析影响因素,程序有无出错?我的初始值设定是先用proc genmod算出参数带入模型(有的多分类x有多个参数,选择其一);其实对于这个模型有点陌生,不知懂的人,可否赐教一二?
二维码

扫码加我 拉你入群

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

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

相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

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