data boc; input days boc; cards; 1 105 1 97 1 104 1 106 2 136 2 161 2 151 2 153 3 173 3 179 3 174 3 174 5 195 5 182 5 201 5 172 7 207 7 194 7 206 7 213 10 218 10 193 10 235 10 229 ; run; /* Make a scatter plot of boc vs. days and comment on the adequacy for making a linear regression on this scale. */ proc gplot data=boc; plot boc*days; symbol v=dot i=rl; run; quit; /* boc = GAMMA exp(-BETA/days) log(boc) = log(GAMMA) - (1/days) * BETA logboc = LGAMMA - invdays * BETA */ data new; set boc; logboc=log(boc); invdays=1/days; run; /* picture using simulated data - what is the interpretation of the parameters */ data simulate; do days=1 to 70 by 0.1; gamma1=200; beta1=1; boc1=gamma1*exp(-beta1/days); gamma2=250; beta2=1; boc2=gamma2*exp(-beta2/days); gamma3=200; beta3=2; boc3=gamma3*exp(-beta3/days); gamma4=250; beta4=2; boc4=gamma4*exp(-beta4/days); output; end; run; proc gplot data=simulate; plot (boc1 boc2 boc3 boc4)*days/overlay; symbol v=dot i=none; run; quit; /* linear relation ? */ proc gplot data=new; plot logboc*invdays; symbol v=dot i=rl; run; quit; /* linear regression */ proc glm data=new; model logboc=invdays/solution clparm; run; quit; /* test linearity */ proc glm data=new; class days; model logboc=invdays days/solution clparm; run; quit; /* parameters */ data gamma; int=5.431253025; lower_int=5.391973906; upper_int=5.470532143; gamma=exp(int); lower_gamma=exp(lower_int); upper_gamma=exp(upper_int); keep gamma lower_gamma upper_gamma; run; proc print data=gamma noobs; run; /* two parameters - linear regression */ proc glm data=new; model logboc=invdays/solution clparm; output out=pred_lin p=p; run; quit; proc gplot data=pred_lin; title 'predicted values - two parameters'; plot p*days; symbol v=none i=join; run; quit; /* six parameters - anova */ proc glm data=new; class days; model logboc=invdays days/solution clparm; output out=pred_nonlin p=p; run; quit; proc gplot data=pred_nonlin; title 'predicted values - six parameters'; plot p*days; symbol v=none i=join; run; quit; /*********************************************************/ /* nu begynder programmet til opgave 2 om height vs. age */ /*********************************************************/ options nocenter ps = 60 ls = 75; goptions device=Win hby=2; data juul; infile 'T:\juul2.txt'; input age height menarche sexnr sigf1 tanner testvol weight; if sexnr=2 then sex='female'; if sexnr=1 then sex='male'; lheight=log10(height); y=log(log(200/height)); oplyst=0; if sigf1 ne . and age ne . and sexnr ne . and tanner ne . then oplyst=1; extra_age10=max(age-10,0); extra_age12=max(age-12,0); extra_age13=max(age-13,0); extra_age15=max(age-15,0); if age<10 then age_group=1; if 10<=age<12 then age_group=2; if 12<=age<13 then age_group=3; if 13<=age<15 then age_group=4; if 15<=age then age_group=5; age5=age-5; age14=age-14; age2=age**2; age3=age**3; run; proc reg data=juul; where age ge 5 and age le 20 and sex='male'; model height=age; run; proc gplot data=juul; where age ge 5 and age le 20 and sex='male'; plot height*age / haxis=axis1 vaxis=axis2 frame; axis1 value=(H=3) minor=NONE label=(H=3); axis2 value=(H=3) minor=NONE label=(A=90 R=0 H=3); symbol1 v=circle i=rl c=black h=2 l=1 w=2; run; proc glm data=juul; where age ge 5 and age le 20; class sex; model height=sex sex*age sex*age2 sex*age3 / noint solution; run; proc glm data=juul; where age ge 5 and age le 20; class sex; model height=age age2 age3 sex sex*age sex*age2 sex*age3 / solution; run; proc sort data=juul; by sex; run; proc glm data=juul; where age ge 5 and age le 20; by sex; model height=age5 extra_age10 extra_age12 extra_age13 extra_age15 / solution; contrast 'all' extra_age10 1, extra_age12 1, extra_age13 1, extra_age15 1; run; proc glm data=juul; where age ge 5 and age le 20; class sex; model height=age5 extra_age10 extra_age12 extra_age13 extra_age15 sex sex*extra_age15 sex*extra_age13 sex*extra_age12 sex*extra_age10 sex*age5 / solution; output out=fit p=predicted r=residual; run; proc nlin data=juul; where age ge 5 and age le 20; by sex; parms alfa=200 beta=1.5 gamma=0.15; model height=alfa*exp(-beta*exp(-gamma*age)); output out=ny p=yhat r=resid; run; proc nlin data=juul; where age ge 5 and age le 20; by sex; parms alfa=200 beta=1.5 gamma=0.15; model lheight=log(alfa)-beta*exp(-gamma*age); output out=ny p=yhat r=resid; run; proc nlin data=juul; where age ge 5 and age le 20; parms alfa=200 beta=1.5 gamma=0.15 alfadif=0 betadif=0 gammadif=0; if sex='male' then model lheight=log(alfa)-beta*exp(-gamma*age); if sex='female' then model lheight=log(alfa+alfadif)- (beta+betadif)*exp(-(gamma+gammadif)*age); run; data fit; set ny; art=1; output; art=2; height=exp(yhat); output; run; proc gplot data=fit uniform; by sex; plot height*age=art / nolegend haxis=axis1 vaxis=axis2 frame; axis1 offset=(3,3) value=(H=3) minor=NONE label=(H=3); axis2 value=(H=3) minor=NONE label=(A=90 R=0 H=3); symbol1 v=circle i=none c=BLACK l=1 w=2 h=3; symbol2 v=none i=spline c=BLACK l=1 w=2 h=3; run; proc reg data=juul; where age ge 5 and age le 20 and sex='male'; model height=age5 extra_age10 extra_age12 extra_age13 extra_age15; output out=fitreg p=predicted r=residual; run; proc gplot data=fitreg; where age ge 5 and age le 20 and sex='male'; plot height*age predicted*age / overlay href=10 lh=33 href=12 lh=33 href=13 lh=33 href=15 lh=33 haxis=axis1 vaxis=axis2 frame; axis1 value=(H=3) minor=NONE label=(H=3); axis2 value=(H=3) minor=NONE label=(A=90 R=0 H=3 'linear spline for height'); symbol1 v=circle i=none c=black h=2 l=1 w=2; symbol2 v=none i=join c=black h=2 l=1 w=2; run; proc glm data=juul; where age ge 5 and age le 20; class sex; model height=age5 extra_age10 extra_age12 extra_age13 extra_age15 sex sex*extra_age15 sex*extra_age13 / solution; * output out=fit p=predicted r=residual; run; data fit; set fit; proc sort data=fit; by sex age; run; proc gplot data=fit; plot predicted*age=sex / haxis=axis1 vaxis=axis2 frame; axis1 value=(H=3) minor=NONE label=(H=3); axis2 value=(H=3) minor=NONE label=(A=90 R=0 H=3); symbol1 v=none i=join c=black h=2 l=1 w=2; symbol2 v=none i=join c=black h=2 l=33 w=2; run;