(* Mathematica Routines for Statistical Models for Behavioral Observations David Rogosa and Ghassan Ghandour November 1991 Copyright (c) 1991 This file contains Mathematica Expressions Corresponding to most of the equations in Rogosa and Ghandour 1991, Statistical Models for Behavioral Observations *) (* Some Notational Compromises Equation numbering: Equation 1.2' p.174 is indicated as eq1p2prime , and so forth Conditioning: \$ stands in for the | character in conditional expections and variances (e.g. eq1p3 below) A rendition of the Glossary pp. 158-9, adapted for the notation of these routines, is appended to the end of this file in comment form. *) (* Useful Constants Vectors To keep in touch with your settings for parameter values, we provide 4 vectors of commonly used values that can be called up during calculations-- indConst, popConst, recConst, timeConst *) indConst := {"alpha = ",alpha, "inc = ",inc, "phi = ",phi, "Iprocess = ",Iprocess, "mu1 = ",mu1, "tau2 = ",tau2, "T = ",T} popConst := {"mualpha = ",mualpha, "sig2alpha = ",sig2alpha, "muinc = ",muinc, "sig2inc = ",sig2inc, "corrincalpha = ",corrincalpha, "Iprocess = ",Iprocess, "muphi = ",muphi,"sig2phi = ",sig2phi, "corrincphi = ",corrincphi, "mumu1 = ",mumu1, "sig2mu1 = ",sig2mu1} recConst := {"R = ",R, "kap = ",kap, "gam = ",gam, "om = ",om, "psi = ", psi, "minkap = ",minkap, "rangekap = ",rangekap, "mingam = ",mingam, "rangegam = ",rangegam, "corrkapgam = ",corrkapgam, "muom = ",muom, "sig2om = ",sig2om, "mupsi = ",mupsi, "mupsif = ",mupsif} timeConst := {"occ = ",occ, "T = ",T, "s = ",s, "s2 = ",s2, "al = ",al, "tau2 = ",tau2, "deltaInMinutes = ",delta*60, "segs = ",T/delta} eq1p2prime ::= vpT vpT := (inc^3 sig2w)/T + (1/6 + (inc^4 sig2w^2)/2 - (inc^3 mu3w)/3)/(T^2) (* For some calculations use a function form for vpT *) funcvpT[T_, inc_, sig2w_, mu3w_] := (inc^3 sig2w)/T + (1/6 + (inc^4 sig2w^2)/2 - (inc^3 mu3w)/3)/(T^2) (* For gamma point process sig2w is 1/(alpha inc^2) mu3w is 2/(alpha^2 inc^3) *) (* For alternating Poisson sig2w is (1 - 2*phi + 2*phi^2)/(inc^2) mu3w is (2 - 6*phi + 6*phi^2)/(inc^3) *) (* Recorder parameters and derived quantities, single recorder *) om := 1 - kap + gam psi := kap*(1 - kap) + gam (* Moments for a fixed recorder *) eq1p3a ::= EMprDivT\$pr EMprDivT\$pr := inc om eq1p3b ::= VMprDivT\$pr VMprDivT\$pr := (inc psi)/T + om^2 vpT (* Recorder parameters and derived quantities for a population of recorders, assuming uniform distributions for kap and gam *) (* Input Recorder population parameters *) minkap = . rangekap = . mingam = . rangegam = . corrkapgam = . (* Standard Omega/Psi Moments from kappa gamma Uniform *) Clear[muunf] Clear[varunf] Clear[meanpsi] Clear[meanom] Clear[varomg] (* Uniform Moments *) muunf[min_, range_] := (range + 2 min) / 2 varunf[range_] := range^2 / 12 (* Theoretical Moments Page 176 Applied to Uniform Distributions *) meanpsi[minkap_, rangekap_, mingam_, rangegam_] := muunf[minkap, rangekap] - varunf[rangekap] - muunf[minkap, rangekap]^2 + muunf[mingam, rangegam] meanom[minkap_, rangekap_, mingam_, rangegam_] := 1 - muunf[minkap, rangekap] + muunf[mingam, rangegam] varom[rangekap_, rangegam_, corrkapgam_] := varunf[rangekap] + varunf[rangegam] - 2 corrkapgam Sqrt[varunf[rangekap] varunf[rangegam]] muom := meanom[minkap, rangekap, mingam, rangegam] sig2om := varom[rangekap, rangegam,corrkapgam] mupsi := meanpsi[minkap, rangekap, mingam, rangegam] (* For fixed recorder we use muspif to represent degenerate kappa gamma distribution *) mupsif := meanpsi[minkap, rangekap, mingam, rangegam] + varunf[rangekap] (* Equation 1.4 *) eq1p4a ::= EMprDivT\$p EMprDivT\$p := inc muom eq1p4b ::= VMprDivT\$p VMprDivT\$p := (inc mupsi/T) + (muom^2 + sig2om)*vpT + inc^2 sig2om (* Equation 1.6 *) eq1p6a ::= EMbarpDivT\$p EMbarpDivT\$p := inc muom eq1p6b ::= VMbarpDivT\$p VMbarpDivT\$p := (muom^2 vpT) + (1/R)((inc mupsi/T) + sig2om*(inc^2 + vpT)) eq1p7b ::= VNpoDivT\$p VNpoDivT\$p := tau2 + vpT (* Defining Observation Schedule S and Allocation Matrix A *) Clear[kp] Clear[alKeep] Clear[alDiag] Clear[s] Clear[one] Clear[vpos] Clear[saatst] Clear[sst] Clear[svt] Clear[svt1] Clear[svt1st] Clear[Tplus] (* Note--- dimension of s, i.e. number of occasions O = occ needs to be specified before computation *) occ = 1 (* Allocation matrix for random recorder rosters resampled each Occ *) ADiag := IdentityMatrix[Round[occ]] (* AKeep is allocation matrix for keeping same recorder roster see pages 178, 180 *) kp[ik_,jk_] := If[jk == 1, 1, 0, 0] AKeep := Table[kp[iak, jak], {iak, Round[occ]}, {jak, Round[occ]}] (* Initial Setting--Allocation Matrix resampled recorder *) al := ADiag (* Observation schedule vectors *) (* for equal observation intervals set s := Table[T, {is, Round[occ]}] *) s = {T} s2 := Table[s[[iss]]^2, {iss, Round[occ]}] (* Vpo(s) vector *) vpos := Table[funcvpT[s[[iv]], inc, sig2w, mu3w],{iv, Round[occ]}] (* Matrix products *) saatst := s.al.Transpose[al].s sst := s.s (* Tp+ *) one := Table[i1/i1, {i1, Round[occ]}] Tplus := s.one (* eq1p8b = VNpplusDivTplus\$p is obtained from eq1p11b with single perfect recorder-- minkap = rangekap = mingam = rangegam = 0 -- and R = 1 *) (* eq1p9b = VMpplusrDivTplus\$p is obtained from eq1p11b with single recorder, R = 1 *) (* Equation 1.11: Multiple Recorders, Multiple Occasions *) eq1p11a ::= EMbarpplusDivTplus\$p EMbarpplusDivTplus\$p := inc muom eq1p11b ::= VMbarpplusDivTplus\$p VMbarpplusDivTplus\$p := (inc^2*sig2om*saatst)/(R*Tplus^2) + (inc*mupsi)/(R*Tplus) + (R muom^2 + sig2om) * (sst*tau2 + s2.vpos)/(R*Tplus^2) (* Section 1.3-- Examples of Renewal Processes *) eq1p12 ::= vpTpoisson vpTpoisson := funcvpT[T, inc, 1/inc^2, 2/inc^3] eq1p13 ::= vpTgamma vpTgamma := funcvpT[T, inc, 1/(alpha inc^2), 2/(alpha^2 inc^3)] eq1p14 ::= vpTaltP vpTaltP := funcvpT[T, inc, (1 - 2*phi + 2*phi^2)/(inc^2) , (2 - 6*phi + 6*phi^2)/(inc^3)] (* SECTION 2 *) (* input for parameter distributions over population of individuals *) muinc = . sig2inc = . muphi = . sig2phi = . corrincphi = . covincphi := corrincphi Sqrt[sig2phi sig2inc] mualpha = . sig2alpha = . corrincalpha = . covincalpha := corrincalpha Sqrt[sig2alpha sig2inc] eq2p4 ::= relNDivT relNDivT := sig2inc/(sig2inc + EpvpT) eq2p6 ::= EpvpTgamma EpvpTgamma := (-mualpha^2 + mualpha^4 - 6*covincalpha*mualpha^2*T + 6*mualpha^3*muinc*T - 3*sig2alpha + 6*mualpha*muinc*T*sig2alpha)/(6*mualpha^4*T^2) eq2p7 ::= EpvpTaltP EpvpTaltP := (2*muphi^2 - 4*muphi^3 + 2*muphi^4 - 2*covincphi*T + muinc*T + 4*covincphi*muphi*T - 2*muinc*muphi*T + 2*muinc*muphi^2*T + 2*sig2phi - 12*muphi*sig2phi + 12*muphi^2*sig2phi + 2*muinc*T*sig2phi)/T^2 (* In funcEpvpT declare Iprocess = 1 for gamma; 0 for Alternating Poisson *) funcEpvpT[T_, muinc_, sig2inc_, muphi_, sig2phi_, covincphi_, mualpha_, sig2alpha_, covincalpha_, Iprocess_] := Iprocess*(-mualpha^2 + mualpha^4 - 6*covincalpha*mualpha^2*T + 6*mualpha^3*muinc*T - 3*sig2alpha + 6*mualpha*muinc*T*sig2alpha)/(6*mualpha^4*T^2) + (1 - Iprocess)*(2*muphi^2 - 4*muphi^3 + 2*muphi^4 - 2*covincphi*T + muinc*T + 4*covincphi*muphi*T - 2*muinc*muphi*T + 2*muinc*muphi^2*T + 2*sig2phi - 12*muphi*sig2phi + 12*muphi^2*sig2phi + 2*muinc*T*sig2phi)/T^2 eq2p8 := sig2inc/((muinc/muT)*(1 + (sig2T/muT^2) - covincT/(muinc muT)) + sig2inc) eq2p9 ::= relMprDivT relMprDivT := sig2inc/(sig2inc + (muinc mupsi)/(T muom^2) + EpvpT + (sig2om/muom^2)*(EpvpT + sig2inc + muinc^2)) (* Eq2.10b is 2.12b for occ = 1, R = 1 *) (* Moments over the population for Multiple Occasions *) Epvpos := Table[funcEpvpT[s[[iv]], muinc, sig2inc, muphi, sig2phi, covincphi, mualpha, sig2alpha, covincalpha, Iprocess],{iv, Round[occ]}] eq2p12 ::= EpVMbarpplusDivTplus EpVMbarpplusDivTplus := (muinc*mupsi)/(R*Tplus) + (muinc^2 + sig2inc)*sig2om*saatst/ (R*Tplus^2) + (R muom^2 + sig2om) * (sst*tau2 + s2.Epvpos)/(R*Tplus^2) (* Tables 3 through 5-- Spearman-Brown Illustrations *) actualReliability := sig2inc/(sig2inc + eq2p12/muom^2) (* Compute seedReliability for k = 1 for use in Spearman-Brown *) SBReliability := k*seedReliability/((k - 1)*seedReliability + 1) (* Note in Table 5, bottom leftmost entry should be .312, not .296 *) eq3p3 ::= rhoMM' rhoMM' := (sig2inc + EpvpT)/(sig2inc + (muinc mupsi)/(T muom^2) + EpvpT + (sig2om/muom^2)*(EpvpT + sig2inc + muinc^2)) eq3p8 ::= prVC prVC := {prVCp, prVCr, prVCpr, prVCeps} prVCp := muom^2 sig2inc + muom^2 EpvpT prVCr := sig2om muinc^2 prVCpr := sig2inc sig2om + sig2om EpvpT prVCeps := muinc mupsi/T (* pxoxr G theory quantities *) eq3p11 ::= porVC porVC := {porVCp, porVCo, porVCr, porVCpo, porVCpr, porVCro, porVCpor, porVCeps} porVCp := muom^2 sig2inc porVCo = 0 porVCr := sig2om muinc^2 porVCpo := muom^2*(EpvpT + tau2) porVCpr := sig2inc sig2om porVCro = 0 porVCpor := sig2om*(EpvpT + tau2) porVCeps := muinc mupsi/T eq3p12 ::= porVCabs porVCabs := porVCo + porVCr + porVCpo + porVCpr + porVCro + (porVCpor + porVCeps) eq3p13 ::= porVCrel porVCrel := porVCpo + porVCpr + (porVCpor + porVCeps) eq3p15 ::= Varpdel\$r Varpdel\$r := sig2inc*(1 - muom)^2 + (muinc mupsif)/T + muom^2 * (EpvpT + tau2) (* D study projections-- Multiple Occasions and Recorders *) eq3p16 ::= porVCabsOR porVCabsOR := porVCo/occ + porVCr/R + porVCpo/occ + porVCpr/R + (porVCro + (porVCpor + porVCeps))/(occ R) eq3p17 ::= porVCrelOR porVCrelOR := porVCo/occ + porVCpo/occ + porVCpr/R + (porVCro + (porVCpor + porVCeps))/(occ R) (* Eq3.18a is 2.12 with al = ADiag and Eq3.18b is 2.12 with al = AKeep *) eq3p19 ::= VarpdelOR\$scriptR VarpdelOR\$scriptR := sig2inc*(1 - muom)^2 + (muinc mupsif)/ (R T occ) + muom^2*(EpvpT + tau2)/occ (* Section 4 *) (* Equation 4.1 *) eq4p1 ::= EC EC := T muc inc (* Equation 4.2 *) eq4p2 ::= VC VC := T inc sig2c*(1 - rhocw^2) + T^2*funcvpT[T, inc, sig2w, mu3w]* (muc - betacw/inc)^2 (* Equation 4.3 *) eq4p3 ::= CovCN CovCN := T^2*(muc - betacw/inc)*funcvpT[T, inc, sig2w, mu3w] (* Equation 4.4a *) eq4p4a ::= ECDivN ECDivN := muc + (betacw*funcvpT[T, inc, sig2w, mu3w])/inc^3 (* Equation 4.4b *) eq4p4b ::= VCDivN VCDivN := (sig2c/(T inc))*((1 - rhocw^2) + T* funcvpT[T, inc, sig2w, mu3w]*(rhocw^2/(sig2w inc^3))) (* For proportions muc = pi sig2c = pi*(1-pi) rhocw = rho betacw = rhocw*Sqrt[sig2c/sig2w] *) (* Equation 4.5 is obtained from 4.4b using the proportion settings for moments of c,w *) eq4p5 ::= VZDivN\$p VZDivN\$p := -((-1 + pi)*pi*(rho^2 - 2*inc^3*mu3w*rho^2 + 6*T*inc^3*sig2w + 3*inc^4*rho^2*sig2w^2))/ (6*T^2*inc^4*sig2w) (* Equation 4.6 can be expressed more generally for either an Alternating Poisson Process or for a Gamma Process Here we give 4.6 for the Gamma Process with alphap = alpha for all p Either one can be reduced to Poisson case given in p.228 *) eq4p6Gam ::= EpvpZDivNG EpvpZDivNG := (2*covincpi*muinc*rho^2 - 2*alpha^2*covincpi*muinc*rho^2 - 4*covincpi*muinc*mupi*rho^2 + 4*alpha^2*covincpi*muinc*mupi*rho^2 - muinc^2*mupi*rho^2 + alpha^2*muinc^2*mupi*rho^2 + muinc^2*mupi^2*rho^2 - alpha^2*muinc^2*mupi^2*rho^2 - 3*mupi*rho^2*sig2inc + 3*alpha^2*mupi*rho^2*sig2inc + 3*mupi^2*rho^2*sig2inc - 3*alpha^2*mupi^2*rho^2*sig2inc + muinc^2*rho^2*sig2pi - alpha^2*muinc^2*rho^2*sig2pi - 6*alpha*covincpi*muinc^2*T + 12*alpha*covincpi*muinc^2*mupi*T + 6*alpha*muinc^3*mupi*T - 6*alpha*muinc^3*mupi^2*T + 6*alpha*muinc*mupi*sig2inc*T - 6*alpha*muinc*mupi^2*sig2inc*T - 6*alpha*muinc^3*sig2pi*T)/ (6*alpha*muinc^4*T^2) (* For Prevalence and Event Duration From Alternating Poisson muc := phi/inc sig2c := muc^2 betacw := phi^2 / ((1-phi)^2 + phi^2) rhocw := Sqrt[betacw] *) eq4p7 ::= VLDivT\$p VLDivT\$p := ((phi^2*(1 - phi^2/((1 - phi)^2 + phi^2))*T)/inc + (phi/inc - phi^2/(inc*((1 - phi)^2 + phi^2)))^2* ((2*(1 - phi)^2*phi^2)/T^2 + (inc*((1 - phi)^2 + phi^2))/T)*T^2)/T^2 (* Section 4.4.3, Psychometric Properties of Empirical Prevalence *) page231Epvpprev := (2*T*muinc^3*muphi^2 - 20*T*muinc^3*muphi^3 + 2*muinc^2*muphi^4 + 98*T*muinc^3*muphi^4 - 24*muinc^2*muphi^5 - 304*T*muinc^3*muphi^5 + 132*muinc^2*muphi^6 + 656*T*muinc^3*muphi^6 - 440*muinc^2*muphi^7 - 1024*T*muinc^3*muphi^7 + 986*muinc^2*muphi^8 + 1168*T*muinc^3*muphi^8 - 1552*muinc^2*muphi^9 - 960*T*muinc^3*muphi^9 + 1736*muinc^2*muphi^10 + 544*T*muinc^3*muphi^10 - 1360*muinc^2*muphi^11 - 192*T*muinc^3*muphi^11 + 712*muinc^2*muphi^12 + 32*T*muinc^3*muphi^12 - 224*muinc^2*muphi^13 + 32*muinc^2*muphi^14 + 2*T*muinc*muphi^2*sig2inc - 20*T*muinc*muphi^3*sig2inc + 6*muphi^4*sig2inc + 98*T*muinc*muphi^4*sig2inc - 72*muphi^5*sig2inc - 304*T*muinc*muphi^5*sig2inc + 396*muphi^6*sig2inc + 656*T*muinc*muphi^6*sig2inc - 1320*muphi^7*sig2inc - 1024*T*muinc*muphi^7*sig2inc + 2958*muphi^8*sig2inc + 1168*T*muinc*muphi^8*sig2inc - 4656*muphi^9*sig2inc - 960*T*muinc*muphi^9*sig2inc + 5208*muphi^10*sig2inc + 544*T*muinc*muphi^10*sig2inc - 4080*muphi^11*sig2inc - 192*T*muinc*muphi^11*sig2inc + 2136*muphi^12*sig2inc + 32*T*muinc*muphi^12*sig2inc - 672*muphi^13*sig2inc + 96*muphi^14*sig2inc - 4*T*corrincphi*muinc^2*muphi*sig2inc^(1/2)* sig2phi^(1/2) + 44*T*corrincphi*muinc^2*muphi^2*sig2inc^(1/2)* sig2phi^(1/2) - 16*corrincphi*muinc*muphi^3*sig2inc^(1/2)* sig2phi^(1/2) - 232*T*corrincphi*muinc^2*muphi^3*sig2inc^(1/2)* sig2phi^(1/2) + 208*corrincphi*muinc*muphi^4*sig2inc^(1/2)* sig2phi^(1/2) + 768*T*corrincphi*muinc^2*muphi^4*sig2inc^(1/2)* sig2phi^(1/2) - 1200*corrincphi*muinc*muphi^5*sig2inc^(1/2)* sig2phi^(1/2) - 1760*T*corrincphi*muinc^2*muphi^5*sig2inc^(1/2)* sig2phi^(1/2) + 4112*corrincphi*muinc*muphi^6*sig2inc^(1/2)* sig2phi^(1/2) + 2912*T*corrincphi*muinc^2*muphi^6*sig2inc^(1/2)* sig2phi^(1/2) - 9376*corrincphi*muinc*muphi^7*sig2inc^(1/2)* sig2phi^(1/2) - 3520*T*corrincphi*muinc^2*muphi^7*sig2inc^(1/2)* sig2phi^(1/2) + 14976*corrincphi*muinc*muphi^8*sig2inc^(1/2)* sig2phi^(1/2) + 3072*T*corrincphi*muinc^2*muphi^8*sig2inc^(1/2)* sig2phi^(1/2) - 17056*corrincphi*muinc*muphi^9*sig2inc^(1/2)* sig2phi^(1/2) - 1856*T*corrincphi*muinc^2*muphi^9*sig2inc^(1/2)* sig2phi^(1/2) + 13728*corrincphi*muinc*muphi^10*sig2inc^(1/2)* sig2phi^(1/2) + 704*T*corrincphi*muinc^2*muphi^10*sig2inc^(1/2)* sig2phi^(1/2) - 7488*corrincphi*muinc*muphi^11*sig2inc^(1/2)* sig2phi^(1/2) - 128*T*corrincphi*muinc^2*muphi^11*sig2inc^(1/2)* sig2phi^(1/2) + 2496*corrincphi*muinc*muphi^12*sig2inc^(1/2)* sig2phi^(1/2) - 384*corrincphi*muinc*muphi^13*sig2inc^(1/2)* sig2phi^(1/2) + 2*T*muinc^3*sig2phi - 28*T*muinc^3*muphi*sig2phi + 12*muinc^2*muphi^2*sig2phi + 172*T*muinc^3*muphi^2*sig2phi - 176*muinc^2*muphi^3*sig2phi - 640*T*muinc^3*muphi^3*sig2phi + 1084*muinc^2*muphi^4*sig2phi + 1616*T*muinc^3*muphi^4*sig2phi - 3832*muinc^2*muphi^5*sig2phi - 2912*T*muinc^3*muphi^5*sig2phi + 8888*muinc^2*muphi^6*sig2phi + 3808*T*muinc^3*muphi^6*sig2phi - 14464*muinc^2*muphi^7*sig2phi - 3584*T*muinc^3*muphi^7*sig2phi + 16936*muinc^2*muphi^8*sig2phi + 2336*T*muinc^3*muphi^8*sig2phi - 14160*muinc^2*muphi^9*sig2phi - 960*T*muinc^3*muphi^9*sig2phi + 8112*muinc^2*muphi^10*sig2phi + 192*T*muinc^3*muphi^10*sig2phi - 2880*muinc^2*muphi^11*sig2phi + 480*muinc^2*muphi^12*sig2phi)/ (T^2*muinc^4*(1 - 2*muphi + 2*muphi^2)^4) relprev := sig2phi/(page231Epvpprev + sig2phi) (* Section 4.5 Event Duration *) (* parameter value *) mu1 := phi/inc page232biasDur ::= biasLDivN\$p biasLDivN\$p := (phi^2*(2*phi^2 - 4*phi^3 + 2*phi^4 + inc*T - 2*inc*phi*T + 2*inc*phi^2*T))/(inc^3*(1 - 2*phi + 2*phi^2)*T^2) eq4p8 ::= VLDivN\$p VLDivN\$p := (phi^2*(2*phi^4 - 4*phi^5 + 2*phi^6 + inc*T - 4*inc*phi*T + 8*inc*phi^2*T - 8*inc*phi^3*T + 4*inc*phi^4*T))/(inc^4*(1 - 2*phi + 2*phi^2)^2*T^2) page232Epvpdur := (T*muinc^3*muphi^2 - 8*T*muinc^3*muphi^3 + 32*T*muinc^3*muphi^4 - 80*T*muinc^3*muphi^5 + 2*muinc^2*muphi^6 + 136*T*muinc^3*muphi^6 - 12*muinc^2*muphi^7 - 160*T*muinc^3*muphi^7 + 34*muinc^2*muphi^8 + 128*T*muinc^3*muphi^8 - 56*muinc^2*muphi^9 - 64*T*muinc^3*muphi^9 + 56*muinc^2*muphi^10 + 16*T*muinc^3*muphi^10 - 32*muinc^2*muphi^11 + 8*muinc^2*muphi^12 + 6*T*muinc*muphi^2*sig2inc - 48*T*muinc*muphi^3*sig2inc + 192*T*muinc*muphi^4*sig2inc - 480*T*muinc*muphi^5*sig2inc + 20*muphi^6*sig2inc + 816*T*muinc*muphi^6*sig2inc - 120*muphi^7*sig2inc - 960*T*muinc*muphi^7*sig2inc + 340*muphi^8*sig2inc + 768*T*muinc*muphi^8*sig2inc - 560*muphi^9*sig2inc - 384*T*muinc*muphi^9*sig2inc + 560*muphi^10*sig2inc + 96*T*muinc*muphi^10*sig2inc - 320*muphi^11*sig2inc + 80*muphi^12*sig2inc - 6*T*corrincphi*muinc^2*muphi*sig2inc^(1/2)* sig2phi^(1/2) + 48*T*corrincphi*muinc^2*muphi^2*sig2inc^(1/2)* sig2phi^(1/2) - 192*T*corrincphi*muinc^2*muphi^3*sig2inc^(1/2)* sig2phi^(1/2) + 480*T*corrincphi*muinc^2*muphi^4*sig2inc^(1/2)* sig2phi^(1/2) - 48*corrincphi*muinc*muphi^5*sig2inc^(1/2)* sig2phi^(1/2) - 816*T*corrincphi*muinc^2*muphi^5*sig2inc^(1/2)* sig2phi^(1/2) + 272*corrincphi*muinc*muphi^6*sig2inc^(1/2)* sig2phi^(1/2) + 960*T*corrincphi*muinc^2*muphi^6*sig2inc^(1/2)* sig2phi^(1/2) - 704*corrincphi*muinc*muphi^7*sig2inc^(1/2)* sig2phi^(1/2) - 768*T*corrincphi*muinc^2*muphi^7*sig2inc^(1/2)* sig2phi^(1/2) + 1056*corrincphi*muinc*muphi^8*sig2inc^(1/2)* sig2phi^(1/2) + 384*T*corrincphi*muinc^2*muphi^8*sig2inc^(1/2)* sig2phi^(1/2) - 960*corrincphi*muinc*muphi^9*sig2inc^(1/2)* sig2phi^(1/2) - 96*T*corrincphi*muinc^2*muphi^9*sig2inc^(1/2)* sig2phi^(1/2) + 512*corrincphi*muinc*muphi^10*sig2inc^(1/2)* sig2phi^(1/2) - 128*corrincphi*muinc*muphi^11*sig2inc^(1/2)* sig2phi^(1/2) + T*muinc^3*sig2phi - 8*T*muinc^3*muphi*sig2phi + 32*T*muinc^3*muphi^2*sig2phi - 80*T*muinc^3*muphi^3*sig2phi + 30*muinc^2*muphi^4*sig2phi + 136*T*muinc^3*muphi^4*sig2phi - 156*muinc^2*muphi^5*sig2phi - 160*T*muinc^3*muphi^5*sig2phi + 344*muinc^2*muphi^6*sig2phi + 128*T*muinc^3*muphi^6*sig2phi - 448*muinc^2*muphi^7*sig2phi - 64*T*muinc^3*muphi^7*sig2phi + 376*muinc^2*muphi^8*sig2phi + 16*T*muinc^3*muphi^8*sig2phi - 192*muinc^2*muphi^9*sig2phi + 48*muinc^2*muphi^10*sig2phi)/ (T^2*muinc^6*(1 - 2*muphi + 2*muphi^2)^4) (* Moments of Event Duration for the {Inc, Phi} distribution *) mumu1 := muphi/muinc - (1/muinc^2)*covincphi + (muphi/muinc^3)*sig2inc sig2mu1 := (muphi/muinc)^2*(sig2phi/muphi^2 + sig2inc/muinc^2 - 2 covincphi/(muphi muinc)) reldur := sig2mu1/(page232Epvpdur + sig2mu1) (* Section 4.6 Relative Precisions *) (* CV for individual p , Section 4.6.1 *) incCV\$p := Sqrt[vpT]/inc prevCV\$p := Sqrt[eq4p7]/phi durCV\$p := Sqrt[eq4p8]/mu1 (* Section 5 *) (* define segs = T/delta *) eq5p1 ::= VelDivsegs\$p VelDivsegs\$p := (delta phi*(1 - phi)/T)* (2/(1 - Exp[-inc delta/(phi*(1 - phi))]) - 1) (* Entries for Table 8 *) table8entries := ((-1 + 2/(1 - E^((delta*inc)/((-1 + phi)*phi))))*delta* inc^2*(1 - phi)*(1 - 2*phi + 2*phi^2)^2*T)/ (2*(-1 + phi)^2*phi*(phi^2 - 6*phi^3 + 13*phi^4 - 12*phi^5 + 4*phi^6 + inc*T - 4*inc*phi*T + 8*inc*phi^2*T - 8*inc*phi^3*T + 4*inc*phi^4*T)) (* Ratio of Instantaneous to Continuous reliabilities, page 237 *) page237RelRatio := (sig2phi + (2*T*muinc^3*muphi^2 - 20*T*muinc^3*muphi^3 + 2*muinc^2*muphi^4 + 98*T*muinc^3*muphi^4 - 24*muinc^2*muphi^5 - 304*T*muinc^3*muphi^5 + 132*muinc^2*muphi^6 + 656*T*muinc^3*muphi^6 - 440*muinc^2*muphi^7 - 1024*T*muinc^3*muphi^7 + 986*muinc^2*muphi^8 + 1168*T*muinc^3*muphi^8 - 1552*muinc^2*muphi^9 - 960*T*muinc^3*muphi^9 + 1736*muinc^2*muphi^10 + 544*T*muinc^3*muphi^10 - 1360*muinc^2*muphi^11 - 192*T*muinc^3*muphi^11 + 712*muinc^2*muphi^12 + 32*T*muinc^3*muphi^12 - 224*muinc^2*muphi^13 + 32*muinc^2*muphi^14 + 2*T*muinc*muphi^2*sig2inc - 20*T*muinc*muphi^3*sig2inc + 6*muphi^4*sig2inc + 98*T*muinc*muphi^4*sig2inc - 72*muphi^5*sig2inc - 304*T*muinc*muphi^5*sig2inc + 396*muphi^6*sig2inc + 656*T*muinc*muphi^6*sig2inc - 1320*muphi^7*sig2inc - 1024*T*muinc*muphi^7*sig2inc + 2958*muphi^8*sig2inc + 1168*T*muinc*muphi^8*sig2inc - 4656*muphi^9*sig2inc - 960*T*muinc*muphi^9*sig2inc + 5208*muphi^10*sig2inc + 544*T*muinc*muphi^10*sig2inc - 4080*muphi^11*sig2inc - 192*T*muinc*muphi^11*sig2inc + 2136*muphi^12*sig2inc + 32*T*muinc*muphi^12*sig2inc - 672*muphi^13*sig2inc + 96*muphi^14*sig2inc - 4*T*corrincphi*muinc^2*muphi*sig2inc^(1/2)* sig2phi^(1/2) + 44*T*corrincphi*muinc^2*muphi^2*sig2inc^(1/2)* sig2phi^(1/2) - 16*corrincphi*muinc*muphi^3*sig2inc^(1/2)* sig2phi^(1/2) - 232*T*corrincphi*muinc^2*muphi^3*sig2inc^(1/2)* sig2phi^(1/2) + 208*corrincphi*muinc*muphi^4*sig2inc^(1/2)* sig2phi^(1/2) + 768*T*corrincphi*muinc^2*muphi^4*sig2inc^(1/2)* sig2phi^(1/2) - 1200*corrincphi*muinc*muphi^5*sig2inc^(1/2)* sig2phi^(1/2) - 1760*T*corrincphi*muinc^2*muphi^5*sig2inc^(1/2)* sig2phi^(1/2) + 4112*corrincphi*muinc*muphi^6*sig2inc^(1/2)* sig2phi^(1/2) + 2912*T*corrincphi*muinc^2*muphi^6*sig2inc^(1/2)* sig2phi^(1/2) - 9376*corrincphi*muinc*muphi^7*sig2inc^(1/2)* sig2phi^(1/2) - 3520*T*corrincphi*muinc^2*muphi^7*sig2inc^(1/2)* sig2phi^(1/2) + 14976*corrincphi*muinc*muphi^8*sig2inc^(1/2)* sig2phi^(1/2) + 3072*T*corrincphi*muinc^2*muphi^8*sig2inc^(1/2)* sig2phi^(1/2) - 17056*corrincphi*muinc*muphi^9*sig2inc^(1/2)* sig2phi^(1/2) - 1856*T*corrincphi*muinc^2*muphi^9*sig2inc^(1/2)* sig2phi^(1/2) + 13728*corrincphi*muinc*muphi^10*sig2inc^(1/2)* sig2phi^(1/2) + 704*T*corrincphi*muinc^2*muphi^10*sig2inc^(1/2)* sig2phi^(1/2) - 7488*corrincphi*muinc*muphi^11*sig2inc^(1/2)* sig2phi^(1/2) - 128*T*corrincphi*muinc^2*muphi^11*sig2inc^(1/2)* sig2phi^(1/2) + 2496*corrincphi*muinc*muphi^12*sig2inc^(1/2)* sig2phi^(1/2) - 384*corrincphi*muinc*muphi^13*sig2inc^(1/2)* sig2phi^(1/2) + 2*T*muinc^3*sig2phi - 28*T*muinc^3*muphi*sig2phi + 12*muinc^2*muphi^2*sig2phi + 172*T*muinc^3*muphi^2*sig2phi - 176*muinc^2*muphi^3*sig2phi - 640*T*muinc^3*muphi^3*sig2phi + 1084*muinc^2*muphi^4*sig2phi + 1616*T*muinc^3*muphi^4*sig2phi - 3832*muinc^2*muphi^5*sig2phi - 2912*T*muinc^3*muphi^5*sig2phi + 8888*muinc^2*muphi^6*sig2phi + 3808*T*muinc^3*muphi^6*sig2phi - 14464*muinc^2*muphi^7*sig2phi - 3584*T*muinc^3*muphi^7*sig2phi + 16936*muinc^2*muphi^8*sig2phi + 2336*T*muinc^3*muphi^8*sig2phi - 14160*muinc^2*muphi^9*sig2phi - 960*T*muinc^3*muphi^9*sig2phi + 8112*muinc^2*muphi^10*sig2phi + 192*T*muinc^3*muphi^10*sig2phi - 2880*muinc^2*muphi^11*sig2phi + 480*muinc^2*muphi^12*sig2phi)/ (T^2*muinc^4*(1 - 2*muphi + 2*muphi^2)^4))/ (((-1 + 2/(1 - E^((delta*muinc)/((-1 + muphi)*muphi))))*delta*(1 - muphi)* muphi)/T + (((4*E^((2*delta*muinc)/((-1 + muphi)*muphi))*delta^3)/ ((1 - E^((delta*muinc)/((-1 + muphi)*muphi)))^3*T*(1 - muphi)* muphi) + (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta^3)/ ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T*(1 - muphi)* muphi))*sig2inc)/2 + corrincphi*((-4*E^((2*delta*muinc)/((-1 + muphi)*muphi))*delta^2* ((delta*muinc)/((1 - muphi)*muphi^2) - (delta*muinc)/((-1 + muphi)^2*muphi)))/ ((1 - E^((delta*muinc)/((-1 + muphi)*muphi)))^3*T) - (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta^2* ((delta*muinc)/((1 - muphi)*muphi^2) - (delta*muinc)/((-1 + muphi)^2*muphi)))/ ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T))*sig2inc^(1/2)* sig2phi^(1/2) + sig2phi + (((-2* (-1 + 2/(1 - E^((delta*muinc)/((-1 + muphi)*muphi))))*delta)/T + (4*E^((delta*muinc)/((-1 + muphi)*muphi))*delta* ((delta*muinc)/((1 - muphi)*muphi^2) - (delta*muinc)/((-1 + muphi)^2*muphi))*(1 - muphi))/ ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T) - (4*E^((delta*muinc)/((-1 + muphi)*muphi))*delta* ((delta*muinc)/((1 - muphi)*muphi^2) - (delta*muinc)/((-1 + muphi)^2*muphi))*muphi)/ ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T) + (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta* ((-2*delta*muinc)/((1 - muphi)*muphi^3) + (2*delta*muinc)/((-1 + muphi)^2*muphi^2) - (2*delta*muinc)/((1 - muphi)^3*muphi))*(1 - muphi)*muphi)/ ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T) + (4*E^((2*delta*muinc)/((-1 + muphi)*muphi))*delta* ((delta*muinc)/((1 - muphi)*muphi^2) - (delta*muinc)/((-1 + muphi)^2*muphi))^2*(1 - muphi)*muphi)/ ((1 - E^((delta*muinc)/((-1 + muphi)*muphi)))^3*T) + (2*E^((delta*muinc)/((-1 + muphi)*muphi))*delta* ((delta*muinc)/((1 - muphi)*muphi^2) - (delta*muinc)/((-1 + muphi)^2*muphi))^2*(1 - muphi)*muphi)/ ((-1 + E^((delta*muinc)/((-1 + muphi)*muphi)))^2*T))*sig2phi)/2) eq5p2Minusinc ::= biaselDivT\$pIninc biaselDivT\$pIninc := phi/delta - inc eq5p2a ::= biasel01DivT\$pIninc biasel01DivT\$pIninc := ((1 - E^((delta*inc)/((-1 + phi)*phi)))* (1 - phi)*phi)/delta - inc eq5p4 ::= En\$p En\$p := (T*((1 - E^((delta*inc)/(-1 + phi)))*(1 - phi) + phi))/delta eq5p5 ::= biasnDivsegsInprev biasnDivsegsInprev := (1 - E^((delta*inc)/(-1 + phi)))*(1 - phi) eq5p5a ::= biasn11DivsegsInprev biasn11DivsegsInprev := -((-1 + E^((delta*inc)/(-1 + phi)))^2* (-1 + phi)) eq5p6 ::= biasnDivsegsIninc biasnDivsegsIninc := -inc + ((1 - E^((delta*inc)/(-1 + phi)))*(1 - phi) + phi)/delta eq5p6a ::= biasn01DivsegsIninc biasn01DivsegsIninc := -inc + ((1 - E^((delta*inc)/(-1 + phi)))* (T/delta - (T*((1 - E^((delta*inc)/(-1 + phi)))*(1 - phi) + phi))/delta))/T (* ---------------------------------------------------------------------- Glossary Adapted for Mathematica Notational Compromises Single Renewal Process d0 Duration in On state for qth occurrence (mean mu1) d1 Duration in Off state for qth occurrence (mean mu0) w = Waiting time between Off-to-On transitions (moments muw,sig2w,mu3w) fw Waiting time probability density function hw Hazard function corresponding to f(w) Behavioral Parameters inc = 1/muw, Incidence phi = mu1/(mu1 + mu0), Prevalence mu1 Event Duration pi Probability an event is of a specified Type, Pr{"success"} Continuous Observation T Length of observation interval for observing individual p N Count of Off-to-On transitions in interval Tp NDivT Empirical rate of Off-to-On transitions in interval Tp vpT Variance-time curve, sampling variance of Np(T)/Tp L Time accumulated in On state in interval Tp LDivT Empirical prevalence in interval Tp Z Accumulated "successes" in interval Tp ZDivN Empirical proportion Multiple Occasions occ; plus subscript indicates a sum over occasions index Npo Counts in interval Tpo on occasion o s = {Tp1 , ..., Tpocc}, Observation schedule (1 x occ vector) s2 1 x occ vector with elements {(Tp1)2 , ..., (Tpocc)2} vpos A 1 x occ vector with elements vpo(T) Recorder Errors Mpr Indicates reported count for individual p "Mangled" by errors from recorder r scriptR Denotes a roster of multiple (size R) recorders Mbar Indicates average of reported counts from a roster of recorders al Recorder allocation matrix over multiple occasions (occ x occ) Representations for Recorder Errors, Heterogeneity Over Occasions kap Thinning parameter for recorder r error process gam Superposition parameter for recorder r error process om Imbalance in recorder r error processes psi Proportional variability from recorder r error processes muom, sig2om, mupsi Moments taken over the population of recorders tau2 variance of Ipo about Ip. Distortion of Incidence on occasion o from heterogeneity over occasions Examples of Renewal Processes Gamma point process: Index, alpha ; Incidence, inc Alternating Poisson: Incidence, inc ; Prevalence, phi (mu1 = phi/inc , mu0 = (1 - phi)/inc) Population of Individual Processes and Psychometric Quantities muinc Mean Incidence over the population of individuals, Ep(Ip) sig2inc Variance of Incidence over population of individuals, Varp(Ip) EpvpT Average sampling variance of Np(T)/Tp over the population of individuals muphi, sig2phi Moments of Prevalence over the population of individuals rhoMM' Inter-recorder correlation; correlation between counts reported by recorders r and r' G Theory quantities for pxoxr Design porVC = {porVCp, porVCo, porVCr, porVCpo, porVCpr, porVCro, porVCpor, porVCeps} Variance components for pxoxr design DEL , del G theory absolute and relative errors porVCabs = porVCo + porVCr + porVCpo + porVCpr + porVCro + (porVCpor + porVCeps) porVCrel = porVCpo + porVCpr + (porVCpor + porVCeps) Definitions of absolute and relative variances Alternatives to Continuous Observation delta Time segment between observations in Instantaneous sampling Time segment for recording in One-zero (interval) sampling segs = T/delta, Number of time segments el Count of 1's for spacing delta and observation interval T in Instantaneous sampling el01 Count of 0-to-1 transitions for spacing delta and observation interval T in Instantaneous sampling record n Count of 1's in One-zero sampling, recording segment delta and observation interval T n01 , n11 Counts of 0-to-1 transitions and of 1-to-1 transitions in the One-zero sampling record, recording segment delta and observation interval T *)