CLS
PRINT "PROGRAM FOR MULTIPLE PVR (SIMULATED DATA)"
PRINT "Version 4.0"
PRINT "by Alexandre S.G.Coelho & J. A. F. Diniz-Filho"
PRINT
PRINT
DEFINT I-K
DEFDBL B, M, R, T, X-Y
INPUT "Filename with EIGENVECTORS (X matrix)"; ax$
OPEN ax$ FOR INPUT AS #1
PRINT
INPUT "Filename with multiple DATA SETS (with two traits)"; ay$
PRINT
INPUT "OUTPUT file"; outcor$
PRINT
PRINT "******************************************************************"
PRINT
INPUT "Number of simulations"; f
PRINT
PRINT
INPUT "Number of species"; n
INPUT "Number of eigenvectors (+ constant)"; p

DIM x(n, p), xx(p, p), xxi(p, p), y1(n), y2(n), xy1(p), xy2(p), b1(p), b2(p), ye1(n), ye2(n), e1(n), e2(n), m1(n), m2(n)
PRINT
INPUT "Parameter (rho) for mean squared-deviation (MDS)"; rho
INPUT "Critical F-value at 5% level"; fcrit
LET gg$ = "fc.out"
OPEN gg$ FOR OUTPUT AS #7
PRINT #7, fcrit

LET nn$ = "n.out"
OPEN nn$ FOR OUTPUT AS #6
PRINT #6, f

FOR i = 1 TO n
FOR j = 1 TO p
INPUT #1, x(i, j)
NEXT j
NEXT i
CLOSE

FOR i = 1 TO p
FOR j = 1 TO p
FOR k = 1 TO n
xx(i, j) = xx(i, j) + x(k, i) * x(k, j)
NEXT k
NEXT j
NEXT i

FOR i = 1 TO p
FOR j = 1 TO p
IF i = j THEN xxi(i, j) = 1 ELSE xxi(i, j) = 0
NEXT j
NEXT i

FOR D = 1 TO p
m = xx(D, D)
FOR j = 1 TO p
xxi(D, j) = xxi(D, j) / m
xx(D, j) = xx(D, j) / m
NEXT j

FOR i = D + 1 TO p
m = xx(i, D)
FOR j = 1 TO p
xxi(i, j) = xxi(i, j) - m * xxi(D, j)
xx(i, j) = xx(i, j) - m * xx(D, j)
NEXT j
NEXT i
NEXT D

FOR D = p TO 1 STEP -1
FOR i = D - 1 TO 1 STEP -1
m = xx(i, D)
FOR j = 1 TO p
xxi(i, j) = xxi(i, j) - m * xxi(D, j)
xx(i, j) = xx(i, j) - m * xx(D, j)
NEXT j
NEXT i
NEXT D


LET outt$ = "t.dat"
LET outr$ = "r.dat"
LET outrq$ = "r2.dat"

OPEN ay$ FOR INPUT AS #1
OPEN outcor$ FOR OUTPUT AS #2
OPEN outr$ FOR OUTPUT AS #3
OPEN outt$ FOR OUTPUT AS #4
OPEN outrq$ FOR OUTPUT AS #5

FOR s1 = 1 TO f

CLS
LOCATE 10, 25
REM PRINT LEFT$(ay$, 7); "  "; "simulation :"; s1
PRINT "Analyzing simulation number..."; s1

FOR s2 = 1 TO n
INPUT #1, y1(s2), y2(s2)
ssy1 = ssy1 + y1(s2)
ssy2 = ssy2 + y2(s2)
NEXT s2
xby1 = ssy1 / n
xby2 = ssy2 / n


FOR i = 1 TO p
FOR k = 1 TO n
xy1(i) = xy1(i) + x(k, i) * y1(k)
xy2(i) = xy2(i) + x(k, i) * y2(k)
NEXT k
NEXT i

FOR i = 1 TO p
FOR k = 1 TO p
b1(i) = b1(i) + xxi(k, i) * xy1(k)
b2(i) = b2(i) + xxi(k, i) * xy2(k)
NEXT k
NEXT i


