<?php
require_once dirname(__DIR__) . '/base/regression.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/basic/describe.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/distribution/studentT.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/distribution/f.php';
require_once dirname(dirname(__DIR__)) . '/zendasmath/tools.php';
require_once __DIR__ . '/regressionMatrix.php';

class GeneralLinearModel extends RegressionBase
{
    public $factors;

    public $covariate;

    public $response;

    public $observeN;

    public $hasFactor;

    public $hasCovariate;

    public $factorLevel;

    public $dataframe;

    public $freedoms;

    public $coefN;

    public function __construct($dataframe, $factors, $covariate, $response)
    {
        $this->factors   = $factors;
        $this->covariate = $covariate;
        $this->response  = $response;

        $this->observeN = count($response);

        $this->hasFactor    = count($factors) > 0;
        $this->hasCovariate = count($covariate) > 0;

        if($this->hasFactor) $this->factorLevel = $this->getGroup();

        parent::__construct($dataframe);
    }

    public function getSingleFormula($yaxis)
    {
        $params     = $this->coef();
        $common     = current($params);
        $covariates = array();
        $factors    = array();

        foreach($this->covariate as $name => $values)
        {
            $covariates[$name] = next($params);
        }

        foreach($this->factorLevel as $name => $fLv)
        {
            $factors[$name] = array();
            $items = $fLv->items;
            $last  = array_pop($items);
            foreach($items as $item)
            {
                $factors[$name][$name . '_' . $item] = next($params);
            }
            $factors[$name][$name . '_' . $last] = -1 * array_sum($factors[$name]);
        }

        $part1 = $this->dataframe->roundDigit($common);

        $part2 = '';

        foreach($covariates as $name => $c)
        {
            $mark = $c >= 0 ? ' + ' : ' - ';
            $c = abs($c);
            $c = $this->dataframe->roundDigit($c);
            $part2 .= ' ' . $mark . $c . ' ' . $name;
        }

        foreach($factors as $f)
        {
            foreach($f as $name => $c)
            {
                $mark = $c >= 0 ? ' + ' : ' - ';
                $c = abs($c);
                $c = $this->dataframe->roundDigit($c);
                $part2 .= ' ' . $mark . $c . ' ' . $name;
            }
        }

        return $yaxis . ' = ' . $part1 . $part2;
    }

    public function getFormula($yaxis)
    {
        $params     = $this->coef();
        $common     = current($params);
        $covariates = array();
        $factors    = array();

        if($this->hasCovariate)
        {
            foreach($this->covariate as $name => $values)
            {
                $covariates[$name] = next($params);
            }
        }

        if($this->hasFactor)
        {
            foreach($this->factorLevel as $name => $fLv)
            {
                $factors[$name] = array();
                $items = $fLv->items;
                $last  = array_pop($items);
                foreach($items as $item)
                {
                    $factors[$name][$name . '_' . $item] = next($params);
                }
                $factors[$name][$name . '_' . $last] = -1 * array_sum($factors[$name]);
            }
        }

        $formulas = array();
        if($this->hasCovariate and $this->hasFactor)
        {
            $part1 = '';
            foreach($covariates as $name => $c)
            {
                $mark = $c >= 0 ? ' + ' : ' - ';
                $c = abs($c);
                $c = $this->dataframe->roundDigit($c);
                $part1 .= ' ' . $mark . $c . ' ' . $name;
            }

            $factorCoef = current($factors);
            $nextCoef   = next($factors);
            while($nextCoef)
            {
                $newCoef = array();
                foreach($factorCoef as $name1 => $c1)
                {
                    foreach($nextCoef as $name2 => $c2)
                    {
                        $newCoef[$name1 . ' ' . $name2] = $c1 + $c2;
                    }
                }
                $factorCoef = $newCoef;
                $nextCoef   = next($factors);
            }

            foreach($factorCoef as $name => $c)
            {
                $c += $common;
                $mark = $c >= 0 ? '' : ' - ';
                $c = abs($c);
                $c = $this->dataframe->roundDigit($c);
                $formulas[$name] = $yaxis . ' = ' . $mark . $c . $part1;
            }
        }
        else if($this->hasFactor)
        {
            $formula = '';
            foreach($factors as $factor)
            {
                foreach($factor as $name => $c)
                {
                    $mark = $c >= 0 ? ' + ' : ' - ';
                    $c = abs($c);
                    $c = $this->dataframe->roundDigit($c);
                    $formula .= $mark . $c . ' ' . $name;
                }
            }
            $common = $this->dataframe->roundDigit($common);
            $formulas[] = $yaxis . ' = ' . $common .  $formula;
        }
        else
        {
            $formula = '';
            foreach($covariates as $name => $c)
            {
                $mark = $c >= 0 ? ' + ' : ' - ';
                $c = abs($c);
                $c = $this->dataframe->roundDigit($c);
                $formula .= $mark . $c . ' ' . $name;
            }
            $common = $this->dataframe->roundDigit($common);
            $formulas[] = $yaxis . ' = ' . $common .  $formula;
        }

        return $formulas;
    }

