line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
package Statistics::Distributions::Ancova; |
2
|
1
|
|
|
1
|
|
25440
|
use 5.008; |
|
1
|
|
|
|
|
4
|
|
|
1
|
|
|
|
|
42
|
|
3
|
|
|
|
|
|
|
|
4
|
1
|
|
|
1
|
|
6
|
use strict; |
|
1
|
|
|
|
|
2
|
|
|
1
|
|
|
|
|
34
|
|
5
|
1
|
|
|
1
|
|
5
|
use warnings; |
|
1
|
|
|
|
|
6
|
|
|
1
|
|
|
|
|
31
|
|
6
|
1
|
|
|
1
|
|
13
|
use Carp; |
|
1
|
|
|
|
|
2
|
|
|
1
|
|
|
|
|
101
|
|
7
|
1
|
|
|
1
|
|
6
|
use List::Util; |
|
1
|
|
|
|
|
1
|
|
|
1
|
|
|
|
|
84
|
|
8
|
1
|
|
|
1
|
|
930
|
use Math::Cephes qw(:utils); |
|
1
|
|
|
|
|
8391
|
|
|
1
|
|
|
|
|
584
|
|
9
|
1
|
|
|
1
|
|
4898
|
use Contextual::Return; |
|
1
|
|
|
|
|
37331
|
|
|
1
|
|
|
|
|
9
|
|
10
|
1
|
|
|
1
|
|
8933
|
use Perl6::Form; |
|
1
|
|
|
|
|
138264
|
|
|
1
|
|
|
|
|
8
|
|
11
|
1
|
|
|
1
|
|
1178
|
use Statistics::Distributions qw( fprob fdistr); |
|
1
|
|
|
|
|
3620
|
|
|
1
|
|
|
|
|
94
|
|
12
|
|
|
|
|
|
|
|
13
|
|
|
|
|
|
|
=head1 NAME |
14
|
|
|
|
|
|
|
|
15
|
|
|
|
|
|
|
Statistics::Distributions::Ancova - Perl implementation of One-Way Analysis of Covariance for Independent Samples. |
16
|
|
|
|
|
|
|
|
17
|
|
|
|
|
|
|
=cut |
18
|
|
|
|
|
|
|
=head1 VERSION |
19
|
|
|
|
|
|
|
|
20
|
|
|
|
|
|
|
This document describes Statistics::Distributions::Ancova version 0.32.2. |
21
|
|
|
|
|
|
|
|
22
|
|
|
|
|
|
|
=cut |
23
|
1
|
|
|
1
|
|
906
|
use version; our $VERSION = qv('0.32.2'); |
|
1
|
|
|
|
|
2161
|
|
|
1
|
|
|
|
|
6
|
|
24
|
|
|
|
|
|
|
=head1 SYNOPSIS |
25
|
|
|
|
|
|
|
|
26
|
|
|
|
|
|
|
use Statistics::Distributions::Ancova; |
27
|
|
|
|
|
|
|
|
28
|
|
|
|
|
|
|
# Create an Ancova object and set significance value of p = 0.05 for statistical test. See METHODS for optional named arguments and default values. |
29
|
|
|
|
|
|
|
my $anc = Statistics::Distributions::Ancova->new ( { significance => 0.005, input_verbosity => 1, output_verbosity => 1 } ); |
30
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
# Example using k=3 groups. Data includes our dependent variable of interest (Y) and covariant data (X) that is used to eliminate obscuring effects of covariance. |
32
|
|
|
|
|
|
|
my @Drug_A_Y = ('29','27','31','33','32','24','16'); |
33
|
|
|
|
|
|
|
my @Drug_A_X = ('53','64','55','67','55','45','35'); |
34
|
|
|
|
|
|
|
my @Drug_B_Y = ('39','34','20','35','57','28','32','17'); |
35
|
|
|
|
|
|
|
my @Drug_B_X = ('24','19','13','18','25','16','16','13'); |
36
|
|
|
|
|
|
|
my @Drug_C_Y = ('12','21','26','17','25','9','12'); |
37
|
|
|
|
|
|
|
my @Drug_C_X = ('5','12','12','9','12','3','3'); |
38
|
|
|
|
|
|
|
|
39
|
|
|
|
|
|
|
# Data is sent to object as nested HASH reference. Individual group names are option, but to distinguish IV/DV, the names Y and X for the variables are compulsory. |
40
|
|
|
|
|
|
|
my $h_ref = { 'group_A' => { |
41
|
|
|
|
|
|
|
Y => \@Drug_A_Y, |
42
|
|
|
|
|
|
|
X => \@Drug_A_X, |
43
|
|
|
|
|
|
|
}, |
44
|
|
|
|
|
|
|
'group_B' => { |
45
|
|
|
|
|
|
|
Y => \@Drug_B_Y, |
46
|
|
|
|
|
|
|
X => \@Drug_B_X, |
47
|
|
|
|
|
|
|
}, |
48
|
|
|
|
|
|
|
'group_C' => { |
49
|
|
|
|
|
|
|
Y => \@Drug_C_Y, |
50
|
|
|
|
|
|
|
X => \@Drug_C_X, |
51
|
|
|
|
|
|
|
}, |
52
|
|
|
|
|
|
|
}; |
53
|
|
|
|
|
|
|
|
54
|
|
|
|
|
|
|
# Feed the object the data pass data HASH reference with named argument 'data'. |
55
|
|
|
|
|
|
|
$anc->load_data ( { data => $h_ref } ); |
56
|
|
|
|
|
|
|
|
57
|
|
|
|
|
|
|
# Perform analysis |
58
|
|
|
|
|
|
|
$anc->ancova_analysis; |
59
|
|
|
|
|
|
|
|
60
|
|
|
|
|
|
|
# To access results use results method. The return of this method is context dependent (see METHODS). |
61
|
|
|
|
|
|
|
# To print a report to STDOUT call results in VOID context. |
62
|
|
|
|
|
|
|
$anc->results(); |
63
|
|
|
|
|
|
|
|
64
|
|
|
|
|
|
|
=cut |
65
|
|
|
|
|
|
|
=head1 DESCRIPTION |
66
|
|
|
|
|
|
|
|
67
|
|
|
|
|
|
|
ANCOVA is a merger of ANOVA and regression for continuous variables. As with paired t-test and repeated-measures ANOVA |
68
|
|
|
|
|
|
|
this test removes the obscuring effects of pre-existing individual differences among subjects and thus may increase |
69
|
|
|
|
|
|
|
statistical power. In cases where a substantial portion of the variability that occurs within each of the set of a dependent variable Y is |
70
|
|
|
|
|
|
|
actually covariance with another concomitant variable X measures, this test removes the covariance with X from Y thus |
71
|
|
|
|
|
|
|
removing a portion of the irrelevant variability of individual differences. See http://en.wikipedia.org/wiki/Analysis_of_covariance for more info. |
72
|
|
|
|
|
|
|
|
73
|
|
|
|
|
|
|
=cut |
74
|
|
|
|
|
|
|
=head1 Methods |
75
|
|
|
|
|
|
|
|
76
|
|
|
|
|
|
|
=cut |
77
|
|
|
|
|
|
|
|
78
|
|
|
|
|
|
|
####################################################################################################################### |
79
|
|
|
|
|
|
|
|
80
|
|
|
|
|
|
|
sub new { |
81
|
0
|
|
|
0
|
0
|
|
my ($class, $args_h_ref ) = @_; |
82
|
0
|
0
|
0
|
|
|
|
croak qq{\nArguments must be passed as HASH reference.} if ( ( $args_h_ref ) && ( ref $args_h_ref ne q{HASH} ) ); |
83
|
|
|
|
|
|
|
|
84
|
0
|
|
|
|
|
|
my $self = {}; |
85
|
0
|
|
|
|
|
|
bless $self, $class; |
86
|
|
|
|
|
|
|
|
87
|
0
|
|
|
|
|
|
$self->_set_significance($args_h_ref); |
88
|
0
|
|
|
|
|
|
$self->_set_verbosity($args_h_ref); |
89
|
|
|
|
|
|
|
|
90
|
0
|
|
|
|
|
|
return $self; |
91
|
|
|
|
|
|
|
} |
92
|
|
|
|
|
|
|
=head2 new |
93
|
|
|
|
|
|
|
|
94
|
|
|
|
|
|
|
Creates new Statistics::Distributions::Ancova object. Without arguments defaults to a significance test value of p = 0.05. |
95
|
|
|
|
|
|
|
|
96
|
|
|
|
|
|
|
my $anc = Statistics::Distributions::Ancova->new (); |
97
|
|
|
|
|
|
|
|
98
|
|
|
|
|
|
|
Use significance option to set the significance level for the test to values other than 0.05. |
99
|
|
|
|
|
|
|
|
100
|
|
|
|
|
|
|
my $anc = Statistics::Distributions::Ancova->new ( { significance => 0.005 } ); |
101
|
|
|
|
|
|
|
|
102
|
|
|
|
|
|
|
To print data-checking step messages (upon data loading with C) to STDOUT set input_verbosity to 1. |
103
|
|
|
|
|
|
|
|
104
|
|
|
|
|
|
|
my $anc = Statistics::Distributions::Ancova->new ( { input_verbosity => 1 } ); |
105
|
|
|
|
|
|
|
|
106
|
|
|
|
|
|
|
To print a detailed report when C method is called in VOID context to STDOUT set output_verbosity to 1. |
107
|
|
|
|
|
|
|
|
108
|
|
|
|
|
|
|
my $anc = Statistics::Distributions::Ancova->new ( { output_verbosity => 1 } ); |
109
|
|
|
|
|
|
|
|
110
|
|
|
|
|
|
|
=cut |
111
|
|
|
|
|
|
|
|
112
|
|
|
|
|
|
|
#/ now made this private |
113
|
|
|
|
|
|
|
sub _set_significance { |
114
|
0
|
|
|
0
|
|
|
my ($self, $args_h_ref) = @_; |
115
|
|
|
|
|
|
|
|
116
|
0
|
0
|
0
|
|
|
|
croak qq{\nArguments must be passed as HASH reference.} if ( ( $args_h_ref ) && ( ref $args_h_ref ne q{HASH} ) ); |
117
|
0
|
0
|
|
|
|
|
if (!exists $args_h_ref->{significance}) { print qq{\n\nFalling back on default 0.05 significance value.\n} } |
|
0
|
|
|
|
|
|
|
118
|
0
|
0
|
|
|
|
|
my $sig = exists $args_h_ref->{significance} ? $args_h_ref->{significance} : q{0.05}; |
119
|
|
|
|
|
|
|
|
120
|
|
|
|
|
|
|
# included exponential number check |
121
|
0
|
0
|
0
|
|
|
|
croak qq{\nThe p value must be numeric and in the range > 0 and < 1.} if ( $sig !~ /\A \d* \.? \d+ ([eE][+-]?\d+)? \z/xms || $sig <= 0 || $sig >= 1) ; |
|
|
|
0
|
|
|
|
|
122
|
|
|
|
|
|
|
#croak qq{\nThe p value must be numeric and in the range > 0 and < 1.} if ( $sig !~ /\A \d{1,7} \.? \d+ ([eE][+-]?\d+)? \z/xms || $sig <= 0 || $sig >= 1) ; |
123
|
|
|
|
|
|
|
|
124
|
0
|
|
|
|
|
|
$self->{significance} = $sig; |
125
|
0
|
|
|
|
|
|
return; |
126
|
|
|
|
|
|
|
} |
127
|
|
|
|
|
|
|
|
128
|
|
|
|
|
|
|
sub set_significance { |
129
|
|
|
|
|
|
|
|
130
|
0
|
|
|
0
|
0
|
|
my ( $self, $sig ) = @_; |
131
|
|
|
|
|
|
|
#/ in this case we must distinguish between 0 and no arugement! |
132
|
|
|
|
|
|
|
#y the 0 detector |
133
|
0
|
0
|
0
|
|
|
|
croak qq{\nThe p value cannot be 0.} if ( defined $sig && $sig == 0 ) ; |
134
|
|
|
|
|
|
|
#y the no arg detector/empty string - already forced numeric |
135
|
0
|
0
|
|
|
|
|
$sig || print qq{\n\nFalling back on default 0.05 significance value.\n}; |
136
|
0
|
|
0
|
|
|
|
$sig ||= 0.05; |
137
|
|
|
|
|
|
|
#y no need to check for 0 again |
138
|
|
|
|
|
|
|
#/ we don´t need to check the regexp part as we´ve already forced only numeric args with the above ==0... |
139
|
|
|
|
|
|
|
|
140
|
|
|
|
|
|
|
# check for exponentials |
141
|
0
|
0
|
0
|
|
|
|
croak qq{\nThe p value must be numeric and in the range > 0 and < 1.} if ( $sig !~ /\A \d* \.? \d+ ([eE][+-]?\d+)? \z/xms || $sig <= 0 || $sig >= 1) ; |
|
|
|
0
|
|
|
|
|
142
|
|
|
|
|
|
|
#croak qq{\nThe p value must be numeric and in the range > 0 and < 1.} if ( $sig !~ /\A \d{1,7} \.? \d+ ([eE][+-]?\d+)? \z/xms || $sig <= 0 || $sig >= 1) ; |
143
|
|
|
|
|
|
|
#croak qq{\nThe p value must be numeric and in the range > 0 and < 1.} if ( $sig !~ /\A[01]?\.\d+([eE][+-]?\d+)?\z/xms || $sig <= 0 || $sig >= 1) ; |
144
|
|
|
|
|
|
|
|
145
|
0
|
|
|
|
|
|
$self->{significance} = $sig; |
146
|
0
|
|
|
|
|
|
return; |
147
|
|
|
|
|
|
|
|
148
|
|
|
|
|
|
|
} |
149
|
|
|
|
|
|
|
=head2 set_significance |
150
|
|
|
|
|
|
|
|
151
|
|
|
|
|
|
|
Convenience method to reset significance level. Without a value it defaults to p = 0.05 to change this use set_significance. |
152
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
$anc->set_significance(); |
154
|
|
|
|
|
|
|
$anc->set_significance( 0.0005 ); |
155
|
|
|
|
|
|
|
|
156
|
|
|
|
|
|
|
=cut |
157
|
|
|
|
|
|
|
|
158
|
|
|
|
|
|
|
#/ now a private method - only called by new and unload |
159
|
|
|
|
|
|
|
sub _set_verbosity { |
160
|
0
|
|
|
0
|
|
|
my ($self, $args_h_ref) = @_; |
161
|
|
|
|
|
|
|
|
162
|
0
|
0
|
0
|
|
|
|
croak qq{\nArguments must be passed as HASH reference.} if ( ( $args_h_ref ) && ( ref $args_h_ref ne q{HASH} ) ); |
163
|
|
|
|
|
|
|
|
164
|
0
|
0
|
|
|
|
|
my $input_verbosity = exists $args_h_ref->{input_verbosity} ? $args_h_ref->{input_verbosity} : 0 ; |
165
|
0
|
0
|
|
|
|
|
my $output_verbosity = exists $args_h_ref->{output_verbosity} ? $args_h_ref->{output_verbosity} : 0 ; |
166
|
|
|
|
|
|
|
|
167
|
0
|
0
|
|
|
|
|
croak qq{\nInput verbosity must be set to 1 or 0.} if ( $input_verbosity !~ /\A[01]\z/xms ) ; |
168
|
0
|
0
|
|
|
|
|
croak qq{\nOutput verbosity must be set to 1 or 0.} if ( $output_verbosity !~ /\A[01]\z/xms ) ; |
169
|
|
|
|
|
|
|
|
170
|
0
|
|
|
|
|
|
$self->{verbosity} = { input => $input_verbosity, |
171
|
|
|
|
|
|
|
output => $output_verbosity }; |
172
|
|
|
|
|
|
|
|
173
|
|
|
|
|
|
|
#$self->{verbosity} = %{$verbosity}; |
174
|
|
|
|
|
|
|
#$self->{verbosity}{input} = $input_verbosity; |
175
|
|
|
|
|
|
|
#$self->{verbosity}{output} = $output_verbosity; |
176
|
|
|
|
|
|
|
|
177
|
0
|
|
|
|
|
|
return; |
178
|
|
|
|
|
|
|
} |
179
|
|
|
|
|
|
|
|
180
|
|
|
|
|
|
|
sub set_input_verbosity { |
181
|
|
|
|
|
|
|
#/ convinience method to reset output verbosity |
182
|
0
|
|
|
0
|
0
|
|
my ( $self, $verb ) = @_; |
183
|
|
|
|
|
|
|
#y don´t care about distinguishing default no arg and 0 here - unlike in set_significance - so just existince |
184
|
0
|
|
0
|
|
|
|
$verb ||= 0; |
185
|
0
|
0
|
|
|
|
|
croak qq{\nYou must pass set_output_verbosity 1 or 0 (without an arguement it defaults to 0).} if ( $verb !~ /\A[01]\z/xms ) ; |
186
|
0
|
|
|
|
|
|
$self->{verbosity}{input} = $verb; |
187
|
0
|
|
|
|
|
|
return; |
188
|
|
|
|
|
|
|
} |
189
|
|
|
|
|
|
|
=head2 set_input_verbosity |
190
|
|
|
|
|
|
|
|
191
|
|
|
|
|
|
|
Convenience method to reset the input verbosity level. Pass it 1 for verbose and 0 or no argument to leave default |
192
|
|
|
|
|
|
|
silent state. |
193
|
|
|
|
|
|
|
|
194
|
|
|
|
|
|
|
$anc->set_input_verbosity (1); # Turns on verbosity |
195
|
|
|
|
|
|
|
$anc->set_input_verbosity (0); |
196
|
|
|
|
|
|
|
$anc->set_input_verbosity (); |
197
|
|
|
|
|
|
|
|
198
|
|
|
|
|
|
|
=cut |
199
|
|
|
|
|
|
|
|
200
|
|
|
|
|
|
|
sub set_output_verbosity { |
201
|
|
|
|
|
|
|
#/ convinience method to reset output verbosity |
202
|
0
|
|
|
0
|
0
|
|
my ( $self, $verb ) = @_; |
203
|
|
|
|
|
|
|
#y don´t care about distinguishing default no arg and 0 here - unlike in set_significance - so just existince |
204
|
0
|
|
0
|
|
|
|
$verb ||= 0; |
205
|
0
|
0
|
|
|
|
|
croak qq{\nYou must pass set_output_verbosity 1 or 0 (without an arguement it defaults to 0).} if ( $verb !~ /\A[01]\z/xms ) ; |
206
|
0
|
|
|
|
|
|
$self->{verbosity}{output} = $verb; |
207
|
0
|
|
|
|
|
|
return; |
208
|
|
|
|
|
|
|
} |
209
|
|
|
|
|
|
|
=head2 set_output_verbosity |
210
|
|
|
|
|
|
|
|
211
|
|
|
|
|
|
|
Convinience method to reset the output verbosity level. Pass it 1 for verbose and 0 or no argument to leave default |
212
|
|
|
|
|
|
|
silent state. |
213
|
|
|
|
|
|
|
|
214
|
|
|
|
|
|
|
$anc->set_output_verbosity (1); # Turns on verbosity |
215
|
|
|
|
|
|
|
$anc->set_output_verbosity (0); |
216
|
|
|
|
|
|
|
$anc->set_output_verbosity (); |
217
|
|
|
|
|
|
|
|
218
|
|
|
|
|
|
|
=cut |
219
|
|
|
|
|
|
|
|
220
|
|
|
|
|
|
|
sub load_data { |
221
|
0
|
|
|
0
|
0
|
|
my ($self, $h_ref) = @_; |
222
|
|
|
|
|
|
|
|
223
|
0
|
|
|
|
|
|
$self->_pre_check($h_ref); |
224
|
|
|
|
|
|
|
#y unpack the data |
225
|
0
|
|
|
|
|
|
my $data_ref = $h_ref->{data}; |
226
|
|
|
|
|
|
|
|
227
|
0
|
0
|
|
|
|
|
$data_ref or croak qq{\nkey \'data\' points to nothing}; |
228
|
0
|
0
|
|
|
|
|
croak qq{\nThe data pointed to by key \'data\' must be passed as HASH reference.} if ( ref $data_ref ne q{HASH} ); |
229
|
|
|
|
|
|
|
|
230
|
0
|
|
|
|
|
|
$self->_groups_info($data_ref); |
231
|
|
|
|
|
|
|
|
232
|
|
|
|
|
|
|
#/ there is no need to deep copy the data - we just use it... |
233
|
0
|
|
|
|
|
|
$self->_data_check($data_ref); |
234
|
|
|
|
|
|
|
|
235
|
|
|
|
|
|
|
#y construct the array consisting off ALL the data |
236
|
0
|
|
|
|
|
|
$self->_all_array($data_ref); |
237
|
|
|
|
|
|
|
|
238
|
|
|
|
|
|
|
#y set flag |
239
|
0
|
|
|
|
|
|
$self->{analysis_state}{load} = 1; |
240
|
0
|
|
|
|
|
|
return; |
241
|
|
|
|
|
|
|
} |
242
|
|
|
|
|
|
|
=head2 load_data |
243
|
|
|
|
|
|
|
|
244
|
|
|
|
|
|
|
To load or re-load data. Pass the data as named arguement 'data' within an anonymous HASH pointing to nested HASH |
245
|
|
|
|
|
|
|
reference containing the data. Within this HASH reference each subsequent nested HASH corresponds to a separate |
246
|
|
|
|
|
|
|
individual/group. The names of these groups are arbitrary. Within each nested group HASH there must be exactly to keys. |
247
|
|
|
|
|
|
|
One called 'Y' (corresponding to the Dependent Variable that we wish to adjust using covariance) that points to an |
248
|
|
|
|
|
|
|
array ref or directly as an anonymous array of the corresponding data. The other key must be termed 'X' and corresponds |
249
|
|
|
|
|
|
|
to the concomitant variable whose covariation will be used to adjust Y. X is also passed as an array ref/anonymous |
250
|
|
|
|
|
|
|
array. |
251
|
|
|
|
|
|
|
|
252
|
|
|
|
|
|
|
$anc->load_data ( { data => { 'GroupA' => { Y => [qw/ 29 27 31 33 32 24 16 /], X => [qw/ 53 64 55 67 55 45 35 /], }, |
253
|
|
|
|
|
|
|
'GroupB' => { Y => [qw/ 39 34 20 35 57 28 32 17 /], X => [qw/ 24 19 13 18 25 16 16 13 /], }, |
254
|
|
|
|
|
|
|
'GroupC' => { Y => [qw/ 12 21 26 17 25 9 12 /], X => [qw/ 5 12 12 9 12 3 3 /], }, }, |
255
|
|
|
|
|
|
|
} ); |
256
|
|
|
|
|
|
|
|
257
|
|
|
|
|
|
|
=cut |
258
|
|
|
|
|
|
|
|
259
|
|
|
|
|
|
|
sub _pre_check { |
260
|
|
|
|
|
|
|
|
261
|
0
|
|
|
0
|
|
|
my ($self, $h_ref) = @_; |
262
|
|
|
|
|
|
|
|
263
|
|
|
|
|
|
|
#croak qq{\nThe data must be passed as HASH reference.} if ( ( $h_ref ) && ( ref $h_ref ne q{HASH} ) ); |
264
|
|
|
|
|
|
|
#croak qq{\nThe data must be passed as HASH reference pointed to by key \'data\'.} if (!exists $h_ref->{data}); |
265
|
|
|
|
|
|
|
#croak qq{\nYou must pass me some data} if ( !$h_ref ); |
266
|
|
|
|
|
|
|
#croak qq{\nThe data must be passed as HASH reference.} if ( ref $h_ref ne q{HASH} ); |
267
|
|
|
|
|
|
|
#croak qq{\nThe data must be passed as HASH reference pointed to by key \'data\'.} if ( !exists $h_ref->{data}); |
268
|
|
|
|
|
|
|
|
269
|
0
|
0
|
0
|
|
|
|
if ( !$h_ref ) { croak qq{\nYou must pass me some data}; } |
|
0
|
0
|
|
|
|
|
|
270
|
|
|
|
|
|
|
#elsif ( ( ref $h_ref ne q{HASH} ) || (!exists $h_ref->{data}) ) { croak qq{\nThe data must be passed as HASH |
271
|
|
|
|
|
|
|
#elsif ( ref $h_ref ne q{HASH} ) { croak qq{\nThe data must be passed as HASH reference.}; } |
272
|
|
|
|
|
|
|
#elsif ( !exists $h_ref->{data}) { croak qq{\nThe data must be passed as HASH reference pointed to by key \'data\'.}; } |
273
|
0
|
|
|
|
|
|
elsif ( ( ref $h_ref ne q{HASH} ) || ( !exists $h_ref->{data} ) ) { croak qq{\nThe data must be passed within a HASH reference pointed to by key \'data\'.}; } |
274
|
|
|
|
|
|
|
|
275
|
0
|
|
|
|
|
|
return; |
276
|
|
|
|
|
|
|
|
277
|
|
|
|
|
|
|
} |
278
|
|
|
|
|
|
|
|
279
|
|
|
|
|
|
|
sub _groups_info { |
280
|
0
|
|
|
0
|
|
|
my ($self, $h_ref ) = @_; |
281
|
0
|
|
|
|
|
|
my @groups = (keys %{$h_ref}); |
|
0
|
|
|
|
|
|
|
282
|
|
|
|
|
|
|
# exist syntax is exists $hash{blah} thus use ${$h_ref->...}{blah} |
283
|
0
|
0
|
|
|
|
|
croak qq{\n\'T\' is not a permitted name for a group.\n} if ( exists ${$h_ref}{q{T}} ); |
|
0
|
|
|
|
|
|
|
284
|
|
|
|
|
|
|
#y do things this way for safety |
285
|
0
|
|
|
|
|
|
my $k = scalar(@groups); |
286
|
0
|
0
|
0
|
|
|
|
croak qq{\nI need at least 2 groups of data.\n} if ( !$k || $k == 1 ); |
287
|
0
|
|
|
|
|
|
$self->{groups} = [@groups]; |
288
|
0
|
|
|
|
|
|
$self->{k} = $k; |
289
|
0
|
|
|
|
|
|
return; |
290
|
|
|
|
|
|
|
} |
291
|
|
|
|
|
|
|
|
292
|
|
|
|
|
|
|
sub _data_check { |
293
|
0
|
|
|
0
|
|
|
my ($self, $h_ref) = @_; |
294
|
0
|
|
|
|
|
|
my $verbose = $self->{verbosity}{input}; |
295
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
296
|
0
|
|
|
|
|
|
my %group_lengths; |
297
|
0
|
|
|
|
|
|
my $k = $self->{k}; |
298
|
0
|
0
|
|
|
|
|
print qq{\n\nData has k = $k (group number).\n} if $verbose; |
299
|
0
|
|
|
|
|
|
for my $group (@groups) { |
300
|
0
|
0
|
|
|
|
|
croak qq{\n\nEach group must have two sets of data - one for DV and one for IV.\n\n} if ( ( scalar ( keys %{$h_ref->{$group}} ) ) != 2 ); |
|
0
|
|
|
|
|
|
|
301
|
0
|
0
|
|
|
|
|
print qq{\n* Group $group has: }, 0+(keys %{$h_ref->{$group}}), q{ sets of data.} if $verbose; |
|
0
|
|
|
|
|
|
|
302
|
0
|
0
|
|
|
|
|
croak qq{\n\nWe need to distinguish independent and dependent variables so force names of data sets to \x27X/y\047 and \x27Y/y\047.\n} if ( !exists ${$h_ref->{$group}}{q/X/} ); |
|
0
|
|
|
|
|
|
|
303
|
0
|
0
|
|
|
|
|
print qq{\n* Group $group has independent variable X} if $verbose; |
304
|
0
|
0
|
|
|
|
|
croak qq{\n\nWe need to distinguish independent and dependent variables so force names of data sets to \x27X/y\047 and \x27Y/y\047.\n} if ( !exists ${$h_ref->{$group}}{q/Y/} ); |
|
0
|
|
|
|
|
|
|
305
|
0
|
0
|
|
|
|
|
print qq{\n* Group $group has dependent variable Y.} if $verbose; |
306
|
0
|
0
|
0
|
|
|
|
croak qq{\n\nData set must be passed as ARRAY references.\n} if ( ( ref $h_ref->{$group}{q/Y/} ne q{ARRAY} ) || ( ref $h_ref->{$group}{q/X/} ne q{ARRAY} ) ); |
307
|
0
|
0
|
|
|
|
|
print qq{\n* Group $group Y and X are both ARRAY references.} if $verbose; |
308
|
|
|
|
|
|
|
|
309
|
0
|
|
|
|
|
|
my $n_check = scalar(@{$h_ref->{$group}{q/Y/}}); |
|
0
|
|
|
|
|
|
|
310
|
0
|
0
|
0
|
|
|
|
croak qq{\nI need some actual data - sample number is too low.\n} if ( !$n_check || $n_check == 1 ); |
311
|
0
|
0
|
|
|
|
|
print qq{\n* Group $group Y has $n_check data points.} if $verbose; |
312
|
0
|
0
|
|
|
|
|
croak qq{\n\nBoth X and Y data sets must have equal length.\n} if scalar(@{$h_ref->{$group}{q/X/}}) != $n_check; |
|
0
|
|
|
|
|
|
|
313
|
0
|
0
|
|
|
|
|
print qq{\n* Group $group Y also has $n_check data points.} if $verbose; |
314
|
|
|
|
|
|
|
|
315
|
0
|
|
|
|
|
|
$group_lengths{$group} = $n_check; |
316
|
|
|
|
|
|
|
|
317
|
0
|
0
|
|
|
|
|
print qq{\n\nData for group $group looks good.\n} if $verbose; |
318
|
|
|
|
|
|
|
} |
319
|
|
|
|
|
|
|
|
320
|
0
|
0
|
|
|
|
|
print qq{\nData passed. Feeding it to Ancova object.} if $verbose; |
321
|
|
|
|
|
|
|
|
322
|
0
|
|
|
|
|
|
$self->{lengths} = {%group_lengths}; |
323
|
|
|
|
|
|
|
|
324
|
|
|
|
|
|
|
#/ we haven´t deep copied so this is pointless! |
325
|
|
|
|
|
|
|
##s we aren´t actually using that hash passed at all - we are copying it - that way they can use that same hash name again later - i.e. we allocate NEW memory location |
326
|
|
|
|
|
|
|
#y point is data passed checks so we put it into object |
327
|
|
|
|
|
|
|
|
328
|
|
|
|
|
|
|
#y that is while we create a new higher level copy we don´t deep copy so its pointless using this syntax and not |
329
|
|
|
|
|
|
|
#y simply $self->{data} = $h_ref - if we deep copy then we are safe from this issue - clearly T is new data... |
330
|
0
|
|
|
|
|
|
$self->{data} = {%{$h_ref}}; |
|
0
|
|
|
|
|
|
|
331
|
|
|
|
|
|
|
|
332
|
0
|
|
|
|
|
|
return; |
333
|
|
|
|
|
|
|
} |
334
|
|
|
|
|
|
|
|
335
|
|
|
|
|
|
|
sub _all_array { |
336
|
0
|
|
|
0
|
|
|
my ($self, $h_ref ) = @_; |
337
|
|
|
|
|
|
|
#my @groups = (keys %{$h_ref}) == (keys %{$h_ref}) == (keys %{$self->{data}}) |
338
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
339
|
|
|
|
|
|
|
|
340
|
|
|
|
|
|
|
#my $T_list = {}; |
341
|
|
|
|
|
|
|
#@{$T_list->{X}} = (); |
342
|
|
|
|
|
|
|
#@{$T_list->{Y}} = (); |
343
|
0
|
|
|
|
|
|
my $T = {}; |
344
|
|
|
|
|
|
|
|
345
|
0
|
|
|
|
|
|
for my $xy( qw/ X Y / ) { |
346
|
|
|
|
|
|
|
|
347
|
0
|
|
|
|
|
|
for my $group (@groups) { |
348
|
|
|
|
|
|
|
|
349
|
|
|
|
|
|
|
#y needs pre-initialisation of everything! |
350
|
|
|
|
|
|
|
#@{$T_list->{$xy}} = (@{$T_list->{$xy}}, @{$h_ref->{$group}{$xy}}); |
351
|
0
|
|
|
|
|
|
push @{$T->{$xy}}, @{$h_ref->{$group}{$xy}}; |
|
0
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
352
|
|
|
|
|
|
|
|
353
|
|
|
|
|
|
|
} |
354
|
|
|
|
|
|
|
} |
355
|
|
|
|
|
|
|
|
356
|
0
|
|
|
|
|
|
$self->{data}{T} = {%{$T}}; |
|
0
|
|
|
|
|
|
|
357
|
0
|
|
|
|
|
|
return; |
358
|
|
|
|
|
|
|
} |
359
|
|
|
|
|
|
|
|
360
|
|
|
|
|
|
|
sub print_data { |
361
|
0
|
|
|
0
|
0
|
|
my $self = shift; |
362
|
0
|
0
|
|
|
|
|
croak qq{\nYou have to load some data first.} if !defined ${$self}{groups}; |
|
0
|
|
|
|
|
|
|
363
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
364
|
0
|
|
|
|
|
|
for my $group (@groups) { |
365
|
0
|
|
|
|
|
|
for my $xy ( qw / X Y / ) { |
366
|
0
|
|
|
|
|
|
my @array = @{$self->{data}{$group}{$xy}}; |
|
0
|
|
|
|
|
|
|
367
|
0
|
|
|
|
|
|
print qq{\n\nGroup $group - data set $xy\n@array.}; |
368
|
|
|
|
|
|
|
} |
369
|
|
|
|
|
|
|
} |
370
|
0
|
|
|
|
|
|
return; |
371
|
|
|
|
|
|
|
} |
372
|
|
|
|
|
|
|
|
373
|
|
|
|
|
|
|
sub unload { |
374
|
0
|
|
|
0
|
0
|
|
my $self = shift; |
375
|
0
|
0
|
|
|
|
|
croak qq{\nYou have to load some data before calling this method} if ( !exists $self->{analysis_state}{load} ); |
376
|
0
|
|
|
|
|
|
my @object_keys = keys %{$self}; |
|
0
|
|
|
|
|
|
|
377
|
|
|
|
|
|
|
OBJECT: |
378
|
0
|
|
|
|
|
|
foreach (@object_keys) { |
379
|
0
|
0
|
|
|
|
|
next OBJECT if $_ eq q{data}; |
380
|
0
|
|
|
|
|
|
$self->{$_} = undef; |
381
|
|
|
|
|
|
|
} |
382
|
0
|
|
|
|
|
|
$self->{data} = {}; # empty h_ref - thus wipe out old data. |
383
|
0
|
|
|
|
|
|
$self->{significance} = 0.05; |
384
|
0
|
|
|
|
|
|
$self->_set_verbosity; |
385
|
0
|
|
|
|
|
|
return; |
386
|
|
|
|
|
|
|
} |
387
|
|
|
|
|
|
|
=head2 unload |
388
|
|
|
|
|
|
|
|
389
|
|
|
|
|
|
|
To clear the object use unload. |
390
|
|
|
|
|
|
|
|
391
|
|
|
|
|
|
|
$anc->unload; |
392
|
|
|
|
|
|
|
|
393
|
|
|
|
|
|
|
=cut |
394
|
|
|
|
|
|
|
|
395
|
|
|
|
|
|
|
sub load_data_old { |
396
|
0
|
|
|
0
|
0
|
|
my ($self, $h_ref ) = @_; |
397
|
0
|
|
|
|
|
|
my $T_y_ref = [ (@{$h_ref->{A}{Y}}, @{$h_ref->{B}{Y}}) ]; |
|
0
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
398
|
0
|
|
|
|
|
|
my $T_x_ref = [ (@{$h_ref->{A}{X}}, @{$h_ref->{B}{X}}) ]; |
|
0
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
399
|
0
|
|
|
|
|
|
$self->{data} = $h_ref; |
400
|
0
|
|
|
|
|
|
$self->{data}{T} = { X => $T_x_ref, |
401
|
|
|
|
|
|
|
Y => $T_y_ref, |
402
|
|
|
|
|
|
|
}; |
403
|
0
|
|
|
|
|
|
return; |
404
|
|
|
|
|
|
|
} |
405
|
|
|
|
|
|
|
|
406
|
|
|
|
|
|
|
sub ancova_analysis { |
407
|
0
|
|
|
0
|
0
|
|
my $self = shift; |
408
|
0
|
0
|
|
|
|
|
croak qq{\nYou have to load some data before calling this method} if ( !exists $self->{analysis_state}{load} ); |
409
|
0
|
|
|
|
|
|
$self->all_SS; |
410
|
0
|
|
|
|
|
|
$self->all_SC; |
411
|
0
|
|
|
|
|
|
$self->adjustments_for_correlation; |
412
|
|
|
|
|
|
|
|
413
|
|
|
|
|
|
|
#y set flag |
414
|
0
|
|
|
|
|
|
$self->{analysis_state}{analysis} = 1; |
415
|
0
|
|
|
|
|
|
return; |
416
|
|
|
|
|
|
|
} |
417
|
|
|
|
|
|
|
=head2 ancova_analysis |
418
|
|
|
|
|
|
|
|
419
|
|
|
|
|
|
|
To perform the analysis. |
420
|
|
|
|
|
|
|
|
421
|
|
|
|
|
|
|
$anc->ancova_analysis; |
422
|
|
|
|
|
|
|
|
423
|
|
|
|
|
|
|
=cut |
424
|
|
|
|
|
|
|
|
425
|
|
|
|
|
|
|
sub all_SS { |
426
|
0
|
|
|
0
|
0
|
|
my $self = shift; |
427
|
0
|
|
|
|
|
|
for my $xy ( qw / X Y / ) { |
428
|
0
|
|
|
|
|
|
$self->group_SS($xy); |
429
|
0
|
|
|
|
|
|
$self->SS_variants($xy); |
430
|
|
|
|
|
|
|
} |
431
|
0
|
|
|
|
|
|
return; |
432
|
|
|
|
|
|
|
} |
433
|
|
|
|
|
|
|
|
434
|
|
|
|
|
|
|
sub group_SS { |
435
|
0
|
|
|
0
|
0
|
|
my ($self, $xy) = @_; |
436
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
437
|
0
|
|
|
|
|
|
for my $group ( @groups, qq{T} ) { |
438
|
0
|
|
|
|
|
|
$self->SS( $group, $xy ); |
439
|
|
|
|
|
|
|
} |
440
|
0
|
|
|
|
|
|
return; |
441
|
|
|
|
|
|
|
} |
442
|
|
|
|
|
|
|
|
443
|
|
|
|
|
|
|
sub SS { |
444
|
0
|
|
|
0
|
0
|
|
my ($self, $subject, $variable) = @_; |
445
|
0
|
|
|
|
|
|
my $a_ref = $self->{data}{$subject}{$variable}; |
446
|
|
|
|
|
|
|
#my $n = @{$a_ref}; |
447
|
0
|
|
|
|
|
|
my $n = scalar(@{$a_ref}); |
|
0
|
|
|
|
|
|
|
448
|
0
|
|
|
|
|
|
my $sum = List::Util::sum @{$a_ref}; |
|
0
|
|
|
|
|
|
|
449
|
|
|
|
|
|
|
#my $mean = ( $sum / @{$a_ref} ); |
450
|
0
|
|
|
|
|
|
my $mean = ( $sum / scalar(@{$a_ref}) ); |
|
0
|
|
|
|
|
|
|
451
|
0
|
|
|
|
|
|
my $square_of_sum = Math::Cephes::pow ( $sum, 2 ); |
452
|
|
|
|
|
|
|
##e SS = ( sum ( Xi**2 ) ) - ( sum ( Xi ) )**2 / n |
453
|
0
|
|
|
|
|
|
my $sum_of_squares = List::Util::sum map { Math::Cephes::pow ( $_, 2 ) } @{$a_ref}; |
|
0
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
454
|
0
|
|
|
|
|
|
my $SS = $sum_of_squares - ( $square_of_sum / @{$a_ref} ); |
|
0
|
|
|
|
|
|
|
455
|
|
|
|
|
|
|
#my $_SS = $_sum_of_squares - ( $_square_of_sum / $#{$a_ref}+1 ); |
456
|
|
|
|
|
|
|
|
457
|
|
|
|
|
|
|
# feed the object |
458
|
0
|
|
|
|
|
|
$self->{SS}{$subject}{$variable}{sum} = $sum; |
459
|
0
|
|
|
|
|
|
$self->{SS}{$subject}{$variable}{mean} = $mean; |
460
|
0
|
|
|
|
|
|
$self->{SS}{$subject}{$variable}{square_of_sum} = $square_of_sum; |
461
|
0
|
|
|
|
|
|
$self->{SS}{$subject}{$variable}{sum_of_squares} = $sum_of_squares; |
462
|
0
|
|
|
|
|
|
$self->{SS}{$subject}{$variable}{SS} = $SS; |
463
|
0
|
|
|
|
|
|
$self->{SS}{$subject}{$variable}{n} = $n; |
464
|
|
|
|
|
|
|
# return $sum, $square_of_sum, $sum_of_squares, $SS; |
465
|
0
|
|
|
|
|
|
return; |
466
|
|
|
|
|
|
|
} |
467
|
|
|
|
|
|
|
|
468
|
|
|
|
|
|
|
sub SS_variants { |
469
|
0
|
|
|
0
|
0
|
|
my ($self, $variable) = @_; |
470
|
|
|
|
|
|
|
|
471
|
|
|
|
|
|
|
##o the calculation involves getting SS_total - i.e. just the SS applied to ALL values of Y or X. then requires SS_within_group - just the sum of each groups SS and then |
472
|
|
|
|
|
|
|
##o finally the SS_back_ground (though not for X - just Y) - this is just the SS_toatl - SS_within_group |
473
|
|
|
|
|
|
|
|
474
|
|
|
|
|
|
|
##s pull the SS_total - i.e. SS method was called on T array containing all X or Y entries |
475
|
0
|
|
|
|
|
|
my $SS_Total = $self->{SS}{T}{$variable}{SS}; |
476
|
|
|
|
|
|
|
#y call wg method - needs to know which variable we´re using |
477
|
0
|
|
|
|
|
|
my $SS_wg = $self->_SS_wg($variable); |
478
|
0
|
|
|
|
|
|
$self->{SS}{$variable}{between_group} = ( $SS_Total - $SS_wg ); # parenthesis are just to make it easier to see what´s happening |
479
|
0
|
|
|
|
|
|
$self->{SS}{$variable}{within_group} = $SS_wg; |
480
|
0
|
|
|
|
|
|
$self->{SS}{$variable}{total} = $SS_Total; |
481
|
0
|
|
|
|
|
|
return; |
482
|
|
|
|
|
|
|
} |
483
|
|
|
|
|
|
|
|
484
|
|
|
|
|
|
|
sub _SS_wg { |
485
|
0
|
|
|
0
|
|
|
my ($self, $variable) = @_; |
486
|
|
|
|
|
|
|
##s we need to sum the group SS scores for SS_within_group |
487
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
488
|
0
|
|
|
|
|
|
my $SS_wg = 0; |
489
|
0
|
|
|
|
|
|
for my $group (@groups) { |
490
|
|
|
|
|
|
|
|
491
|
0
|
|
|
|
|
|
my $SS_group = $self->{SS}{$group}{$variable}{SS}; |
492
|
0
|
|
|
|
|
|
$SS_wg += $SS_group; |
493
|
|
|
|
|
|
|
|
494
|
|
|
|
|
|
|
} |
495
|
0
|
|
|
|
|
|
return $SS_wg; |
496
|
|
|
|
|
|
|
} |
497
|
|
|
|
|
|
|
|
498
|
|
|
|
|
|
|
sub all_SC { |
499
|
0
|
|
|
0
|
0
|
|
my ($self, $xy) = @_; |
500
|
0
|
|
|
|
|
|
$self->group_SC; |
501
|
0
|
|
|
|
|
|
$self->{SC}{T}{sum_of_X_and_Y_SC_within_group} = $self->_SC_wg; |
502
|
0
|
|
|
|
|
|
return; |
503
|
|
|
|
|
|
|
} |
504
|
|
|
|
|
|
|
|
505
|
|
|
|
|
|
|
sub group_SC { |
506
|
0
|
|
|
0
|
0
|
|
my $self = shift; |
507
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
508
|
|
|
|
|
|
|
#y loop through all the groups and the Total array |
509
|
0
|
|
|
|
|
|
for my $group ( @groups, qq{T} ) { |
510
|
0
|
|
|
|
|
|
$self->SC ( $group ); |
511
|
|
|
|
|
|
|
} |
512
|
0
|
|
|
|
|
|
return; |
513
|
|
|
|
|
|
|
} |
514
|
|
|
|
|
|
|
|
515
|
|
|
|
|
|
|
sub _SC_wg { |
516
|
0
|
|
|
0
|
|
|
my ($self) = @_; |
517
|
|
|
|
|
|
|
##s we need to sum the group SC scores for SS_within_group |
518
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
519
|
0
|
|
|
|
|
|
my $SC_wg = 0; |
520
|
|
|
|
|
|
|
|
521
|
0
|
|
|
|
|
|
for my $group (@groups) { |
522
|
0
|
|
|
|
|
|
my $SC_group = $self->{SC}{$group}{SC_within_group}; |
523
|
0
|
|
|
|
|
|
$SC_wg += $SC_group; |
524
|
|
|
|
|
|
|
} |
525
|
0
|
|
|
|
|
|
return $SC_wg; |
526
|
|
|
|
|
|
|
} |
527
|
|
|
|
|
|
|
|
528
|
|
|
|
|
|
|
sub SC { |
529
|
|
|
|
|
|
|
##o just calculate covariates - i.e. X * Y in place of X**2... |
530
|
0
|
|
|
0
|
0
|
|
my ( $self, $subject ) = @_; |
531
|
0
|
|
|
|
|
|
my $subject_y = $self->{data}{$subject}{Y}; |
532
|
0
|
|
|
|
|
|
my $subject_x = $self->{data}{$subject}{X}; |
533
|
0
|
|
|
|
|
|
my $product_xy_sum; |
534
|
|
|
|
|
|
|
|
535
|
0
|
|
|
|
|
|
for (0..$#{$subject_x}) { |
|
0
|
|
|
|
|
|
|
536
|
0
|
|
|
|
|
|
my $val = $subject_x->[$_] * $subject_y->[$_]; |
537
|
0
|
|
|
|
|
|
$product_xy_sum += $val; |
538
|
|
|
|
|
|
|
} |
539
|
|
|
|
|
|
|
|
540
|
0
|
|
|
|
|
|
my $subject_x_sum = List::Util::sum @{$subject_x}; |
|
0
|
|
|
|
|
|
|
541
|
0
|
|
|
|
|
|
my $subject_y_sum = List::Util::sum @{$subject_y}; |
|
0
|
|
|
|
|
|
|
542
|
|
|
|
|
|
|
|
543
|
0
|
|
|
|
|
|
my $SC = $product_xy_sum - ( $subject_x_sum * $subject_y_sum ) / @{$subject_x}; |
|
0
|
|
|
|
|
|
|
544
|
|
|
|
|
|
|
|
545
|
|
|
|
|
|
|
# feed object - probably ought to always use this syntax to feed multiple values |
546
|
0
|
|
|
|
|
|
$self->{SC}{$subject} = { sum_of_xy_products => $product_xy_sum, |
547
|
|
|
|
|
|
|
sum_of_x => $subject_x_sum, |
548
|
|
|
|
|
|
|
sum_of_y => $subject_y_sum, |
549
|
|
|
|
|
|
|
SC_within_group => $SC |
550
|
|
|
|
|
|
|
}; |
551
|
0
|
|
|
|
|
|
return; |
552
|
|
|
|
|
|
|
} |
553
|
|
|
|
|
|
|
|
554
|
|
|
|
|
|
|
sub adjustments_for_correlation { |
555
|
|
|
|
|
|
|
##o this runs all the adjustment methods as nested private methods |
556
|
0
|
|
|
0
|
0
|
|
my $self = shift; |
557
|
0
|
|
|
|
|
|
$self->_adjust_SS_Y_total; |
558
|
0
|
|
|
|
|
|
$self->_adjust_SS_Y_wg; |
559
|
0
|
|
|
|
|
|
$self->_adjust_SS_Y_bg; |
560
|
0
|
|
|
|
|
|
$self->_adjust_Y_means; |
561
|
0
|
|
|
|
|
|
$self->_analysis_covariance_with_adjusted_SS; |
562
|
0
|
|
|
|
|
|
return; |
563
|
|
|
|
|
|
|
} |
564
|
|
|
|
|
|
|
|
565
|
|
|
|
|
|
|
sub _adjust_SS_Y_total { |
566
|
|
|
|
|
|
|
##o Adjusting SS_Y_total in light of covariance with X - 4a |
567
|
0
|
|
|
0
|
|
|
my $self = shift; |
568
|
|
|
|
|
|
|
|
569
|
|
|
|
|
|
|
##e r_T = SC_T (this is the within group measure for all data = SC_T_SC_within_group) / ( sqrt ( SS_(X)_Total * SS_Y_Total ) |
570
|
|
|
|
|
|
|
|
571
|
0
|
|
|
|
|
|
my $SS_Y_total = $self->{SS}{Y}{total}; |
572
|
0
|
|
|
|
|
|
my $r_T = $self->{SC}{T}{SC_within_group} / sqrt ( $self->{SS}{X}{total} * $SS_Y_total ); |
573
|
|
|
|
|
|
|
|
574
|
|
|
|
|
|
|
##e The proportion of the total variability of Y attributable to its covariance with X / r_T_sq = r_T**2 |
575
|
|
|
|
|
|
|
|
576
|
0
|
|
|
|
|
|
my $r_T_sq = Math::Cephes::pow ( $r_T, 2); |
577
|
|
|
|
|
|
|
|
578
|
|
|
|
|
|
|
##s we adjust SS_Y_Total by removing from it this proportion of covariance. (1) we get this proportion of covariance proportion_of_SS_Y_Total = SS_Y_Total * r_T_sq. |
579
|
|
|
|
|
|
|
##s (2) we subtract that from SS_Y_Total to get SS_Y_Total_Adj |
580
|
|
|
|
|
|
|
|
581
|
0
|
|
|
|
|
|
my $SS_Y_total_Adj = $SS_Y_total - ( $SS_Y_total * $r_T_sq ); |
582
|
|
|
|
|
|
|
|
583
|
|
|
|
|
|
|
##e - this is an algerbraic equivalent to prevent excessive rounding of r__T_sq: SS_Y_Total_Adj = SS_Y_Total - ( SC_Total / SS_X_Total ) |
584
|
|
|
|
|
|
|
|
585
|
|
|
|
|
|
|
# send to object |
586
|
0
|
|
|
|
|
|
$self->{SS}{Y}{total_adjusted} = $SS_Y_total_Adj; |
587
|
0
|
|
|
|
|
|
$self->{output}{r_T} = $r_T; |
588
|
0
|
|
|
|
|
|
$self->{output}{r_T_sq} = $r_T_sq; |
589
|
|
|
|
|
|
|
|
590
|
0
|
|
|
|
|
|
return; |
591
|
|
|
|
|
|
|
} |
592
|
|
|
|
|
|
|
|
593
|
|
|
|
|
|
|
sub _adjust_SS_Y_wg { |
594
|
|
|
|
|
|
|
##o Adjusting SS_Y_wg on basis of covariance - 4b |
595
|
0
|
|
|
0
|
|
|
my $self = shift; |
596
|
|
|
|
|
|
|
|
597
|
|
|
|
|
|
|
##e r_wg = SC_Total_wg / sqrt (SS_X_wg * SS_Y_wg ) |
598
|
|
|
|
|
|
|
|
599
|
0
|
|
|
|
|
|
my $SS_Y_wg = $self->{SS}{Y}{within_group}; |
600
|
0
|
|
|
|
|
|
my $r_wg = $self->{SC}{T}{sum_of_X_and_Y_SC_within_group} / sqrt ($self->{SS}{X}{within_group} * $SS_Y_wg ); |
601
|
|
|
|
|
|
|
|
602
|
|
|
|
|
|
|
##e The proportion of the within-groups variability of Y attributable to covariance with X / r_wg_sq = r_wg**2 |
603
|
|
|
|
|
|
|
|
604
|
0
|
|
|
|
|
|
my $r_wg_sq = Math::Cephes::pow ( $r_wg, 2); |
605
|
|
|
|
|
|
|
|
606
|
0
|
|
|
|
|
|
my $SS_Y_wg_Adj = $SS_Y_wg - ( $SS_Y_wg * $r_wg_sq ); |
607
|
|
|
|
|
|
|
|
608
|
|
|
|
|
|
|
# send to object |
609
|
0
|
|
|
|
|
|
$self->{SS}{Y}{within_group_adjusted} = $SS_Y_wg_Adj; |
610
|
0
|
|
|
|
|
|
$self->{output}{r_wg} = $r_wg; |
611
|
0
|
|
|
|
|
|
$self->{output}{r_wg_sq} = $r_wg_sq; |
612
|
|
|
|
|
|
|
|
613
|
0
|
|
|
|
|
|
return; |
614
|
|
|
|
|
|
|
|
615
|
|
|
|
|
|
|
} |
616
|
|
|
|
|
|
|
|
617
|
|
|
|
|
|
|
sub _adjust_SS_Y_bg { |
618
|
|
|
|
|
|
|
##o Adjustment of SS_Y_bg - 4c |
619
|
0
|
|
|
0
|
|
|
my $self = shift; |
620
|
|
|
|
|
|
|
|
621
|
|
|
|
|
|
|
##e SS_Y_bg_Adj = SS_Y_Total_Adj — SS_Y_wg_Adj |
622
|
|
|
|
|
|
|
|
623
|
0
|
|
|
|
|
|
my $SS_Y_bg_Adj = $self->{SS}{Y}{total_adjusted} - $self->{SS}{Y}{within_group_adjusted}; |
624
|
|
|
|
|
|
|
|
625
|
|
|
|
|
|
|
# send to object |
626
|
0
|
|
|
|
|
|
$self->{SS}{Y}{between_group_adjusted} = $SS_Y_bg_Adj; |
627
|
0
|
|
|
|
|
|
return; |
628
|
|
|
|
|
|
|
} |
629
|
|
|
|
|
|
|
|
630
|
|
|
|
|
|
|
sub _adjust_Y_means { |
631
|
|
|
|
|
|
|
##o Adjustment of the Means of Y for Groups A and B - 4d |
632
|
0
|
|
|
0
|
|
|
my $self = shift; |
633
|
|
|
|
|
|
|
|
634
|
|
|
|
|
|
|
##e bwg / slope_aggreage_wg = SC_wg / SS_X_wg |
635
|
|
|
|
|
|
|
|
636
|
0
|
|
|
|
|
|
my $slope_aggregate_wg = $self->{SC}{T}{sum_of_X_and_Y_SC_within_group} / $self->{SS}{X}{within_group}; |
637
|
|
|
|
|
|
|
|
638
|
0
|
|
|
|
|
|
$self->_adjust_each_mean($slope_aggregate_wg); |
639
|
0
|
|
|
|
|
|
return; |
640
|
|
|
|
|
|
|
} |
641
|
|
|
|
|
|
|
|
642
|
|
|
|
|
|
|
sub _adjust_each_mean { |
643
|
0
|
|
|
0
|
|
|
my ($self, $slope_aggregate_wg) = @_; |
644
|
|
|
|
|
|
|
|
645
|
|
|
|
|
|
|
##s $self->{SS}{T}{X}{mean} will be used for each group |
646
|
0
|
|
|
|
|
|
my $Mean_X_for_all_samples = $self->{SS}{T}{X}{mean}; |
647
|
|
|
|
|
|
|
|
648
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
649
|
|
|
|
|
|
|
|
650
|
0
|
|
|
|
|
|
for my $group (@groups) { |
651
|
|
|
|
|
|
|
|
652
|
|
|
|
|
|
|
##e Mean_Y_for_A_Adj = Mean_Y_for_A — slope_aggregate_wg (Mean_X_for_A — Mean_X_for_Total - i.e. all samples) |
653
|
0
|
|
|
|
|
|
my $Mean_group_Y_Adj = $self->{SS}{$group}{Y}{mean} - ( $slope_aggregate_wg * ( $self->{SS}{$group}{X}{mean} - $Mean_X_for_all_samples ) ); |
654
|
|
|
|
|
|
|
|
655
|
|
|
|
|
|
|
#s send to object |
656
|
0
|
|
|
|
|
|
$self->{SS}{$group}{Y}{mean_adjusted} = $Mean_group_Y_Adj; |
657
|
|
|
|
|
|
|
|
658
|
|
|
|
|
|
|
} |
659
|
0
|
|
|
|
|
|
return; |
660
|
|
|
|
|
|
|
} |
661
|
|
|
|
|
|
|
|
662
|
|
|
|
|
|
|
sub _analysis_covariance_with_adjusted_SS { |
663
|
|
|
|
|
|
|
##o Analysis of Covariance Using Adjusted Values of SS - calculating F - 4e |
664
|
0
|
|
|
0
|
|
|
my $self = shift; |
665
|
0
|
|
|
|
|
|
my $k = $self->{k}; |
666
|
|
|
|
|
|
|
|
667
|
|
|
|
|
|
|
##o In ANOVA the within-group variance df is: Nt (total number of subjects — k (number of groups). |
668
|
|
|
|
|
|
|
##o In ANCOVA the within-groups df is reduced by 1 due accomodate the fact that the CV portion of within-groups variability has been removed from the analysis. |
669
|
|
|
|
|
|
|
|
670
|
|
|
|
|
|
|
##e df_Y_wg_Adj = df_Y_wg - 1; = NT (total number of measurements) — k (total number of groups/individuals/things) — 1 - here = 20 - 2 - 1) = 17 |
671
|
|
|
|
|
|
|
|
672
|
0
|
|
|
|
|
|
my $df_wg_Y_Adj = $self->{SS}{T}{X}{n} - $k - 1; |
673
|
|
|
|
|
|
|
|
674
|
|
|
|
|
|
|
##o The df for between-groups remains the same as for one-way ANOVA |
675
|
|
|
|
|
|
|
|
676
|
|
|
|
|
|
|
##e df_Y_bg = k — 1 - here = 2 — 1 = 1 |
677
|
|
|
|
|
|
|
|
678
|
0
|
|
|
|
|
|
my $df_bg_Y = $k - 1; |
679
|
|
|
|
|
|
|
|
680
|
|
|
|
|
|
|
##e we use F = ( SS_bg_Y_Adj / df_bg_Y ) / ( SS_wg_Y_Adj / df_wg_Y_Adj ) |
681
|
|
|
|
|
|
|
|
682
|
|
|
|
|
|
|
##e MS_bg is SS_bg_Y_Adj / df_bg_Y and |
683
|
|
|
|
|
|
|
|
684
|
0
|
|
|
|
|
|
my $MS_bg = ( $self->{SS}{Y}{between_group_adjusted} / $df_bg_Y ); |
685
|
|
|
|
|
|
|
|
686
|
|
|
|
|
|
|
##e MS_wg is SS_wg_Y_Adj / df_wg_Y_Adj |
687
|
|
|
|
|
|
|
|
688
|
0
|
|
|
|
|
|
my $MS_wg = ( $self->{SS}{Y}{within_group_adjusted} / $df_wg_Y_Adj ); |
689
|
|
|
|
|
|
|
|
690
|
|
|
|
|
|
|
##e F is usually expressed as MS_bg / MS_wg - |
691
|
|
|
|
|
|
|
|
692
|
0
|
|
|
|
|
|
my $F = ( $MS_bg / $MS_wg ); |
693
|
|
|
|
|
|
|
|
694
|
|
|
|
|
|
|
# feed to $self |
695
|
0
|
|
|
|
|
|
$self->{output}{df_Y_wg_Adj} = $df_wg_Y_Adj; |
696
|
0
|
|
|
|
|
|
$self->{output}{df_Y_bg} = $df_bg_Y; |
697
|
0
|
|
|
|
|
|
$self->{output}{MS_bg} = $MS_bg; |
698
|
0
|
|
|
|
|
|
$self->{output}{MS_wg} = $MS_wg; |
699
|
0
|
|
|
|
|
|
$self->{output}{F_score} = $F; |
700
|
|
|
|
|
|
|
|
701
|
|
|
|
|
|
|
# $self->{output} = { df_Y_wg_Adj => $df_wg_Y_Adj, |
702
|
|
|
|
|
|
|
# df_Y_bg => $df_bg_Y, |
703
|
|
|
|
|
|
|
# MS_bg => $MS_bg, |
704
|
|
|
|
|
|
|
# MS_wg => $MS_wg, |
705
|
|
|
|
|
|
|
# F_score => $F, |
706
|
|
|
|
|
|
|
# }; |
707
|
0
|
|
|
|
|
|
return; |
708
|
|
|
|
|
|
|
|
709
|
|
|
|
|
|
|
} |
710
|
|
|
|
|
|
|
|
711
|
|
|
|
|
|
|
sub results { |
712
|
|
|
|
|
|
|
|
713
|
|
|
|
|
|
|
# unpack rest of @_ - may go to verbose or list printing |
714
|
0
|
|
|
0
|
0
|
|
my @other_args = @_; |
715
|
0
|
|
|
|
|
|
my $self = shift @other_args; |
716
|
|
|
|
|
|
|
|
717
|
0
|
0
|
|
|
|
|
croak qq{\nYou have to load some data before calling this method} if ( !exists $self->{analysis_state}{load} ); |
718
|
0
|
0
|
|
|
|
|
croak qq{\nYou have to run ancova_analysis before calling this method} if ( !exists $self->{analysis_state}{analysis} ); |
719
|
|
|
|
|
|
|
##o get standard F values and generate messages |
720
|
|
|
|
|
|
|
|
721
|
|
|
|
|
|
|
#my ( $self, $verbose ) = @_; |
722
|
|
|
|
|
|
|
#$verbose ||= 0; |
723
|
|
|
|
|
|
|
#my $self = shift; |
724
|
|
|
|
|
|
|
#my $verbose = shift ||= 0; |
725
|
|
|
|
|
|
|
#$verbose = $verbose eq q{verbose} ? 1 : 0 ; |
726
|
|
|
|
|
|
|
|
727
|
0
|
|
|
|
|
|
my $df_wg_Y_Adj = $self->{output}{df_Y_wg_Adj}; |
728
|
0
|
|
|
|
|
|
my $df_bg_Y = $self->{output}{df_Y_bg}; |
729
|
0
|
|
|
|
|
|
my $F = $self->{output}{F_score}; |
730
|
|
|
|
|
|
|
|
731
|
|
|
|
|
|
|
#@{$self->{output}{standard_F_values}} = map { my $standard_F = fdistr ( $df_bg_Y, $df_wg_Y_Adj, $_ ) ; { standard_F => $standard_F, p_val => $_ } } |
732
|
|
|
|
|
|
|
# (0.005, 0.01, 0.05, 0.1); # using standard values of p |
733
|
|
|
|
|
|
|
|
734
|
|
|
|
|
|
|
#if ( $F > Statistics::Distributions::fdistr ($df_bg_Y,$df_wg_Y_Adj,0.01) ) { print qq{\n\nthis value of F is significant at the p=0.01 level} } |
735
|
|
|
|
|
|
|
#elsif ( $F > Statistics::Distributions::fdistr ($df_bg_Y,$df_wg_Y_Adj,0.05) ) { print qq{\n\nthis value of F is significant at the p=0.05 level} } |
736
|
|
|
|
|
|
|
#else { print qq{\n\nthis value of F is not significant } } |
737
|
|
|
|
|
|
|
|
738
|
0
|
|
|
|
|
|
my $chosen_p_val = $self->{significance}; |
739
|
0
|
|
|
|
|
|
my $standard_F = fdistr ($df_bg_Y,$df_wg_Y_Adj,$chosen_p_val); |
740
|
|
|
|
|
|
|
|
741
|
|
|
|
|
|
|
#/ this approach confuses people! have it simply pass of fail at their selected p_value - it already has a default value |
742
|
|
|
|
|
|
|
# my $message = $F > $standard_F ? qq{This value of F is significant at your chosen p = $chosen_p_val level. } |
743
|
|
|
|
|
|
|
# : $F > fdistr ($df_bg_Y,$df_wg_Y_Adj,0.01) ? qq{This value of F is significant at the p = 0.01 level. } |
744
|
|
|
|
|
|
|
# : $F > fdistr ($df_bg_Y,$df_wg_Y_Adj,0.05) ? qq{This value of F is significant at the p = 0.05 level. } |
745
|
|
|
|
|
|
|
# : qq{This is not a significant value of F. } # default behaviour |
746
|
|
|
|
|
|
|
# ; |
747
|
|
|
|
|
|
|
|
748
|
0
|
0
|
|
|
|
|
my $message = $F > $standard_F ? qq{This value of F is significant at the p = $chosen_p_val level. } : |
749
|
|
|
|
|
|
|
qq{This is not a significant value of F at the p = $chosen_p_val level. }; |
750
|
|
|
|
|
|
|
|
751
|
|
|
|
|
|
|
|
752
|
0
|
|
|
|
|
|
my $p_for_F = fprob ( $df_bg_Y, $df_wg_Y_Adj, $F ) ; |
753
|
0
|
|
|
|
|
|
$self->{output}{message} = $message; |
754
|
0
|
|
|
|
|
|
$self->{output}{standard_F} = $standard_F; |
755
|
0
|
|
|
|
|
|
$self->{output}{p_for_F} = $p_for_F; |
756
|
|
|
|
|
|
|
|
757
|
|
|
|
|
|
|
return ( |
758
|
|
|
|
|
|
|
#VOID { $self->_print_form(@other_args) } |
759
|
0
|
|
|
0
|
|
|
VOID { $self->_print_form() } |
760
|
|
|
|
|
|
|
|
761
|
|
|
|
|
|
|
# LIST { ( sprintf (qq{%.3f},$F), sprintf (qq{%.3f},$p_for_F), sprintf (qq{%.3f},$self->{output}{MS_bg}), |
762
|
|
|
|
|
|
|
# sprintf (qq{%.3f},$self->{SS}{Y}{between_group_adjusted}), $df_bg_Y, |
763
|
|
|
|
|
|
|
# sprintf (qq{%.3f},$self->{output}{MS_wg}), sprintf (qq{%.3f},$self->{SS}{Y}{within_group_adjusted}), |
764
|
|
|
|
|
|
|
# $df_wg_Y_Adj, sprintf (qq{%.3f},$self->{SS}{Y}{total_adjusted}), ) } |
765
|
|
|
|
|
|
|
|
766
|
0
|
|
|
0
|
|
|
LIST { $self->_return_list(@other_args) } |
767
|
0
|
0
|
|
0
|
|
|
BOOL { $F > $standard_F ? 1 : undef; } |
768
|
0
|
|
|
0
|
|
|
NUM { $F ; } |
769
|
0
|
|
|
0
|
|
|
STR { $message } |
770
|
0
|
|
|
|
|
|
); |
771
|
|
|
|
|
|
|
} |
772
|
|
|
|
|
|
|
=head2 results |
773
|
|
|
|
|
|
|
|
774
|
|
|
|
|
|
|
Used to access the results of the ANCOVA analysis. This method is context-dependent and will return a variety of |
775
|
|
|
|
|
|
|
different values depending on its calling context. In VOID context prints a report to STDOUT (use |
776
|
|
|
|
|
|
|
C to print more detailed report). |
777
|
|
|
|
|
|
|
|
778
|
|
|
|
|
|
|
# To print a short report to STDOUT |
779
|
|
|
|
|
|
|
$anc->results(); |
780
|
|
|
|
|
|
|
# To print a detailed report set output_verbosity to 1 on object creation or using the set_output_verbosity> method. |
781
|
|
|
|
|
|
|
$anc->set_output_verbosity(1); |
782
|
|
|
|
|
|
|
$anc->results(); |
783
|
|
|
|
|
|
|
|
784
|
|
|
|
|
|
|
|
785
|
|
|
|
|
|
|
In LIST context it either returns the full list of all relevant values of F, p, df, MS... or returns an ordered subset of the values |
786
|
|
|
|
|
|
|
depending on whether you call it without or with numbered arguments respectively (see below). |
787
|
|
|
|
|
|
|
|
788
|
|
|
|
|
|
|
# Calling results in LIST without arguments returns the full list of relevant values of F, p, df, MS... |
789
|
|
|
|
|
|
|
my %hash; |
790
|
|
|
|
|
|
|
@hash{qw($F_score, $p_value, $MS_bg, $SS_bg_Adj, $df_bg_Y, $MS_wg, $SS_wg_Adj, $df_wg_Y_Adj, $SS_total_Adj)} = $anc->results(); |
791
|
|
|
|
|
|
|
for (keys %hash) { print qq{\n$_ = $hash{$_} } }; |
792
|
|
|
|
|
|
|
|
793
|
|
|
|
|
|
|
However, calling C in LIST context with numbered arguments corresponding to those below returns those arguments |
794
|
|
|
|
|
|
|
in the order passed to the method. |
795
|
|
|
|
|
|
|
|
796
|
|
|
|
|
|
|
# 0 1 2 3 4 5 6 7 8 |
797
|
|
|
|
|
|
|
# ($F_score, $p_value, $MS_bg, $SS_bg_Adj, $df_bg_Y, $MS_wg, $SS_wg_Adj, $df_wg_Y_Adj, $SS_total_Adj) = $anc->results(2,3,5) |
798
|
|
|
|
|
|
|
print qq{\n\nCalling in LIST context. The F value, p_value, MS_bg and MS_wg are: @{$anc->results(0,1,2,,5)}}; |
799
|
|
|
|
|
|
|
|
800
|
|
|
|
|
|
|
In BOOLEAN context it returns true or false depending on whether the obtained F score was significant at the p_value chosen |
801
|
|
|
|
|
|
|
upon object creation or set using the C method (defaults to p = 0.05). |
802
|
|
|
|
|
|
|
|
803
|
|
|
|
|
|
|
if ($anc->results) { print qq{\nThis result is significant.} } else { print qq{\nThis result is not significant.} } |
804
|
|
|
|
|
|
|
|
805
|
|
|
|
|
|
|
In STRING context it returns a string message about whether the obtained F score was significant at the chosen p_value. |
806
|
|
|
|
|
|
|
|
807
|
|
|
|
|
|
|
print qq{\n\nCall result in string returns a message : }, ''.$anc->results; # Prints 'This value of F is significant at your chosen .05 level'... |
808
|
|
|
|
|
|
|
|
809
|
|
|
|
|
|
|
=cut |
810
|
|
|
|
|
|
|
|
811
|
|
|
|
|
|
|
sub _print_form { |
812
|
|
|
|
|
|
|
#my ( $self, $verbose ) = @_; |
813
|
0
|
|
|
0
|
|
|
my $self = shift; |
814
|
|
|
|
|
|
|
|
815
|
|
|
|
|
|
|
#$verbose ||= 0; |
816
|
|
|
|
|
|
|
#my $self = shift; |
817
|
|
|
|
|
|
|
#my $verbose = shift ||= 0; |
818
|
|
|
|
|
|
|
#$verbose = $verbose eq q{verbose} ? 1 : 0 ; |
819
|
|
|
|
|
|
|
|
820
|
0
|
|
|
|
|
|
my $verbose = $self->{verbosity}{output}; |
821
|
|
|
|
|
|
|
|
822
|
0
|
|
|
|
|
|
my @groups = @{$self->{groups}}; |
|
0
|
|
|
|
|
|
|
823
|
|
|
|
|
|
|
|
824
|
|
|
|
|
|
|
#$verbose or print form { bullet => q{*} }, |
825
|
0
|
0
|
|
|
|
|
$verbose and print form { bullet => q{*} }, |
826
|
|
|
|
|
|
|
qq{\n\n ============================================================================= }, |
827
|
|
|
|
|
|
|
qq{| Sum of squared | Sum of squared | Sum of co-deviates |}, |
828
|
|
|
|
|
|
|
qq{| deviates for X | deviates for Y | |}, |
829
|
|
|
|
|
|
|
qq{|-------------------------|-------------------------|-------------------------|}, |
830
|
|
|
|
|
|
|
qq{| SS_T_x = {<<<<<<<<<<<} | SS_T_y = {<<<<<<<<<<<} | SC_T = {<<<<<<<<<<<<<} |}, |
831
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{X}{total}), sprintf (qq{%.3f},$self->{SS}{Y}{total}), sprintf (qq{%.3f},$self->{SC}{T}{SC_within_group}), |
832
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
833
|
|
|
|
|
|
|
qq{| SS_wg_x = {<<<<<<<<<<<} | SS_wg_y = {<<<<<<<<<<<} | SC_wg = {<<<<<<<<<<<<<} |}, |
834
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{X}{within_group}), sprintf (qq{%.3f},$self->{SS}{Y}{within_group}), sprintf (qq{%.3f},$self->{SC}{T}{sum_of_X_and_Y_SC_within_group}), |
835
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
836
|
|
|
|
|
|
|
qq{| | SS_bg_y = {<<<<<<<<<<<} | |}, |
837
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{Y}{between_group}), |
838
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
839
|
|
|
|
|
|
|
qq{ }, |
840
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
841
|
|
|
|
|
|
|
qq{| Overall correlation and adjustment of SS_T_Y |}, |
842
|
|
|
|
|
|
|
qq{|-------------------------|-------------------------|-------------------------|}, |
843
|
|
|
|
|
|
|
qq{| r_T = {<<<<<<<<<<<<<<<} | r_T_sq = {<<<<<<<<<<<<} | SS_T_y_Adj = {<<<<<<<} |}, |
844
|
|
|
|
|
|
|
sprintf (qq{%.9f},$self->{output}{r_T}), sprintf (qq{%.9f},$self->{output}{r_T_sq}), sprintf (qq{%.3f},$self->{SS}{Y}{total_adjusted}), |
845
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
846
|
|
|
|
|
|
|
qq{| Aggregate correlation and adjustment of SS_wg_y |}, |
847
|
|
|
|
|
|
|
qq{|-------------------------|-------------------------|-------------------------|}, |
848
|
|
|
|
|
|
|
qq{| r_wg = {<<<<<<<<<<<<<<} | r_wg_sq = {<<<<<<<<<<<} | SS_wg_y_Adj = {<<<<<<<} |}, |
849
|
|
|
|
|
|
|
sprintf (qq{%.9f},$self->{output}{r_wg}), sprintf (qq{%.9f},$self->{output}{r_wg_sq}), sprintf (qq{%.3f},$self->{SS}{Y}{within_group_adjusted}), |
850
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
851
|
|
|
|
|
|
|
qq{| Adjustment of SS_bg_y |}, |
852
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
853
|
|
|
|
|
|
|
qq{| | SS_bg_y_Adj = {<<<<<<<} | |}, # this specifies locations and formating of variables |
854
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{Y}{between_group_adjusted}), |
855
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
856
|
|
|
|
|
|
|
qq{ }, |
857
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
858
|
|
|
|
|
|
|
qq{| Overall Means for X and Y variables |}, |
859
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
860
|
|
|
|
|
|
|
qq{| Mean_X_overall = {<<<<<<<<<<<<<<<<<<} | Mean_Y_overall = {<<<<<<<<<<<<<<<<} |}, |
861
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{T}{X}{mean}), sprintf (qq{%.3f},$self->{SS}{T}{Y}{mean}), |
862
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
863
|
|
|
|
|
|
|
; |
864
|
|
|
|
|
|
|
|
865
|
|
|
|
|
|
|
#if (!$verbose) { |
866
|
0
|
0
|
|
|
|
|
if ($verbose) { |
867
|
0
|
|
|
|
|
|
for my $group ( (sort {$a cmp $b} @groups) ) { |
|
0
|
|
|
|
|
|
|
868
|
0
|
|
|
|
|
|
print form { bullet => q{*} }, |
869
|
|
|
|
|
|
|
qq{ }, |
870
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
871
|
|
|
|
|
|
|
qq{| Means for group {[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[} |}, |
872
|
|
|
|
|
|
|
$group, |
873
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
874
|
|
|
|
|
|
|
qq{| Mean_X = {<<<<<<<<<<<<} | Mean_Y = {<<<<<<<<<<<<} | Mean_Y_Adj = {<<<<<<<<} |}, |
875
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{$group}{X}{mean}), sprintf (qq{%.3f},$self->{SS}{$group}{Y}{mean}), sprintf (qq{%.3f},$self->{SS}{$group}{Y}{mean_adjusted}), |
876
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
877
|
|
|
|
|
|
|
; |
878
|
|
|
|
|
|
|
} |
879
|
|
|
|
|
|
|
} |
880
|
|
|
|
|
|
|
|
881
|
0
|
|
|
|
|
|
print form { bullet => q{*} }, |
882
|
|
|
|
|
|
|
qq{\n ============================================================================= }, |
883
|
|
|
|
|
|
|
qq{| ANCOVA |}, |
884
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
885
|
|
|
|
|
|
|
qq{| | df | SS | MS | F | p |}, |
886
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
887
|
|
|
|
|
|
|
qq{| Adjusted means be- | {<<<<<<} | {<<<<<<} | {<<<<<<} | {<<<<<<} | {<<<<<<<<} |}, |
888
|
|
|
|
|
|
|
$self->{output}{df_Y_bg}, sprintf (qq{%.3f},$self->{SS}{Y}{between_group_adjusted}), |
889
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{output}{MS_bg}), sprintf (qq{%.3f},$self->{output}{F_score}), $self->{output}{p_for_F}, |
890
|
|
|
|
|
|
|
qq{| teen groups effect | | | | | |}, |
891
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
892
|
|
|
|
|
|
|
qq{| Adjusted error | {<<<<<<} | {<<<<<<} | {<<<<<<} | | |}, |
893
|
|
|
|
|
|
|
$self->{output}{df_Y_wg_Adj}, sprintf (qq{%.3f},$self->{SS}{Y}{within_group_adjusted}), sprintf (qq{%.3f},$self->{output}{MS_wg}), |
894
|
|
|
|
|
|
|
qq{| within groups | | | | | |}, |
895
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
896
|
|
|
|
|
|
|
qq{| Adjusted total | | {<<<<<<} | | | |}, |
897
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{Y}{total_adjusted}), |
898
|
|
|
|
|
|
|
qq{| | | | | | |}, |
899
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
900
|
|
|
|
|
|
|
; |
901
|
|
|
|
|
|
|
|
902
|
0
|
0
|
|
|
|
|
$verbose and print form { bullet => q{*} }, |
903
|
|
|
|
|
|
|
qq{ }, |
904
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
905
|
|
|
|
|
|
|
qq{| Overview |}, |
906
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
907
|
|
|
|
|
|
|
qq{| your chosen p value = {<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<} |}, |
908
|
|
|
|
|
|
|
$self->{significance}, |
909
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
910
|
|
|
|
|
|
|
qq{| standard F value for these df and p = {<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<} |}, |
911
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{output}{standard_F}), |
912
|
|
|
|
|
|
|
qq{|-----------------------------------------------------------------------------|}, |
913
|
|
|
|
|
|
|
qq{| {[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[} |}, |
914
|
|
|
|
|
|
|
$self->{output}{message}, |
915
|
|
|
|
|
|
|
qq{ ============================================================================= }, |
916
|
|
|
|
|
|
|
; |
917
|
|
|
|
|
|
|
|
918
|
0
|
|
|
|
|
|
return; |
919
|
|
|
|
|
|
|
} |
920
|
|
|
|
|
|
|
|
921
|
|
|
|
|
|
|
sub _return_list { |
922
|
|
|
|
|
|
|
|
923
|
0
|
|
|
0
|
|
|
my @list = @_; |
924
|
0
|
|
|
|
|
|
my $self = shift @list; |
925
|
|
|
|
|
|
|
# unpack the rest of @_ |
926
|
|
|
|
|
|
|
#print qq{\n\nlist is: @list}; |
927
|
0
|
|
|
|
|
|
my @returns = (); |
928
|
|
|
|
|
|
|
#if ( scalar(@list) > 0 ) { |
929
|
0
|
0
|
|
|
|
|
if ( scalar(@list) ) { |
930
|
|
|
|
|
|
|
#while (my $parameter = @list) |
931
|
|
|
|
|
|
|
|
932
|
0
|
|
|
|
|
|
my %mapping = ( 0 => q{F_score}, |
933
|
|
|
|
|
|
|
1 => q{p_for_F}, |
934
|
|
|
|
|
|
|
2 => q{MS_bg}, |
935
|
|
|
|
|
|
|
3 => q{between_group_adjusted}, |
936
|
|
|
|
|
|
|
4 => q{df_Y_bg}, |
937
|
|
|
|
|
|
|
5 => q{MS_wg}, |
938
|
|
|
|
|
|
|
6 => q{within_group_adjusted}, |
939
|
|
|
|
|
|
|
7 => q{df_Y_wg_Adj}, |
940
|
|
|
|
|
|
|
8 => q{total_adjusted}, |
941
|
|
|
|
|
|
|
); |
942
|
|
|
|
|
|
|
|
943
|
|
|
|
|
|
|
#print qq{\nhere is the list @list}; |
944
|
|
|
|
|
|
|
|
945
|
0
|
|
|
|
|
|
for my $parameter (@list) { |
946
|
|
|
|
|
|
|
|
947
|
0
|
0
|
|
|
|
|
croak qq{\nThe parameters passed must be numeric corresponding to those documented in the synopsis.} if ($parameter !~ /\A[0-8]\z/xms); |
948
|
|
|
|
|
|
|
|
949
|
|
|
|
|
|
|
#print qq{\npushing $parameter}; |
950
|
|
|
|
|
|
|
|
951
|
0
|
|
|
|
|
|
my $named_param = $mapping{$parameter}; |
952
|
|
|
|
|
|
|
#print qq{\nmy named $named_param}; |
953
|
|
|
|
|
|
|
|
954
|
0
|
0
|
0
|
|
|
|
if ( $parameter == 3 || $parameter == 6 || $parameter == 8 ) { |
|
|
|
0
|
|
|
|
|
955
|
|
|
|
|
|
|
|
956
|
0
|
|
|
|
|
|
push @returns, sprintf(qq{%.3f},$self->{SS}{Y}{$named_param}); |
957
|
|
|
|
|
|
|
} |
958
|
|
|
|
|
|
|
else { |
959
|
|
|
|
|
|
|
|
960
|
0
|
|
|
|
|
|
my $value = $self->{output}{$named_param}; |
961
|
0
|
0
|
|
|
|
|
$value = sprintf(qq{%.3f},$value) if ($parameter != 1); |
962
|
0
|
|
|
|
|
|
push @returns, $value; |
963
|
|
|
|
|
|
|
} |
964
|
|
|
|
|
|
|
|
965
|
|
|
|
|
|
|
} |
966
|
0
|
|
|
|
|
|
return @returns; |
967
|
|
|
|
|
|
|
} |
968
|
|
|
|
|
|
|
|
969
|
|
|
|
|
|
|
else { |
970
|
|
|
|
|
|
|
|
971
|
0
|
|
|
|
|
|
@returns = ( sprintf (qq{%.3f},$self->{output}{F_score}), $self->{output}{p_for_F}, |
972
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{output}{MS_bg}), sprintf (qq{%.3f},$self->{SS}{Y}{between_group_adjusted}), |
973
|
|
|
|
|
|
|
$self->{output}{df_Y_bg}, sprintf (qq{%.3f},$self->{output}{MS_wg}), |
974
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{Y}{within_group_adjusted}), $self->{output}{df_Y_wg_Adj}, |
975
|
|
|
|
|
|
|
sprintf (qq{%.3f},$self->{SS}{Y}{total_adjusted}) ); |
976
|
|
|
|
|
|
|
|
977
|
|
|
|
|
|
|
} |
978
|
0
|
|
|
|
|
|
return @returns; |
979
|
|
|
|
|
|
|
} |
980
|
|
|
|
|
|
|
|
981
|
|
|
|
|
|
|
1; # Magic true value required at end of module |
982
|
|
|
|
|
|
|
__END__ |