FOR i = 1 TO n
FOR k = 1 TO p
ye1(i) = ye1(i) + x(i, k) * b1(k)
ye2(i) = ye2(i) + x(i, k) * b2(k)

NEXT k
NEXT i


FOR i = 1 TO n
m1(i) = ye1(i) - xby1
m2(i) = ye2(i) - xby2
e1(i) = y1(i) - ye1(i)
e2(i) = y2(i) - ye2(i)
NEXT i


FOR i = 1 TO n
sm1 = sm1 + m1(i) ^ 2
sm2 = sm2 + m2(i) ^ 2
te1 = te1 + e1(i)
te2 = te2 + e2(i)
te1q = te1q + e1(i) ^ 2
te2q = te2q + e2(i) ^ 2
te12 = te12 + e1(i) * e2(i)
NEXT i
st1 = sm1 + te1q
st2 = sm2 + te2q

rq1 = sm1 / st1
rq2 = sm2 / st2
f1 = (sm1 / (p - 1)) / (te1q / (n - p))
f2 = (sm2 / (p - 1)) / (te2q / (n - p))
r = (te12 - te1 * te2 / n) / SQR((te1q - te1 * te1 / n) * (te2q - te2 * te2 / n))
msd = (r - rho) * (r - rho)
tp = r * (SQR((n - 2 - (p - 1)) / (1 - (r * r))))

PRINT #2, USING " ###.### "; rq1; rq2; f1; f2; r; tp; msd
PRINT #3, USING " ###.####"; r
PRINT #4, USING " ###.####"; tp
PRINT #5, USING " ###.####"; rq1; f1

FOR i = 1 TO p
xy1(i) = 0
xy2(i) = 0
b1(i) = 0
b2(i) = 0
NEXT i
FOR i = 1 TO n
ye1(i) = 0
m1(i) = 0
m2(i) = 0
ye2(i) = 0
NEXT i
te1 = 0
te2 = 0
te1q = 0
te2q = 0
te12 = 0

ssy1 = 0
ssy2 = 0
xby1 = 0
xby2 = 0
sm1 = 0
sm2 = 0
st1 = 0
st2 = 0
rq1 = 0
rq2 = 0
f1 = 0
f2 = 0
tp = 0
msd = 0

NEXT s1

CLOSE 1
CLOSE 2
CLOSE #3
CLOSE #4
CLOSE #5
CLOSE #6
CLOSE #7
CLS
PRINT
PRINT "OUTPUT FILE: "; outcor$
PRINT
PRINT "The 7 columns in the output file are:"
PRINT
PRINT "1. Squared multiple-R for trait 1"
PRINT "2. Squared multiple-R for trait 2"
PRINT "3. F-value for trait 1"
PRINT "4. F-value for trait 2"
PRINT "5. Partial correlation between traits"
PRINT "6. t-value for testing partial correlation"
PRINT "7. Mean squared-deviation (MSD) from rho"
PRINT
PRINT
PRINT

CLEAR
PRINT "PHYLOGENETIC INERTIA"
nf = 0
PRINT
LET rq2$ = "r2.dat"
LET fc$ = "fc.out"
LET nn$ = "n.out"

OPEN rq2$ FOR INPUT AS #1
OPEN fc$ FOR INPUT AS #2
OPEN nn$ FOR INPUT AS #3

INPUT #2, fcrit
INPUT #3, n
DIM rq(n), f(n)

FOR i = 1 TO n
        INPUT #1, rq(i), f(i)
        IF f(i) > fcrit THEN LET nf = nf + 1
        sr = sr + rq(i)
NEXT i

LET med = sr / n
stq2 = 0
FOR i = 1 TO n
        LET stq2 = stq2 + ((rq(i) - med) * (rq(i) - med))
NEXT i
LET var = stq2 / n
LET dpr2 = SQR(var)

