*! version 1.0.3 PR 13oct2003. program define frac_s3b, rclass * calc derivatives of basis functions for restricted cubic splines. * Gives orthogonalised basis functions if option q() is specified. * Based on version 1.0.2 of frac_s3b, and frac_s3d. version 7 syntax varlist(max=1), K(string) [ BKnots(string) Name(string) Q(string) Second ] local x `varlist' if "`bknots'"!="" { local bk: word count `bknots' if `bk'!=2 { di in red "invalid bknots(), must specify 2 boundary knots" exit 198 } } if "`name'"!="" { if length("`name'")>6 { di in red "name() must be at most 6 characters" exit 198 } } else local name=substr("`x'",1,6) local nk: word count `k' quietly { /* Generate interior knots */ count if `x'!=. local nobs=r(N) sort `x' 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'] } /* Calc derivative of each basis function */ local names local j `nk' while `j' > 0 { local j1=`j'+1 local kj: word `j' of `k' * lambda notation of Royston & Parmar (2002) local lambdaj=(`kN'-`kj')/(`kN'-`k0') cap drop `name'`j' if "`second'"=="" { /* first derivative */ gen double `name'`j'=3*(`x'>`kj')*(`x'-`kj')^2 /* */ -3*`lambdaj'*(`x'>`k0')*(`x'-`k0')^2 /* */ -3*(1-`lambdaj')*(`x'>`kN')*(`x'-`kN')^2 } else { /* second derivative */ gen double `name'`j'=6*(`x'>`kj')*(`x'-`kj') /* */ -6*`lambdaj'*(`x'>`k0')*(`x'-`k0') /* */ -6*(1-`lambdaj')*(`x'>`kN')*(`x'-`kN') lab var `name'`j' "deriv of basis function for knot `kj'" } local names `name'`j' `names' local j=`j'-1 } if "`q'"!="" { /* Convert to orthogonalised basis via inverse of Q matrix !! second derivs are shifted by a constant cf. correct values - unimportant bug not corrected */ tempname Qinv matrix `Qinv'=inv(`q') * j indexes knots, basis fns (1, 2, ... ); * 0 is the linear term lnt - not orthogonalised forvalues j=1/`nk' { tempvar sum`j' local j1=`j'+1 /* j1 indexes cols in the Qinv matrix */ gen double `sum`j''=`Qinv'[1, `j1'] /* for the linear term in ln(_t) */ forvalues i=1/`j' { local i1=`i'+1 /* i1 indexes rows in the Qinv matrix */ replace `sum`j''=`sum`j''+`Qinv'[`i1', `j1']*`name'`i' } } forvalues j=1/`nk' { replace `name'`j'=`sum`j'' local kj: word `j' of `k' lab var `name'`j' "deriv of orthogonalised basis function for knot `kj'" } } } return local names `names' end exit Example: . mat Qinv=inv(Q) . mat list Qinv Qinv[4,4] r0 r1 r2 _cons z0 1.1917998 2.5496751 -1.9985308 0 z1 0 .51212157 -5.6359568 0 z2 0 0 7.1723579 0 _cons -8.0805192 -13.158682 8.4948663 1 Note that unused, zero'th orthog basis function r0=1.1917998*z0-8.0805192 where z0=ln(_t) /* First deriv of first and second orthog basis functions from derivs of untransformed basis functions */ . gen deriv1=2.5496751+I__c1*.51212157 . gen deriv2=-1.9985308+I__c1*(-5.6359568)+I__c2*7.1723579