a_matrix(scan("diabetes.dat"),43,3,byrow=T) a_a[sort.list(a[,1]),] age_a[,1] basedef_a[,2] cpep_a[,3] lncpep_log(cpep) #Smoothing spline estimates with lambda calculated using REML via PROC MIXED b_matrix(scan("spline_mixed.out"), 37,3,byrow=T) uage_b[,1] ssreml_b[,2] bsereml_b[,3] #Bayesian CIs bupreml_ssreml+1.96*bsereml blowreml_ssreml-1.96*bsereml #Smoothing spline estimates with lambda calculated using REML via SPMM b1_matrix(scan("spmm_fit.out"), 37,8,byrow=T) uage1_b1[,1] ssreml1_b1[,2] fsereml1_b1[,3] #Freq CIs: flowreml1_b1[,4] fupreml1_b1[,5] #Bayesian CIs bsereml1_b1[,6] blowreml1_b1[,7] bupreml1_b1[,8] #Compare Results of PROC MIXED and SPMM cbind(ssreml,ssreml1,bsereml,bsereml1,fsereml1,bupreml,bupreml1) ssreml ssreml1 bsereml bsereml1 fsereml1 bupreml bupreml1 [1,] 1.252865 1.252954 0.06405417 0.06404336 0.05894708 1.378411 1.378477 [2,] 1.260305 1.260390 0.06201116 0.06200222 0.05759590 1.381847 1.381912 [3,] 1.319539 1.319588 0.04965587 0.04965358 0.04776965 1.416864 1.416907 [4,] 1.326882 1.326927 0.04859118 0.04858878 0.04669134 1.422121 1.422159 [5,] 1.348763 1.348795 0.04596269 0.04595900 0.04369422 1.438850 1.438873 [6,] 1.480340 1.480302 0.03879836 0.03878621 0.03370211 1.556385 1.556321 [7,] 1.511672 1.511623 0.03766821 0.03765737 0.03306822 1.585502 1.585430 [8,] 1.525580 1.525526 0.03729373 0.03728324 0.03282556 1.598676 1.598600 [9,] 1.529936 1.529881 0.03720238 0.03719191 0.03274263 1.602852 1.602775 [10,] 1.552889 1.552830 0.03689625 0.03688505 0.03217734 1.625206 1.625124 [11,] 1.580810 1.580760 0.03613847 0.03612584 0.03090241 1.651641 1.651566 [12,] 1.593023 1.592993 0.03428145 0.03427012 0.02955948 1.660215 1.660161 [13,] 1.594336 1.594310 0.03382641 0.03381547 0.02925055 1.660635 1.660587 [14,] 1.596089 1.596071 0.03289875 0.03288846 0.02858455 1.660571 1.660531 [15,] 1.596793 1.596780 0.03219176 0.03218193 0.02805686 1.659889 1.659856 [16,] 1.596942 1.596932 0.03195436 0.03194469 0.02787879 1.659573 1.659542 [17,] 1.596785 1.596786 0.02972735 0.02971961 0.02639566 1.655050 1.655035 [18,] 1.596619 1.596621 0.02948081 0.02947336 0.02625658 1.654401 1.654388 [19,] 1.595390 1.595396 0.02838247 0.02837640 0.02566978 1.651019 1.651012 [20,] 1.595064 1.595071 0.02820689 0.02820104 0.02557674 1.650349 1.650344 [21,] 1.594712 1.594720 0.02805275 0.02804709 0.02549272 1.649695 1.649691 [22,] 1.593938 1.593947 0.02781663 0.02781123 0.02535471 1.648459 1.648456 [23,] 1.593089 1.593099 0.02768444 0.02767916 0.02526367 1.647350 1.647349 [24,] 1.592644 1.592655 0.02765974 0.02765445 0.02523877 1.646857 1.646856 [25,] 1.591718 1.591731 0.02769501 0.02768962 0.02523748 1.646000 1.646001 [26,] 1.590748 1.590761 0.02784006 0.02783448 0.02531255 1.645314 1.645316 [27,] 1.589240 1.589256 0.02825539 0.02824948 0.02560113 1.644621 1.644624 [28,] 1.588740 1.588756 0.02844613 0.02844011 0.02575179 1.644495 1.644498 [29,] 1.586421 1.586439 0.02980719 0.02980091 0.02699590 1.644843 1.644848 [30,] 1.586020 1.586038 0.03016894 0.03016267 0.02735420 1.645151 1.645156 [31,] 1.585313 1.585332 0.03099690 0.03099071 0.02819127 1.646067 1.646072 [32,] 1.585007 1.585025 0.03146676 0.03146062 0.02867255 1.646682 1.646687 [33,] 1.584005 1.584022 0.03375530 0.03374953 0.03104209 1.650165 1.650169 [34,] 1.583173 1.583187 0.03859385 0.03858873 0.03600195 1.658817 1.658819 [35,] 1.583477 1.583482 0.04684628 0.04684058 0.04387346 1.675296 1.675288 [36,] 1.584509 1.584501 0.06493905 0.06492353 0.05815425 1.711789 1.711749 [37,] 1.584617 1.584608 0.06721772 0.06720033 0.05971609 1.716364 1.716318 #plot the smoothing spline fit and its Bayesian CIs #using the PROX MIXED output par(mfrow=c(1,1)) plot(age,lncpep,xlim=c(0,16),ylim=c(1.1,1.9)) par(new=T) matplot(cbind(uage,uage,uage,uage), cbind(ssreml,bupreml,blowreml),type="l",lty=c(1,2,2), xlab=" ",ylab=" ",xlim=c(0,16),ylim=c(1.1,1.9)) title("Smoothing Spline Estimate (Lambda=REML) and 95% Frequent CI") #plot the smoothing spline fit and its Frequenst and Bayesian CIs #using the SPMM output par(mfrow=c(1,1)) plot(age,lncpep,xlim=c(0,16),ylim=c(1.1,1.9)) par(new=T) matplot(cbind(uage,uage,uage,uage,uage,uage,uage,uage,uage), cbind(ssreml1,fupreml1,flowreml1,bupreml1,blowreml1),type="l",lty=c(1,2,2,3,3), xlab=" ",ylab=" ",xlim=c(0,16),ylim=c(1.1,1.9)) legend(5,1.3,lty=c(1:3),legend=c("Estimated Curve", "Frequentist CI", "Bayesian CI")) title("SS Estimate (Lambda=REML) and its 95% Frequest and Bayesian CIs") #Smoothing spline estimates with lambda calculated using GCV ssgcv_smooth.spline(age,lncpep,cv=F) par(mfrow=c(1,1)) plot(age,lncpep,xlim=c(0,16),ylim=c(1.1,1.9)) par(new=T) plot(ssgcv$x,ssgcv$y,type="l",lty=1,xlim=c(0,16),ylim=c(1.1,1.9),xlab=" ", ylab=" ") par(new=T) plot(uage,ssreml,type="l",lty=2,xlim=c(0,16),ylim=c(1.1,1.9),xlab=" ", ylab=" ") legend(5,1.3,lty=c(1:2),legend=c("lambda=GCV", "lambda=REML")) title("Smoothing Spline Estimates with Lambda Estimated Using GCV and REML") #Cannot calculate the Frequenst 95% CI of the curve using S-plus #Can calculate the Bayesian 95% CI of the curve using the Splus smooth.spline output res_ssgcv$yin-ssgcv$y sigma2_sum(res^2)/sum(1-ssgcv$lev) #Estimate of residual variance sigma_sqrt(sigma2) countage_table(category(age,unique(age))) bupgcv_ssgcv$y+1.96*sigma*sqrt(ssgcv$lev)/sqrt(countage) #Bayesian upper 95% CI blowgcv_ssgcv$y-1.96*sigma*sqrt(ssgcv$lev)/sqrt(countage) #Bayesian lower 95% CI #Compare the Bayesian CIs calculated from PROC MIXED and from S par(mfrow=c(1,1)) plot(age,lncpep,xlim=c(0,16),ylim=c(1.1,1.9)) par(new=T) plot(ssgcv$x,ssgcv$y,type="l",lty=1,xlim=c(0,16),ylim=c(1.1,1.9),xlab=" ", ylab=" ") matplot(cbind(uage,uage,uage,uage,uage), cbind(ssgcv$y,bupgcv,blowgcv,bupreml,blowreml),type="l",lty=c(1,2,2,3,3), xlab=" ",ylab=" ",xlim=c(0,16),ylim=c(1.1,1.9)) legend(5,1.3,lty=c(1:3),legend=c("Smoothing Spline Estimate", "GCV Bayesian 95% CI(S)", "REML Bayesian 95% CI(SAS)")) title("Smoothing Spline Estimate with its Bayesian 95% CIs")