/****************************************************************************************
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;