*! 1.5 16 Apr 2017 Austin Nichols add new default for age indicator option
*  1.4 31 Mar 2017 Austin Nichols removed lnmvnormal and replaced with (-.5*(x-m)*invsym(v)*(x-m)')-ln(pi()*det(v))
*  1.3 31 Mar 2017 Austin Nichols added round() and missok options
*  1.2 20 Mar 2017 Austin Nichols removed grouplist and gvar options, added agevar and modevar options
*  1.1 19 Mar 2017 Austin Nichols renumbered groups, added check for correct average correlations across items, added round() and missok options, changed variable name defaults
*  1.0 17 Oct 2016 Austin Nichols
prog pfs, rclass
 version 10.2
 syntax name [if] [in] [, modevar(varname numeric) agevar(varname numeric) qlist(varlist numeric min=10 max=10) replace INTPoints(int 31) zlimit(real 4) ghq round(real 1) missok skipcorr]
 marksample touse
 if "`replace'"=="" qui g double `namelist'=.
 else {
  cap g double `namelist'=.
  }
 *check rounding value
 if `round'<0 {
  di as err "Rounding value for scores must be positive"
  error 198
  }
 *check mode var
 if "`modevar'"=="" {
  loc modevar "self"
  }
 cap assert inlist(`modevar',0,1) if `touse'
 if _rc {
  di as err "mode variable `modevar' not defined correctly"
  err 459
  }
 *check age var
 if "`agevar'"=="" {
  loc agevar "age18_61"
  }
 cap assert inlist(`agevar',0,1) if `touse'
 if _rc {
  di as err "age variable `agevar' not defined correctly"
  err 459
  }
 *check group var
 tempvar gvar
 qui g byte `gvar'=1+2*`agevar'+`modevar' if `touse'
 cap conf numeric var `gvar'
 if _rc!=0 {
  di as err "Error: group var `gvar' not found"
  error 198
  }
 cap assert inlist(`gvar',1,2,3,4) if `touse'
 if _rc!=0 {
  di as err "Error: group var not in {1,2,3,4}"
  error 459
  }
 *check ghq option
 loc GHQ=cond("`ghq'"!="",1,0)
 if "`ghq'"!="" {
  if "`zlimit'"!="6" {
   di as err "zlimit option is ignored when option ghq is specified"
   }
  }
 *check questions
if "`qlist'"=="" loc qlist "fs1_complexdecision fs2_goodnewdecision fs3_followthrough fs4_recognizegoodinvestment fs5_keepfromspending fs6_howtosave fs7_findadvice fs8_notenoughinfo fs9_whenadvice fs10_struggleunderstand"
if wordcount("`qlist'")!=10 {
 di as err "Must specify 10 variables holding all 10 items in FS instrument (even if some are everywhere missing)"
 error 198
 }
if "`skipcorr'"=="" {
 tokenize `qlist'
 loc avgcorr=0
 loc ctcorr=0
 loc c 0
 foreach pos in 1 2 3 4 5 6 7 8 9 {
  foreach neg in 10 {
   cap corr ``pos'' ``neg''
   if _rc==0 loc c=r(rho)
   if `c'<. {
    loc avgcorr=`avgcorr'+`c'
    loc ctcorr=`ctcorr'+1
    }
   }
  }
 if `avgcorr'<0 {
  di as err "Reverse-coded items seem to have a negative correlation with other items; check coding of items."
  di as err "If you are sure of your coding of items, use the" as txt " skipcorr " as err "option." _n
  err 459
  }
 }
if "`grouplist'"=="" {
  loc grouplist "1 2 3 4"
  }
tempvar allmiss
g byte `allmiss'=1
 foreach v in `qlist' {
  cap conf numeric var `v'
  if _rc!=0 {
   di as err "Error: question var `v' not found"
   di as err "Every question needs to be present, though all values may be missing."
   error 198
   }
  cap assert inlist(`v',0,1,2,3,4)|mi(`v')
  if _rc!=0 {
   di as err "Error: question var `v' responses not in {0,1,2,3,4}"
   di as err "Every question needs to be coded 0 to 4, or missing."
   error 198
   }
  qui replace `allmiss'=0 if !mi(`v')
  }
 foreach v in `qlist' {
  cap conf numeric var `v'
  if _rc!=0 {
   di as err "Error: question var `v' not found"
   di as err "Every question needs to be present, though all values may be missing."
   error 198
   }
  cap assert inlist(`v',0,1,2,3,4)|mi(`v')
  if _rc!=0 {
   di as err "Error: question var `v' responses not in {0,1,2,3,4}"
   di as err "Every question needs to be coded 0 to 4, or missing."
   error 198
   }
  }
di as res "Computing scores... this may take some time."
 * define parameter vectors
 mat g1m=(.6441367,.0752469,.5430175,.3178859,.0341296)
 mat g1v=(1.2708368,1.0376652,2.6169667,2.7600912,2.3128359)
 mat g1q1a=(2.4110204,1.6608551,0,0,0)
 mat g1q1c=(-3.8782355,-1.7671008,.9215354,3.767286)
 mat g1q2a=(2.6124065,.9456864,0,0,0)
 mat g1q2c=(-4.7960932,-2.5118584,.5013833,3.4591656)
 mat g1q3a=(2.94338,0,.6065566,0,0)
 mat g1q3c=(-5.5042085,-3.321763,-.2910979,2.9441576)
 mat g1q4a=(2.1575548,1.3135908,0,0,0)
 mat g1q4c=(-3.6915979,-1.6576427,.9776602,3.5609863)
 mat g1q5a=(2.7242041,0,1.823336,0,0)
 mat g1q5c=(-6.2985207,-3.9930611,-1.0702193,2.299984)
 mat g1q6a=(2.6266401,0,1.2713007,0,0)
 mat g1q6c=(-5.2970614,-3.2486689,-.6403721,2.2178339)
 mat g1q7a=(2.1436959,1.3510367,0,0,0)
 mat g1q7c=(-4.1650335,-2.3304422,.0861086,2.7708504)
 mat g1q8a=(1.6588753,0,0,1.4824103,0)
 mat g1q8c=(-5.8700343,-3.9450868,-.889881,2.2294792)
 mat g1q9a=(1.4490391,0,0,1.097,0)
 mat g1q9c=(-4.854968,-3.1900848,-.6261177,1.9932152)
 mat g1q10a=(.9935323,0,0,0,1.4649552)
 mat g1q10c=(-3.4590475,-1.677168,.426561,2.8675858)
 mat g2m=(.0681382,-.0824054,.2943358,.1840323,.4392121)
 mat g2v=(1.0099854,.8838703999999999,1.2562739,.9213053,.3354192)
 mat g2q1a=(2.4110204,1.6608551,0,0,0)
 mat g2q1c=(-3.8782355,-1.7671008,.9215354,3.767286)
 mat g2q2a=(2.6124065,.9456864,0,0,0)
 mat g2q2c=(-4.7960932,-2.5118584,.5013833,3.4591656)
 mat g2q3a=(2.94338,0,.6065566,0,0)
 mat g2q3c=(-5.5042085,-3.321763,-.2910979,2.9441576)
 mat g2q4a=(2.1575548,1.3135908,0,0,0)
 mat g2q4c=(-3.6915979,-1.6576427,.9776602,3.5609863)
 mat g2q5a=(2.7242041,0,1.823336,0,0)
 mat g2q5c=(-6.2985207,-3.9930611,-1.0702193,2.299984)
 mat g2q6a=(2.6266401,0,1.2713007,0,0)
 mat g2q6c=(-5.2970614,-3.2486689,-.6403721,2.2178339)
 mat g2q7a=(2.1436959,1.3510367,0,0,0)
 mat g2q7c=(-4.1650335,-2.3304422,.0861086,2.7708504)
 mat g2q8a=(1.6588753,0,0,1.4824103,0)
 mat g2q8c=(-5.8700343,-3.9450868,-.889881,2.2294792)
 mat g2q9a=(1.4490391,0,0,1.097,0)
 mat g2q9c=(-4.854968,-3.1900848,-.6261177,1.9932152)
 mat g2q10a=(.9935323,0,0,0,1.4649552)
 mat g2q10c=(-3.4590475,-1.677168,.426561,2.8675858)
 mat g3m=(.4085444,.2285651,.6113262,.349839,.3778126)
 mat g3v=(1.292261,.9167938,2.2512673,2.6216329,1.5590249)
 mat g3q1a=(2.4110204,1.6608551,0,0,0)
 mat g3q1c=(-3.8782355,-1.7671008,.9215354,3.767286)
 mat g3q2a=(2.6124065,.9456864,0,0,0)
 mat g3q2c=(-4.7960932,-2.5118584,.5013833,3.4591656)
 mat g3q3a=(2.94338,0,.6065566,0,0)
 mat g3q3c=(-5.5042085,-3.321763,-.2910979,2.9441576)
 mat g3q4a=(2.1575548,1.3135908,0,0,0)
 mat g3q4c=(-3.6915979,-1.6576427,.9776602,3.5609863)
 mat g3q5a=(2.7242041,0,1.823336,0,0)
 mat g3q5c=(-6.2985207,-3.9930611,-1.0702193,2.299984)
 mat g3q6a=(2.6266401,0,1.2713007,0,0)
 mat g3q6c=(-5.2970614,-3.2486689,-.6403721,2.2178339)
 mat g3q7a=(2.1436959,1.3510367,0,0,0)
 mat g3q7c=(-4.1650335,-2.3304422,.0861086,2.7708504)
 mat g3q8a=(1.6588753,0,0,1.4824103,0)
 mat g3q8c=(-5.8700343,-3.9450868,-.889881,2.2294792)
 mat g3q9a=(1.4490391,0,0,1.097,0)
 mat g3q9c=(-4.854968,-3.1900848,-.6261177,1.9932152)
 mat g3q10a=(.9935323,0,0,0,1.4649552)
 mat g3q10c=(-3.4590475,-1.677168,.426561,2.8675858)
 mat g4m=(0,0,0,0,0)
 mat g4v=(1,1,1,1,1)
 mat g4q1a=(2.4110204,1.6608551,0,0,0)
 mat g4q1c=(-3.8782355,-1.7671008,.9215354,3.767286)
 mat g4q2a=(2.6124065,.9456864,0,0,0)
 mat g4q2c=(-4.7960932,-2.5118584,.5013833,3.4591656)
 mat g4q3a=(2.94338,0,.6065566,0,0)
 mat g4q3c=(-5.5042085,-3.321763,-.2910979,2.9441576)
 mat g4q4a=(2.1575548,1.3135908,0,0,0)
 mat g4q4c=(-3.6915979,-1.6576427,.9776602,3.5609863)
 mat g4q5a=(2.7242041,0,1.823336,0,0)
 mat g4q5c=(-6.2985207,-3.9930611,-1.0702193,2.299984)
 mat g4q6a=(2.6266401,0,1.2713007,0,0)
 mat g4q6c=(-5.2970614,-3.2486689,-.6403721,2.2178339)
 mat g4q7a=(2.1436959,1.3510367,0,0,0)
 mat g4q7c=(-4.1650335,-2.3304422,.0861086,2.7708504)
 mat g4q8a=(1.6588753,0,0,1.4824103,0)
 mat g4q8c=(-5.8700343,-3.9450868,-.889881,2.2294792)
 mat g4q9a=(1.4490391,0,0,1.097,0)
 mat g4q9c=(-4.854968,-3.1900848,-.6261177,1.9932152)
 mat g4q10a=(.9935323,0,0,0,1.4649552)
 mat g4q10c=(-3.4590475,-1.677168,.426561,2.8675858)
 loc qn 1
 loc qqlist
 foreach q in `qlist' {
  * turn each response into five columns in {0,1} to multiply probability
  tempvar q`qn'c1
  tempvar q`qn'c2
  tempvar q`qn'c3
  tempvar q`qn'c4
  tempvar q`qn'c5
  g byte `q`qn'c1'=(`q'==0)
  g byte `q`qn'c2'=(`q'==1)
  g byte `q`qn'c3'=(`q'==2)
  g byte `q`qn'c4'=(`q'==3)
  g byte `q`qn'c5'=(`q'==4)
  loc qqlist `qqlist' `q`qn'c1' `q`qn'c2' `q`qn'c3' `q`qn'c4' `q`qn'c5'
  loc qn=`qn'+1
  }
 *loop over groups
 if "`grouplist'"=="" {
  loc grouplist "1 2 3 4"
  }
 foreach g of numlist `grouplist' {
  qui count if `touse' & `gvar'==`g'
  if r(N)>0 {
   timer clear 1
   timer on 1
   mat FAmn=g`g'm
   mat FAva=g`g'v
   forv q=1/10 {
    mat FAq`q'a=g`g'q`q'a
    mat FAq`q'c=g`g'q`q'c
    }
   tempname tmp`g'
   g byte `tmp`g''=min(`touse',`gvar'==`g')
   mata:ipr_fa("`namelist'","`qqlist'","`tmp`g''",`intpoints',`zlimit',`GHQ')
   timer off 1
   qui timer list 1
   di as res "Group `g' complete in " r(t1) " seconds."
   }
  }
 qui replace `namelist'=round(`namelist',`round')
 if "`missok'"=="" qui replace `namelist'=. if `allmiss'==1
 ret scalar round=`round'
 ret scalar ghq=`GHQ'
 ret scalar intpoints=`intpoints'
 if "`ghq'"==""  ret scalar zlimit=`zlimit'
 ret local name "`namelist'"