    public function getGroup()
    {
        $result = array();
        $this->freedoms = array();
        foreach($this->factors as $name => $factor)
        {
            $values = array_count_values($factor);
            $keys   = array_keys($values);
            utf8_array_asort($keys);
            $keys = array_values($keys);
            $result[$name] = new stdclass();
            $result[$name]->level = count($keys);
            $result[$name]->items = $keys;

            $this->freedoms[] = count($keys) - 1;
        }

        return $result;
    }


    public function checkInvertible()
    {
        $X      = $this->X();
        $Xt     = $X->transpose(); // Xᵀ

        return $Xt->multiply($X)->isInvertible();
    }

    public function getFactor($round = true)
    {
        $X      = $this->X();
        $Xt     = $X->transpose(); // Xᵀ
        $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        $Y      = RegressionMatrix::yMatrix($this->response);

        $coef = $this->coef($Xt, $XtXinv, $Y);
        $B    = RegressionMatrix::coefMatrix($coef);

        $this->coefN = count($coef);

        $ssError = $this->sse($X, $Xt, $Y, $B);
        $dfError = $this->ssedf();
        $mse     = $this->mse($ssError, $dfError);

        $se  = $this->se($XtXinv, $mse);

        $t = $this->tFactor($coef, $se);
        $p = $this->pFactor($dfError, $t);

        // $vif = $this->vif();
        // array_unshift($vif, '');

        $result = new stdclass();
        $result->coef = array();
        $result->se   = array();
        $result->t    = array();
        $result->p    = array();
        // $result->vif  = array();

        foreach($coef as $index => $c)
        {
            $result->coef[] = $this->dataframe->roundDigit($c);
            $result->se[]   = $this->dataframe->roundDigit($se[$index]);
            $result->t[]    = $this->dataframe->roundDigit($t[$index]);
            $result->p[]    = $this->dataframe->roundDigit($p[$index]);
            // $result->vif[]  = $this->dataframe->roundDigit($vif[$index]);
        }

        return $result;
    }

