File Coverage

blib/lib/Statistics/MVA/MultipleRegression.pm

Criterion	Covered	Total	%
statement	21	71	29.5
branch	0	10	0.0
condition			n/a
subroutine	7	8	87.5
pod	0	1	0.0
total	28	90	31.1

package Statistics::MVA::MultipleRegression;

use warnings;
use strict;
use Carp;
# re-implement with Cepehes?
use Statistics::MVA;
use Math::MatrixReal;
use List::Util qw/sum/;

require Exporter;
our @ISA = qw(Exporter);
our @EXPORT = qw(linear_regression);

=head1 NAME

Statistics::MVA::MultipleRegression - Simple Least Squares Linear Multiple Regression Module.

=cut
=head1 VERSION

This document describes Statistics::MVA::MultipleRegression version 0.0.1

=cut
=head1 SYNOPSIS

    my $lol = [
        [qw/745  36  66/],
        [qw/895  37  68/],
        [qw/442  47  64/],
        [qw/440  32  53/],
        [qw/1598 1   101/],
     ];

    use Statistics::MVA::MultipleRegression;

    # Call exported routine on List-of-List of data (each nested list - row - corresponds to a set of observations. In void-context it prints a report to STDOUT.
    linear_regression($lol);

    # In LIST-context the routine returns an ARRAY reference of the coefficients and R^2.
    my ($Array_ref_of_coefficients, $R_sq) = linear_regression($d);

=cut
=head1 DESCRIPTION

    The general purpose of multiple regression is to gain information about the relationship between several independent
    variables (x_i) and a dependent variable (y). The procedure involves fitting an equation of the form:

        y = b_o + x_1 * b_1 + x_2 * b_2 +... ... x_n * b_n

    This is a minimal implementation of least squares procedure that accepts a List-of-Lists in which each nested list
    corresponds to a single set of observations. The single exported routine returns the coefficients (b_i) and R^2.

=cut

use version; our $VERSION = qv('0.0.1');


#/ at this moment this is a minimal version that just calculates the coefficients. i will expand it there is interest.

sub linear_regression {
    my $lol = shift;

    croak qq{\nAll data must be a matrix of numbers} if &Statistics::MVA::CVMat::_check($lol);

    $lol = &Statistics::MVA::CVMat::transpose($lol);

    croak qq{\nYou have no dependent variables - if you want to calculate the mean of Y do not use this module.} if (scalar @{$lol} == 1);
    carp qq{\nYou have one dependent variables - why use a multivariate procedure?} if (scalar @{$lol} == 2);

    my $y = [@{$lol->[0]}];
    my $n = scalar @{$y};
    my $y_mean = (sum @{$y}) / $n;
    my $x = [@{$lol}[1..$#{$lol}]];

    my $y_mat = Math::MatrixReal->new_from_cols([$lol->[0]]);

    #y can either build a vector of 1´s using x - i.e. 1 x scalar $lol1->[0] and combine it with everything after ->[0]
    #y or can simply in place modify the whole of ->[0] in place now that it´s been copied to a matrix

    # in place modification
    for my $i (0..$#{$lol->[0]}) { $lol->[0][$i] = 1 }

    my $x_mat = Math::MatrixReal->new_from_cols($lol);
    my $b_mat = $y_mat->shadow;

    $b_mat =  (((~$x_mat * $x_mat)**-1) * ~$x_mat) * $y_mat ;

    # constant is $Bs[0] then coefficients are @Bs[1..$#Bs]
    my @Bs = map { $b_mat->[0][$_][0] } (0..$#{$b_mat->[0]});

    my $ss_y = sum map { ($_-$y_mean)**2 } @{$y};

    my @Bs_copy = @Bs;
    my $con = shift @Bs_copy;

    my $y_prime = [];

    #for my $r (0..$#{$x->[0]}) {
    for my $r (0..$#{$y}) {
        my $val = sum $con, map { $Bs_copy[$_] * $x->[$_][$r] } (0..$#Bs_copy);
        push @{$y_prime }, $val;
    }

    #my $ss_error = sum map { ($y->[$_]-$y_prime->[$_])**2 } (0..$#{$y});
    # can use ss_y_prime / ss_y or (ss_y - ss_error) / ss_y
    my $ss_y_prime = sum map { ($_-$y_mean)**2 } @{$y_prime};

    my $R2 = $ss_y_prime/$ss_y;

    if (!wantarray) {
        print qq{\nThe coefficients are: };
        for (0..$#Bs) {
            my $punct = $_ == $#Bs ? q{.} : q{, };
            print qq{B[$_] = $Bs[$_]$punct};
        }
        print qq{\nR^2 is $R2};
        return;
    }
    else { return ([@Bs], $R2); }
}

1; # Magic true value required at end of module
__END__