*! version 1.0.5 PR/GA 18feb2005. program define frac_spl, rclass /* generate spline basis functions */ version 6 gettoken v 0 : 0, parse(" ,") unab v : `v', max(1) local done 0 gettoken nxt : 0, parse(" ,") while ("`nxt'" != "" & "`nxt'" != ",") & !`done' { cap confirm number `nxt' if !_rc { local knots "`knots' `nxt'" gettoken nxt 0 : 0, parse(" ,") gettoken nxt : 0, parse(" ,") } else local done 1 } if "`knots'"=="" { local options "DF(int 0) UNIQUE" } syntax [if] [in] [, `options' all Orthog Q(string) Name(string) Basis(string) /* */ DEGree(int 3) BKnots(string) ] if "`bknots'"!="" { local bk: word count `bknots' if `bk'!=2 { di in red "invalid bknots(), must specify 2 boundary knots" exit 198 } } if "`basis'"!="" { local name `basis' } if "`name'"!="" { if length("`name'")>6 { di in red "name() must be at most 6 characters" exit 198 } } else local name=substr("`v'",1,6) if `degree'>3 | `degree'<0 | `degree'==2 { di in red "degree must be either 0, 1, or 3." exit 198 } if `degree'==0 { local lt "<" /* allows knots at min for deg 0 */ } else local lt "<=" * All is same as in fracgen, to generate basis functions for all obs but knots in generatation sample if "`all'" == "" { local all `"`if' `in'"' } else local all marksample touse markout `touse' `v' tempvar x quietly { gen double `x' = `v' `all' gsort -`touse' `x' if "`unique'"!="" { tempvar xold gen double `xold'=`x' bysort `x': replace `x'=. if _n!=1 sort `x' count if `touse' & !missing(`x') local nobs=r(N) } else { count if `touse' local nobs=r(N) } /* Generate interior knots */ if "`bknots'"!="" { local k0: word 1 of `bknots' confirm num `k0' local kN: word 2 of `bknots' confirm num `kN' } else { local k0=`x'[1] local kN=`x'[`nobs'] } local klast `kN' if "`knots'"=="" { if `df'==0 { local nk = (`degree'==0)+int((`nobs')^.25) } else local nk=`df'-(`degree'!=0) if `nk'!=1 { local s s } noi di in blue "[`nk' knot`s' to be used]" local j `nk' while `j' > 0 { local k`j'=`x'[int((`j'*`nobs'/(`nk'+1))+.5)] if `k`j''==`klast'|`k`j''`lt'`k0' { noisily di in bl "[knot at `k`j'' ignored]" local i `j' while `i'<`nk' { local i2=`i'+1 local k`i' `k`i2'' local i=`i'+1 } local nk=`nk'-1 } else { local knots "`k`j'' `knots'" local klast `k`j'' } local j = `j' -1 } if "`unique'"!="" { replace `x'=`xold' } } else { tokenize `knots' local nk 0 while "`1'"!="" { if `1'`lt'`k0' | `1'>=`kN' { noi di in bl "[knot at `1' ignored]" } else { local nk = `nk' + 1 local k`nk' `1' } mac shift } if `nk'!=1 { local s s } noi di in blue "[`nk' knot`s' to be used]" } /* Create basis variables */ local j `nk' while `j' > 0 { local vn `name' if `j'<10 { local vn "`vn'_" } cap drop `vn'`j' if "`orthog'"=="" { local n `vn'`j' } else { tempvar cov`j' local n `cov`j'' } if `degree'==0 { * care needed for missing x gen byte `n'=(`x'>`k`j'') if !missing(`x') frac_mdp 3 `k`j'' local low: word 1 of `r(mdp)' lab var `n' "`v'>`low'" } else if `degree'==1 { * is correct for missing x gen double `n'=(`x'>`k`j'')*(`x'-`k`j'') } else if `degree'==3 { * now using lambda notation of Royston & Parmar (2001) local lambdaj=(`kN'-`k`j'')/(`kN'-`k0') gen double `n'=(`x'>`k`j'')*(`x'-`k`j'')^3 /* */ -`lambdaj'*(`x'>`k0')*(`x'-`k0')^3 /* */ -(1-`lambdaj')*(`x'>`kN')*(`x'-`kN')^3 } local covars "`n' `covars'" local j = `j' -1 } if "`orthog'"=="" & `degree'!= 0 { if "`v'"!="`name'" { cap gen double `name'_0=`x' if _rc { drop `name'_0 gen double `name'_0=`x' } local covars "`name'_0 `covars'" } else { local covars "`name' `covars'" } } else if "`orthog'"!="" { local nk1=`nk'+(`degree'!=0) local vn `name' if `nk1'<10 { local vn "`vn'_" } if `degree'!=0 { cap drop `vn'`nk1' local xvar `x' } /* If the matrix Q is fed in thru q() option, we use its inverse to compute the transformed (orthogonalised) basis functions. Otherwise, compute them from scratch using orthog.ado. */ tempname Q if "`q'"!="" { /* Convert to orthogonalised basis fns via inverse of Q matrix */ tempname Qinv matrix `Q'=`q' matrix `Qinv'=inv(`Q') * j indexes knots, basis fns (1, 2, ... ); * 0 is the linear term lnt - not orthogonalised local nk2=`nk1'+1 forvalues j=1/`nk' { local j1=`j'+1 /* j1 indexes cols in the Qinv matrix */ /* constant + term in ln(_t) */ gen double `vn'`j1'=`Qinv'[`nk2', `j1']+`Qinv'[1, `j1']*`xvar' forvalues i=1/`j' { local n: word `i' of `covars' /* covars don't include ln(_t) */ local i1=`i'+1 /* i1 indexes rows in the Qinv matrix */ replace `vn'`j1'=`vn'`j1'+`Qinv'[`i1', `j1']*`n' } } } else { orthog `xvar' `covars', gen(`vn'*) matrix(`Q') drop `vn'1 /* discard first orthog var, using original x */ } if "`v'"!="`name'" { cap gen double `name'_0=`x' if _rc { drop `name'_0 gen double `name'_0=`x' } local covars "`name'_0" } else { local covars "`v'" } forvalues j=2/`nk1' { local jm=`j'-1 local vnew `vn'`jm' local vold `vn'`j' rename `vold' `vnew' local covars "`covars' `vnew'" lab var `vnew' "orthogonalized basis `jm'" cap drop `cov`jm'' /* label basis vectors with chars corresponding to knot positions */ if `jm'>1 { cap char `vnew'[fp] `k`jm'' } else cap char `name'[fp] "linear" /* !! this may be incorrect */ } } } global S_1 `covars' global S_2 `nk' global S_3 `knots' return local names `covars' return scalar nk=`nk' return local knots `knots' return local bknots `k0' `kN' if "`orthog'"!="" { return matrix Q `Q' } end program define frac_mdp, rclass * 1=decimal places, rest=numbers sep by spaces * Output in `r(mdp)'. args dp mac shift while "`1'"!="" { cap confirm num `1' if _rc { local ddp `ddp' `1' } else { if int(2*`1')==(2*`1') { local ddp `ddp' `1' /* respect .5 */ } else { frac_ddp `1' `dp' local ddp `ddp' `r(ddp)' } } mac shift } return local mdp `ddp' end *! version 1.0.1 PR 25Feb1999. program define frac_ddp, rclass version 6 * 1=input number, 2=decimal places required. * Output in (string) r(ddp). * `1' and `2' allowed to be scalars. local n=`1' local dp=int(`2') if `dp'<0 | `dp'>20 { ret local ddp `n' exit } local z=int(abs(`n')*(10^`dp')+.5) if `z'>1e15 { /* can't cope with number this large -- E notation messes it up */ ret local ddp `n' exit } local lz=length("`z'") if `lz'<`dp' { local z=substr("00000000000000000000",1,`dp'-`lz')+"`z'" } if `dp'>0 { local f=length("`z'")-`dp' ret local ddp=substr("`z'",1,`f')+"."+substr("`z'",`f'+1,`dp') /* add leading zero */ if abs(`n')<1 { ret local ddp 0`return(ddp)' } } else ret local ddp `z' if `n'<0 { ret local ddp -`return(ddp)' } * failsafe check cap confirm num `return(ddp)' if _rc { ret local ddp `n' } end