    public function getModel($round = true)
    {
        $X      = $this->X();
        $Xt     = $X->transpose(); // Xᵀ
        $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        $Y      = RegressionMatrix::yMatrix($this->response);
        $coef   = $this->coef($Xt, $XtXinv, $Y);
        $B      = RegressionMatrix::coefMatrix($coef);
        $C      = $Xt->multiply($X)->inverse();
        $H      = $X->multiply($C)->multiply($Xt); //  X⟮XᵀX⟯⁻¹Xᵀ
        $fitY   = $this->fitY($coef);

        $dfError = $this->ssedf();
        $dfTotal = $this->sstdf();
        $ssError = $this->sse($X, $Xt, $Y, $B);
        $ssTotal = $this->sst($Y);
        $msError = $this->mse($ssError, $dfError);
        $msTotal = $this->mst($ssTotal, $dfTotal);

        $s      = self::is_num($msError) ? sqrt($msError) : '-';
        $rsq    = $this->rsq($ssError, $ssTotal);
        $rsqAdj = $this->rsqAdj($msError, $msTotal);
        $rsqPre = $this->rsqPre($H, $fitY);

        if(self::is_num($rsq) and $rsq < 0) $rsq = 0;
        if(self::is_num($rsqAdj) and $rsqAdj < 0) $rsqAdj = 0;

        $s      = $this->dataframe->roundDigit($s, $round);
        $rsq    = $this->dataframe->roundDigit($rsq, $round);
        $rsqAdj = $this->dataframe->roundDigit($rsqAdj, $round);
        $rsqPre = $this->dataframe->roundDigit($rsqPre, $round);

        $result  = new stdclass();

        $result->s      = $s;
        $result->rsq    = self::is_num($rsq) ? $rsq * 100 : $rsq;
        $result->rsqAdj = self::is_num($rsqAdj) ? $rsqAdj * 100 : $rsqAdj;
        $result->rsqPre = self::is_num($rsqPre) ? $rsqPre * 100 : $rsqPre;

        return $result;
    }

    public function getAnova($round = true)
    {
        $ssrdf = $this->ssrdf();
        $ssedf = $this->ssedf($ssrdf);
        $sstdf = $this->sstdf();

        $X      = $this->X();
        $Xt     = $X->transpose(); // Xᵀ
        $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        $Y      = RegressionMatrix::yMatrix($this->response);
        $coef   = $this->coef($Xt, $XtXinv, $Y);
        $B      = RegressionMatrix::coefMatrix($coef);

        $ssr = $this->ssr($X, $Xt, $Y, $B);
        $sse = $this->sse($X, $Xt, $Y, $B);
        $sst = $this->sst($Y);

        $msr = $this->msr($ssr, $ssrdf);
        $mse = $this->mse($sse, $ssedf);

        $fr = $this->fr($msr, $mse);
        $pr = $this->pr($ssrdf, $ssedf, $fr);

        $same  = $this->sameItem();
        $pedf  = '-';
        $lofdf = '-';
        $sspe  = '-';
        $sslof = '-';
        $mspe  = '-';
        $mslof = '-';
        $flof  = '-';
        $plof  = '-';
        if(!empty($same))
        {
            $pedf  = $this->pedf($same);
            $lofdf = $this->lofdf($ssedf, $pedf);

            $sspe  = $this->sspe($same);
            $sslof = $this->sslof($sse, $sspe);

            $mspe = $this->mspe($sspe, $pedf);
            $mslof = $this->mslof($sslof, $lofdf);

            if($mspe != 0)
            {
                $flof  = $this->flof($mslof, $mspe);
                $plof  = $this->plof($lofdf, $pedf, $flof);
            }
        }

        $result = new stdclass();

        $result->hasSame = !empty($same);

        $result->df      = $ssrdf;
        $result->dfError = $ssedf;
        $result->dfLOF   = $lofdf;
        $result->dfPE    = $pedf;
        $result->dfTotal = $sstdf;

        $result->ssRegr  = $this->dataframe->roundDigit($ssr, $round);
        $result->ssError = $this->dataframe->roundDigit($sse, $round);
        $result->ssTotal = $this->dataframe->roundDigit($sst, $round);
        $result->ssPE    = $this->dataframe->roundDigit($sspe, $round);
        $result->ssLOF   = $this->dataframe->roundDigit($sslof, $round);

        $result->msRegr  = $this->dataframe->roundDigit($msr, $round);
        $result->msError = $this->dataframe->roundDigit($mse, $round);
        $result->msPE    = $this->dataframe->roundDigit($mspe, $round);
        $result->msLOF   = $this->dataframe->roundDigit($mslof, $round);

        $result->fRegr = $this->dataframe->roundDigit($fr, $round);
        $result->pRegr = $this->dataframe->roundDigit($pr, $round);
        $result->fLOF  = $this->dataframe->roundDigit($flof, $round);
        $result->pLOF  = $this->dataframe->roundDigit($plof, $round);

        $result->itemDf = array();
        $result->itemSS = array();
        $result->itemMS = array();
        $result->itemF = array();
        $result->itemP = array();

        if($this->hasCovariate)
        {
            foreach($this->covariate as $name => $covariate)
            {
                $covariates = $this->covariate;
                if(count($covariates) > 1 or $this->hasFactor) unset($covariates[$name]);
                $X = $this->X($this->factors, $covariates);

                $itemDF    = 1;
                $itemSSR   = $ssr - $this->ssr($X);
                $itemMSR   = $this->msr($itemSSR, $itemDF);
                $itemFR    = (self::is_num($itemMSR) and self::is_num($mse) and $mse != 0) ? $itemMSR / $mse : '-';
                $fdist     = new F($itemDF, $ssedf);
                $cdf       = $fdist->cdf($itemFR);
                $itemPR    = self::is_num($cdf) ? 1 - $cdf : '-';

                $result->itemDf[] = $itemDF;
                $result->itemSS[] = $this->dataframe->roundDigit($itemSSR, $round);
                $result->itemMS[] = $this->dataframe->roundDigit($itemMSR, $round);
                $result->itemF[]  = $this->dataframe->roundDigit($itemFR, $round);
                $result->itemP[]  = $this->dataframe->roundDigit($itemPR, $round);
            }
        }

        if($this->hasFactor)
        {
            foreach($this->factors as $name => $factor)
            {
                $factors = $this->factors;
                if(count($factors)  > 1 or $this->hasCovariate) unset($factors[$name]);
                $X = $this->X($factors, $this->covariate);

                $itemDF    = count(array_count_values($factor)) - 1;
                $itemSSR   = $ssr - $this->ssr($X);
                $itemMSR   = $this->msr($itemSSR, $itemDF);
                $itemFR    = (self::is_num($mse) and $mse != 0 and self::is_num($itemMSR)) ? $itemMSR / $mse : '-';
                $fdist     = self::is_num($itemFR) ? new F($itemDF, $ssedf) : '-';
                $itemPR    = self::is_num($fdist) ? 1 - $fdist->cdf($itemFR) : '-';

                $result->itemDf[] = $itemDF;
                $result->itemSS[] = $this->dataframe->roundDigit($itemSSR, $round);
                $result->itemMS[] = $this->dataframe->roundDigit($itemMSR, $round);
                $result->itemF[]  = $this->dataframe->roundDigit($itemFR, $round);
                $result->itemP[]  = $this->dataframe->roundDigit($itemPR, $round);
            }
        }

        return $result;
    }