end
version 10.2
mata:
void ipr_fa(string scalar f, string scalar x, string scalar tousename, real intpt, real zlim, real ghq) {
  st_view(y, ., tokens(f), tousename)
  q = st_data(., tokens(x), tousename)
  mean = st_matrix("FAmn")
  v = st_matrix("FAva")
  var=diag(v)
q1a= st_matrix("FAq1a")
q1c= st_matrix("FAq1c")
q2a= st_matrix("FAq2a")
q2c= st_matrix("FAq2c")
q3a= st_matrix("FAq3a")
q3c= st_matrix("FAq3c")
q4a= st_matrix("FAq4a")
q4c= st_matrix("FAq4c")
q5a= st_matrix("FAq5a")
q5c= st_matrix("FAq5c")
q6a= st_matrix("FAq6a")
q6c= st_matrix("FAq6c")
q7a= st_matrix("FAq7a")
q7c= st_matrix("FAq7c")
q8a= st_matrix("FAq8a")
q8c= st_matrix("FAq8c")
q9a= st_matrix("FAq9a")
q9c= st_matrix("FAq9c")
q10a=st_matrix("FAq10a")
q10c=st_matrix("FAq10c")
P=J(rows(q),1,0)
Pf=J(rows(q),1,0)
Pft=J(rows(q),1,0)
 Z=_gauss_hermite_nodes(intpt)
 for (i=1; i<=intpt; i++) {
  if (ghq==1) t1=Z[1,i]
  else t1=(i-1)*2*zlim/(intpt-1)-(zlim)
  for (j=1; j<=intpt; j++) {
   if (ghq==1) t2=Z[1,j]
   else t2=(j-1)*2*zlim/(intpt-1)-(zlim)
   for (k=1; k<=intpt; k++) {
    if (ghq==1) t3=Z[1,k]
    else t3=(k-1)*2*zlim/(intpt-1)-(zlim)
    for (l=1; l<=intpt; l++) {
     if (ghq==1) t4=Z[1,l]
     else t4=(l-1)*2*zlim/(intpt-1)-(zlim)
     for (n=1; n<=intpt; n++) {
      if (ghq==1) t5=Z[1,n]
      else t5=(n-1)*2*zlim/(intpt-1)-(zlim)
 theta=(t1,t2,t3,t4,t5)
 a1theta=theta*q1a'
 a2theta=theta*q2a'
 a3theta=theta*q3a'
 a4theta=theta*q4a'
 a5theta=theta*q5a'
 a6theta=theta*q6a'
 a7theta=theta*q7a'
 a8theta=theta*q8a'
 a9theta=theta*q9a'
 a10theta=theta*q10a'
 lp1=(ln(1-invlogit(a1theta-q1c[1]))\ln(invlogit(a1theta-q1c[1])-invlogit(a1theta-q1c[2]))\ln(invlogit(a1theta-q1c[2])-invlogit(a1theta-q1c[3]))\ln(invlogit(a1theta-q1c[3])-invlogit(a1theta-q1c[4]))\ln(invlogit(a1theta-q1c[4])))
 lp2=(ln(1-invlogit(a2theta-q2c[1]))\ln(invlogit(a2theta-q2c[1])-invlogit(a2theta-q2c[2]))\ln(invlogit(a2theta-q2c[2])-invlogit(a2theta-q2c[3]))\ln(invlogit(a2theta-q2c[3])-invlogit(a2theta-q2c[4]))\ln(invlogit(a2theta-q2c[4])))
 lp3=(ln(1-invlogit(a3theta-q3c[1]))\ln(invlogit(a3theta-q3c[1])-invlogit(a3theta-q3c[2]))\ln(invlogit(a3theta-q3c[2])-invlogit(a3theta-q3c[3]))\ln(invlogit(a3theta-q3c[3])-invlogit(a3theta-q3c[4]))\ln(invlogit(a3theta-q3c[4])))
 lp4=(ln(1-invlogit(a4theta-q4c[1]))\ln(invlogit(a4theta-q4c[1])-invlogit(a4theta-q4c[2]))\ln(invlogit(a4theta-q4c[2])-invlogit(a4theta-q4c[3]))\ln(invlogit(a4theta-q4c[3])-invlogit(a4theta-q4c[4]))\ln(invlogit(a4theta-q4c[4])))
 lp5=(ln(1-invlogit(a5theta-q5c[1]))\ln(invlogit(a5theta-q5c[1])-invlogit(a5theta-q5c[2]))\ln(invlogit(a5theta-q5c[2])-invlogit(a5theta-q5c[3]))\ln(invlogit(a5theta-q5c[3])-invlogit(a5theta-q5c[4]))\ln(invlogit(a5theta-q5c[4])))
 lp6=(ln(1-invlogit(a6theta-q6c[1]))\ln(invlogit(a6theta-q6c[1])-invlogit(a6theta-q6c[2]))\ln(invlogit(a6theta-q6c[2])-invlogit(a6theta-q6c[3]))\ln(invlogit(a6theta-q6c[3])-invlogit(a6theta-q6c[4]))\ln(invlogit(a6theta-q6c[4])))
 lp7=(ln(1-invlogit(a7theta-q7c[1]))\ln(invlogit(a7theta-q7c[1])-invlogit(a7theta-q7c[2]))\ln(invlogit(a7theta-q7c[2])-invlogit(a7theta-q7c[3]))\ln(invlogit(a7theta-q7c[3])-invlogit(a7theta-q7c[4]))\ln(invlogit(a7theta-q7c[4])))
 lp8=(ln(1-invlogit(a8theta-q8c[1]))\ln(invlogit(a8theta-q8c[1])-invlogit(a8theta-q8c[2]))\ln(invlogit(a8theta-q8c[2])-invlogit(a8theta-q8c[3]))\ln(invlogit(a8theta-q8c[3])-invlogit(a8theta-q8c[4]))\ln(invlogit(a8theta-q8c[4])))
 lp9=(ln(1-invlogit(a9theta-q9c[1]))\ln(invlogit(a9theta-q9c[1])-invlogit(a9theta-q9c[2]))\ln(invlogit(a9theta-q9c[2])-invlogit(a9theta-q9c[3]))\ln(invlogit(a9theta-q9c[3])-invlogit(a9theta-q9c[4]))\ln(invlogit(a9theta-q9c[4])))
 lp10=(ln(1-invlogit(a10theta-q10c[1]))\ln(invlogit(a10theta-q10c[1])-invlogit(a10theta-q10c[2]))\ln(invlogit(a10theta-q10c[2])-invlogit(a10theta-q10c[3]))\ln(invlogit(a10theta-q10c[3])-invlogit(a10theta-q10c[4]))\ln(invlogit(a10theta-q10c[4])))
 lp=lp1\lp2\lp3\lp4\lp5\lp6\lp7\lp8\lp9\lp10
 P=exp(cross(q',lp):+((-.5*(theta-mean)*invsym(var)*(theta-mean)')-ln(pi()*det(var))))
 Pf=Pf+P
 Pft=Pft+P:*t1
     }
    }
   }
  }
 }
y[.,.]=Pft:/Pf:*15:+50
}
end

exit