PRINT "Mean Squared multiple-R (phylogenetic inertia):", USING "##.###"; med
PRINT "(ñ "; dpr2; ")"
PRINT
PRINT "Proportion of significant F-values (at 5% level) in the simulations:"
PRINT USING " #.###"; nf / n
PRINT
PRINT "........................................................"

CLOSE #1
CLOSE #2
CLOSE #3

INPUT "Do you want to estimate Type I error rates? [S/N]"; b$
IF b$ = "N" OR b$ = "n" THEN RUN


CLS
PRINT
CLEAR
PRINT "TYPE I ERROR RATE ANALYSES"
PRINT
INPUT "Number of classes in the distribution"; k
PRINT
PRINT "Filename of the vector with critical t-values"
INPUT vec$
OPEN vec$ FOR INPUT AS #1
DIM lv1(k), lv2(k)
FOR j = 1 TO k
        INPUT #1, lv1(j), lv2(j)
NEXT j

LET nn$ = "n.out"
OPEN nn$ FOR INPUT AS #6
INPUT #6, n

REM "Name of the vector with distribution of statistics"
dist$ = "t.dat"
OPEN dist$ FOR INPUT AS #2
DIM nu(k), cn(k), ch(k), sch(j), st(n), st2(n), ecrt(k), r(n)
rdist$ = "r.dat"
OPEN rdist$ FOR INPUT AS #3

FOR j = 1 TO k
        nu(k) = 0
NEXT j
PRINT
INPUT "Critical t-value at 5% level"; rc5
PRINT
INPUT "Parameter rho for MSD"; pr

FOR i = 1 TO n
        INPUT #3, r(i)
        INPUT #2, st(i)
        LET st(i) = ABS(st(i))
NEXT i

af5 = 0
st = 0
FOR i = 1 TO n
        FOR j = 1 TO k
                IF st(i) < lv1(j) AND st(i) > lv2(j) THEN LET nu(j) = nu(j) + 1
        NEXT j
        IF st(i) > rc5 THEN LET af5 = af5 + 1
        LET st = st + ((r(i) - pr) * (r(i) - pr))
        LET st2(i) = r(i) * r(i)

NEXT i
md = st / n
sv = 0
FOR i = 1 TO n
        sv = sv + (st2(i) - md) * (st2(i) - md)
NEXT i
var = sv / n
dp = SQR(var)
CLS
PRINT
PRINT "Results from COUNT..."
PRINT "(class   -   count)"
PRINT
FOR j = 1 TO k
        PRINT j, nu(j)
NEXT j
PRINT
PRINT
PRINT "Type I error at 5% level: ", USING "##.####"; af5 / n
PRINT
PRINT "Mean squared-deviation from rho: ", USING "##.####"; st / n
PRINT "(ñ "; dp; ")"
PRINT
PRINT
PRINT
INPUT "Do you want to calculate the Chi-square? [S/N]"; s$
IF s$ = "S" OR s$ = "s" THEN GOTO CHI:

RUN


CHI:
CLS
PRINT "Enter the name of the vector with the expected counts"
INPUT cont$
OPEN cont$ FOR INPUT AS #4

FOR j = 1 TO k
        INPUT #4, ecrt(j)
        LET sch(j) = (nu(j) - ecrt(j)) * (nu(j) - ecrt(j))
        LET ch(j) = sch(j) / ecrt(j)
        chi2 = chi2 + ch(j)
NEXT j
PRINT
CLS
PRINT "Results from CHI-SQUARE..."
PRINT
PRINT "Squared deviations from expectation, by class"
FOR j = 1 TO k
        PRINT USING " ####.###"; ch(j)
NEXT j
PRINT
PRINT
PRINT "Chi-square: "; chi2
PRINT
PRINT "(with ["; k - 1; "]"; " degrees of freedom)"
PRINT
PRINT
PRINT
PRINT "OBS: Significant Chi-square indicates deviations"
PRINT "from expected null distribution!"
PRINT
PRINT "End of processing...[return for another analyses]"
LINE INPUT c$
RUN