    public function getTest($round = true)
    {
        $X      = $this->X();
        $Xt     = $X->transpose(); // Xᵀ
        $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        $Y      = RegressionMatrix::yMatrix($this->response);
        $coef   = $this->coef($Xt, $XtXinv, $Y);
        $B      = RegressionMatrix::coefMatrix($coef);
        $C      = $Xt->multiply($X)->inverse();
        $H      = $X->multiply($C)->multiply($Xt); //  X⟮XᵀX⟯⁻¹Xᵀ
        $h      = $H->getDiagonalElements();

        $fitY     = $this->fitY($coef);
        $y        = $this->response;
        $residual = $this->residual($fitY);
        $dfError  = $this->ssedf();
        $ssError  = $this->sse($X, $Xt, $Y, $B);
        $msError  = $this->mse($ssError, $dfError);

        $p_n = 3 * ($this->coefN + 1) / $this->observeN;
        $p_n = min($p_n, 0.99);

        $result = new stdclass();
        $result->observed = array();
        $result->value    = array();
        $result->fit      = array();
        $result->residual = array();
        $result->stdResid = array();
        $result->hasTest  = false;

        if(!self::is_num($ssError) or !self::is_num($msError) or $msError == 0) return $result;

        foreach($residual as $index => $e)
        {
            $setX = false;
            $setR = false;

            $stdResid = (int)$h[$index] == 1 ? '-' : $e / sqrt($msError * (1 - $h[$index]));
            if($h[$index] > $p_n) $setX = true;
            if($stdResid != '-' && abs($stdResid) > 2) $setR = true;

            if($setX || $setR)
            {
                $result->hasTest    = true;
                $result->observed[] = $index + 1;
                $result->value[]    = $y[$index];
                $result->fit[]      = $this->dataframe->roundDigit($fitY[$index], $round);
                $result->residual[] = $this->dataframe->roundDigit($e, $round);
                $stdResid = $this->dataframe->roundDigit($stdResid, $round);
                if($setR) $stdResid .= ' R';
                if($setX) $stdResid .= ' X';
                $result->stdResid[] = $stdResid;
            }
        }

        return $result;
    }

