/**************************************************************************************** Date: 6/28/2013 Created By: Brandy Ringham Description: This module calculates the expected value of N(m1), N(m2), and N(m9) (Catellier and Muller, 2000; Ringham et al., in submission). Copyright (C) 2010 Regents of the University of Colorado. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. Inputs: Nt number of independent sampling units p number of repeated measures pi probability missing value k index for missing data summary statistic (see index key) print 1 = print on, otherwise printing is turned off Output: ENmk label for mean of the missing data summary statistic error label for error code matrix Error Code Key: The error matrix contains error indicators for the inputs that are checked in the beginning of the module. Error indicators are defined as follows: 0 no error 1 error The error indicator's position in the error matrix defines the referent input. Row Input 1 error state of Nt 2 error state of p 3 error state of pi 4 error state of k Usage: call NMK( 12, 4, 0.10, 1, 1, ENm1, error ); ****************************************************************************************/ /*begin module definition*/ start NMK( Nt, p, pi, k, print, ENmk, error ); /*initialize error indicators to 0*/ error = j( 4, 1, 0 ); /*check that Nt is valid*/ if Nt < 0 then do; error[ 1, 1 ] = 1; if print = 1 then do; print "invalid number of independent sampling units"; print "Nt should be a positive integer"; print "currently, Nt =" Nt; print "no expected values or variances were calculated"; end; end; /*check that p is valid*/ if p < 0 then do; error[ 2, 1 ] = 1; if print = 1 then do; print "invalid number of repeated measures"; print "p should be a positive integer"; print "currently, p =" p; print "no expected values or variances were calculated"; end; end; /*check that pi is valid*/ if pi < 0 | pi > 1 then do; error[ 3, 1 ] = 1; if print = 1 then do; print "invalid value for pi"; print "pi should be a number between 0 and 1"; print "currently, pi =" pi; print "no expected values or variances were calculated"; end; end; /*check that k is valid*/ if k ^= 1 & k ^= 2 & k ^= 9 then do; error[ 4, 1 ] = 1; if print = 1 then do; print "invalid index for the missing data summary statistic"; print "k should be in the set {1, 2, 9}"; print "currently, k = " k; print "No expected values or variances were calculated"; end; end; /*only complete the next set of steps if there were no errors*/ /*otherwise, program skips over the next block and ends with no calculations*/ if sum( error ) = 0 then do; /*calculate expected value of N(m1)*/ if k = 1 then ENmk = Nt * ( 1 - pi )**p; /*calculate expected value of for N(m2)*/ else if k = 2 then do; if pi = 0 then ENmk = Nt; else if pi = 1 then ENmk = 0; else /*regression model for expected value of N(m2)*/ /*see Table 4, Ringham et al. (in submission)*/ ENmk = 62.7318676 - 5.1567678 * ( Nt / 10 ) - 0.1324363 * ( Nt / 10 )**2 - 1.6042196 * p + 0.0640387 * p**2 - 147.7255861 * ( 1 - pi ) + 87.3472243 * ( 1 - pi )**2 - 0.4981166 * ( Nt / 10 ) * p + 0.0019218 * ( ( Nt / 10 ) * p )**2 + 1.2421550 * p * ( 1 - pi ) - 0.0423721 * ( p * ( 1 - pi ) )**2 + 14.8545994 * ( Nt / 10 ) * ( 1 - pi ) + 0.1812043 * ( ( Nt / 10 ) * ( 1 - pi ) )**2 + 0.4137540 * ( Nt / 10 ) * p * ( 1 - pi ) - 0.0013272 * ( ( Nt / 10 ) * p * ( 1 - pi ) )**2; end; /*calculate expected value of Nm9*/ else if k = 9 then ENmk = Nt * ( 1 - pi ); end; /*end module*/ finish;