    public function rsq($ssError = null, $ssTotal = null)
    {
        if(empty($ssError)) $ssError = $this->sse();
        if(empty($ssTotal)) $ssTotal = $this->sst();

        if(!self::is_num($ssError) or !self::is_num($ssTotal) or $ssTotal == 0) return '-';

        return 1 - $ssError / $ssTotal;
    }

    public function rsqAdj($msError = null, $msTotal = null)
    {
        if(empty($msError)) $msError = $this->mse();
        if(empty($msTotal)) $msTotal = $this->mst();

        if(!self::is_num($msError) or !self::is_num($msTotal) or $msTotal == 0) return '-';

        return 1 - $msError / $msTotal;
    }

    public function rsqPre($H = null, $yFit = null)
    {
        if(empty($H))
        {
            $X  = $this->X();
            $Xt = $X->transpose(); // Xᵀ
            $B  = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
            $H  = $X->multiply($B)->multiply($Xt); //  X⟮XᵀX⟯⁻¹Xᵀ
        }
        $e    = $this->residual($yFit);
        $y    = $this->response;
        $mean = Describe::mean($y);

        $sum1 = 0;
        $sum2 = 0;
        foreach($y as $index => $value)
        {
            $h = $H[$index][$index];
            if($h == 1) continue;
            $sum1 += pow($e[$index] / (1 - $h), 2);
            $sum2 += pow($value - $mean, 2);
        }

        if($sum2 == 0) return '-';

        $rsqPre = 1 - ($sum1 / $sum2);
        $rsqPre = $rsqPre < 0 ? 0 : $rsqPre;

        return $rsqPre;
    }

    public function vif()
    {
        $vif = array();
        foreach($this->factors as $index => $factor)
        {
            $arr = $this->factors;
            unset($arr[$index]);

            $rsq = $this->vifFactor($arr, $this->response);
            $vif[] = self::is_num($rsq) ? 1 / (1 - $rsq) : '-';
        }

        return $vif;
    }

    public function vifFactor($factors, $response)
    {
        $X = $this->X($factors);
        $Xt = $X->transpose(); // Xᵀ
        $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        $Y = RegressionMatrix::yMatrix($response);
        $coef = $this->coef($X, $Xt, $XtXinv, $Y);
        $B = RegressionMatrix::coefMatrix($coef);

        $ssError = $this->sse($X, $Xt, $Y, $B);
        $ssTotal = $this->sst($Y);

        if(!self::is_num($ssError) or !self::is_num($ssTotal) or $ssTotal === 0) return '-';

        return 1 - $ssError / $ssTotal;
    }

    public function pFactor($freedom = null, $t = null)
    {
        if(empty($freedom)) $freedom = $this->ssedf();
        if(empty($t)) $t = $this->tFactor();

        $studentT = new StudentT($freedom);

        $p = array();
        foreach($t as $value)
        {
            if(!self::is_num($value))
            {
                $p[] = '-';
                continue;
            }
            $pi = 2 * (1 - $studentT->cdf($value));
            $p[] = $pi >= 1 ? 2 - $pi : $pi;
        }

        return $p;
    }

    public function tFactor($coef = null, $se = null)
    {
        if(empty($coef)) $coef = $this->coef();
        if(empty($se)) $se = $this->se();

        $t = array();
        foreach($coef as $index => $c)
        {
            if(!self::is_num($c) or !self::is_num($se[$index]) or $se[$index] == 0)
            {
                $t[] = '-';
                continue;
            }
            $t[] = $c / $se[$index];
        }

        return $t;
    }

    public function fr($msr = null, $mse = null)
    {
        if(empty($msr)) $msr = $this->msr();
        if(empty($mse)) $mse = $this->mse();

        if(!self::is_num($msr) or !self::is_num($mse) or $mse == 0) return '-';

        return $msr / $mse;
    }

    public function pr($ssrdf = null, $ssedf = null, $fr = null)
    {
        if(empty($ssrdf)) $ssrdf = $this->ssrdf();
        if(empty($ssedf)) $ssedf = $this->ssedf();
        if(empty($fr)) $fr = $this->fr();

        if(!self::is_num($ssedf) or !self::is_num($ssrdf) or !self::is_num($fr)) return '-';

        $fdist = new F($ssrdf, $ssedf);
        $cdf   = $fdist->cdf($fr);

        if(!self::is_num($cdf)) return '-';
        return 1 - $cdf;
    }

    public function ssr($X = null, $Xt = null, $Y = null, $B = null)
    {
        if(empty($X)) $X = $this->X();
        if(empty($Xt)) $Xt = $X->transpose(); // Xᵀ
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->response);
        if(empty($B))
        {
            $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
            $coef = $this->coef($Xt, $XtXinv, $Y);
            $B = RegressionMatrix::coefMatrix($coef);
        }

        $J = RegressionMatrix::oneMatrix($this->observeN, $this->observeN);

        $part1 = $B->transpose()->multiply($Xt)->multiply($Y);
        $part2 = $Y->transpose()->multiply($J)->multiply($Y)->scalarMultiply(1 /  $this->observeN);

        $ss = $part1->subtract($part2);

        return $ss->get(0, 0);
    }

    public function sse($X = null, $Xt = null, $Y = null, $B = null)
    {
        if(empty($X)) $X = $this->X();
        if(empty($Xt)) $Xt = $X->transpose(); // Xᵀ
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->response);
        if(empty($B))
        {
            $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
            $coef = $this->coef($Xt, $XtXinv, $Y);
            $B = RegressionMatrix::coefMatrix($coef);
        }

        $part1 = $Y->transpose()->multiply($Y);
        $part2 = $B->transpose()->multiply($Xt)->multiply($Y);

        $ss = $part1->subtract($part2);

        return $ss->get(0, 0);
    }

    public function sst($Y = null)
    {
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->response);

        $J = RegressionMatrix::oneMatrix($this->observeN, $this->observeN);

        $part1 = $Y->transpose()->multiply($Y);
        $part2 = $Y->transpose()->multiply($J)->multiply($Y)->scalarMultiply(1 /  $this->observeN);

        $ss = $part1->subtract($part2);

        return $ss->get(0, 0);
    }

    public function mst($ssTotal = null, $dfTotal = null)
    {
        if(empty($ssTotal)) $ssTotal = $this->sst();
        if(empty($dfTotal)) $dfTotal = $this->sstdf();

        if(!self::is_num($ssTotal) or !self::is_num($dfTotal) or $dfTotal === 0) return '-';

        return $ssTotal / $dfTotal;
    }

    public function ssrdf()
    {
        $df = 0;
        if($this->hasFactor) $df += array_sum($this->freedoms);
        if($this->hasCovariate) $df += count($this->covariate);

        return $df;
    }

    public function ssedf($ssrdf = null)
    {
        if(empty($ssrdf)) $ssrdf = $this->ssrdf();

        return $this->observeN - $ssrdf - 1;
    }

    public function sstdf()
    {
        return $this->observeN - 1;
    }

    public function pedf($same = null)
    {
        if(empty($same)) $same = $this->sameItem();

        $k = count($same);
        $s = 0;
        foreach($same as $values)
        {
            $s += count($values);
        }
        $m = $this->observeN - $s + $k;

        return $this->observeN - $m;
    }

    public function lofdf($ssedf = null, $pedf = null)
    {
        if(empty($ssedf)) $ssedf = $this->ssedf();
        if(empty($pedf))  $pedf  = $this->pedf();

        if(!self::is_num($ssedf) || !self::is_num($pedf)) return '-';

        return $ssedf - $pedf;
    }

    public function sameItem()
    {
        $row = array_fill(0, $this->observeN, '');

        for($i = 0; $i < $this->observeN; $i++)
        {
            if($this->hasFactor)
            {
                foreach($this->factors as $index => $factor)
                {
                    $row[$i] .= $factor[$i];
                }
            }

            if($this->hasCovariate)
            {
                foreach($this->covariate as $value)
                {
                   $row[$i] .= $value[$i];
                }
            }
        }

        $same = array();
        foreach($row as $index => $key)
        {
            if(!isset($same[$key])) $same[$key] = array();
            $same[$key][] = $this->response[$index];
        }

        foreach($same as $index => $value)
        {
            if(count($value) == 1) unset($same[$index]);
        }

        return $same;
    }

    public function mse($sse = null, $ssedf = null)
    {
        if(empty($sse)) $sse = $this->sse();
        if(empty($ssedf)) $ssedf = $this->ssedf();

        if(!self::is_num($sse) || !self::is_num($ssedf) || $ssedf == 0) return '-';
        return $sse / $ssedf;
    }

    public function se($XtXinv = null, $mse = null)
    {
        if(empty($XtXinv))
        {
            $X      = $this->X();
            $Xt     = $X->transpose(); // Xᵀ
            $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        }
        if(empty($mse)) $mse = $this->mse();

        if(!self::is_num($mse)) return array_fill(0, $this->coefN + 1, '-');

        $se = $XtXinv->scalarMultiply($mse);

        $diagonal = $se->getDiagonalElements();

        $sqrtDiago = array();
        foreach($diagonal as $value)
        {
            $sqrtDiago[] = sqrt($value);
        }

        return $sqrtDiago;
    }

    public function msr($ssr = null, $ssrdf = null)
    {
        if(empty($ssr)) $ssr = $this->ssr();
        if(empty($ssrdf)) $ssrdf = $this->ssrdf();

        if(!self::is_num($ssr) or !self::is_num($ssrdf) or $ssrdf == 0) return '-';

        return $ssr / $ssrdf;
    }

    public function sslof($sse = null, $sspe = null)
    {
        if(empty($sse)) $sse = $this->sse();
        if(empty($sspe)) $sspe = $this->sspe();

        if(!self::is_num($sse) or !self::is_num($sspe)) return '-';

        return $sse - $sspe;
    }

    public function mslof($sslof = null, $lofdf = null)
    {
        if(empty($sslof)) $sslof = $this->sslof();
        if(empty($lofdf)) $lofdf = $this->lofdf();

        if(!self::is_num($sslof) or !self::is_num($lofdf) or $lofdf === 0) return '-';

        return $sslof / $lofdf;
    }

    public function flof($mslof = null, $mspe = null)
    {
        if(empty($mslof)) $mslof = $this->mslof();
        if(empty($mspe)) $mspe = $this->mspe();

        if(!self::is_num($mslof) or !self::is_num($mspe) or $mspe === 0) return '-';

        return $mslof / $mspe;
    }

    public function plof($lofdf = null, $pedf = null, $flof = null)
    {
        if(empty($lofdf)) $lofdf = $this->lofdf();
        if(empty($pedf)) $pedf = $this->pedf();
        if(empty($flof)) $flof = $this->flof();

        if(!self::is_num($lofdf) or !self::is_num($pedf) or !self::is_num($flof)) return '-';

        $fdist = new F($lofdf, $pedf);
        return 1 - $fdist->cdf($flof);
    }

    public function sspe($same = null)
    {
        if(empty($same)) $same = $this->sameItem();
        $sum = 0;
        foreach($same as $values)
        {
            $mean = Describe::mean($values);
            foreach($values as $value)
            {
                $sum += ($value - $mean) ** 2;
            }
        }

        return $sum;
    }

    public function mspe($sspe = null, $pedf = null)
    {
        if(empty($sspe)) $sspe = $this->sspe();
        if(empty($pedf)) $pedf = $this->pedf();

        if(!self::is_num($sspe) or !self::is_num($pedf) or $pedf === 0) return '-';

        return $sspe / $pedf;
    }

    public function coef($Xt = null, $XtXinv = null, $Y = null)
    {
        if(empty($Xt))
        {
            $X  = $this->X();
            $Xt = $X->transpose(); // Xᵀ
        }
        if(empty($XtXinv)) $XtXinv = $Xt->multiply($X)->inverse();// ⟮XᵀX⟯⁻¹
        if(empty($Y)) $Y = RegressionMatrix::yMatrix($this->response);

        return $XtXinv->multiply($Xt)->multiply($Y)->getColumn(0);
    }

    public function residual($yFit = null)
    {
        $e    = array();
        $y    = $this->response;
        if(empty($yFit)) $yFit = $this->fitY();

        foreach($yFit as $index => $value)
        {
           $e[] = $y[$index] - $value;
        }

        return $e;
    }

    public function fitY($coef = null)
    {
        if(empty($coef)) $coef = $this->coef();

        list($common, $factors, $covariates) = $this->fullCoef($coef);

        $yFit = array();
        for($index = 0; $index < $this->observeN; $index++)
        {
            $y = $common;
            if(!empty($factors))
            {
                foreach($this->factors as $level => $values)
                {
                    $key = $values[$index];
                    $y += $factors[$level][$key];
                }
            }

            if(!empty($covariates))
            {
                foreach($this->covariate as $name => $values)
                {
                    $value = $values[$index];
                    $y += $covariates[$name] * $value;
                }
            }

            $yFit[] = $y;
        }

        return $yFit;
    }

    public function fullCoef($coef = null)
    {
        if(empty($coef)) $coef = $this->coef();

        $common     = $coef[0];
        $factors    = array();
        $covariates = array();

        $count = 1;
        if($this->hasCovariate)
        {
            foreach($this->covariate as $name => $item)
            {
                $covariates[$name] = $coef[$count];
                $count += 1;
            }
        }

        if($this->hasFactor)
        {
            foreach($this->factorLevel as $level => $fLv)
            {
                $factors[$level] = array();
                $items = $fLv->items;
                $last = array_pop($items);
                foreach($items as $name)
                {
                    $factors[$level][$name] = $coef[$count];
                    $count += 1;
                }
                $factors[$level][$last] = -1 * array_sum($factors[$level]);
            }
        }

        return array($common, $factors, $covariates);
    }

    public function X($factors = null, $covariate = null)
    {
        if($factors === null) $factors = $this->factors;
        if($covariate === null) $covariate = $this->covariate;

        $hasFactor = ($this->hasFactor and !empty($factors));
        $hasCovariate = ($this->hasCovariate and !empty($covariate));

        if($hasFactor and !$hasCovariate) $X = RegressionMatrix::glmFactorDM($factors);
        if(!$hasFactor and $hasCovariate) $X = RegressionMatrix::glmCovariateDM($covariate);
        if($hasFactor and $hasCovariate)  $X = RegressionMatrix::glmFCDM($factors, $covariate);
        return $X;
    }

    public static function factorValues($factors)
    {
        $data = array();
        foreach($factors as $factor)
        {
            $data[$factor->name] = $factor->data;
        }
        return $data;
    }

    public static function covariateValues($covariates)
    {
        $data = array();
        foreach($covariates as $covariate)
        {
            $data[$covariate->name] = $covariate->data;
        }
        return $data;
    }
}
