Statistical analysis of growth curves with Mathematica

Martin Ricker

Instituto de Biología, Universidad Nacional Autónoma de México, México D.F., Mexico
<mricker@ib.unam.mx>, <martin_tuxtlas@yahoo.com.mx>

PL model:     Ln[(dQ/dt) /Q] = a + b·Q   with (dQ/dt) = annual increment, (dQ/dt) /Q = relative growth

               =>  t = (Subscript[t, 1])+ (e^(-a))·(Ei[-b·Q] - Ei[-b·(Subscript[Q, 1])])   with Ei(x) = ∫ ((e^x)/x)dx
               
                   a, b, (Subscript[t, 1]) or (Subscript[Q, 1]) = regression coefficients
                   
                   a = initial logarithmic relative growth
                   b = -1 /(Turning point quantity of sigmoid curve)
                   (Subscript[Q, 1]) at (Subscript[t, 1]) = one known time-quantity point
                   
                   with b = 0 exponential growth of Q(t)
                   with b < 0 over-exponential growth
                   with a = large negative number, growth becomes almost linear

Two segments with one knot: 7.8 cm at 73.7 years:

t1Seg1 = Mean[AgesSeg1];
SolSeg1 = NonlinearModelFit[InputDataSeg1, t1Seg1 + Exp[-aSeg1] * (ExpIntegralEi[-bSeg1*Q2Seg1] - ExpIntegralEi[-bSeg1*Q1Seg1]), {Q1Seg1, aSeg1, bSeg1}, Q2Seg1];

Q1Seg2 = Mean[QuantitiesSeg2];
SolSeg2 = NonlinearModelFit[InputDataSeg2, t1Seg2 + (t1Seg1 - t1Seg2 + Exp[-aSeg1] * ExpIntegralEi[-bSeg1 * OptKnotQ] - Exp[-aSeg1] * ExpIntegralEi[-bSeg1*Q1Seg1]) /(ExpIntegralEi[-bSeg2*OptKnotQ] - ExpIntegralEi[-bSeg2*Q1Seg2]) * (ExpIntegralEi[-bSeg2*Q2Seg2] - ExpIntegralEi[-bSeg2*Q1Seg2]), {t1Seg2, bSeg2}, Q2Seg2];

Ten segments with nine knots:
3.1 cm at 14.8 years, 4.7 cm at 31.4 years, 5.0 cm at 39.4 years, 7.2 cm at 56.6 years, 7.6 cm at 73.3 years
8.5 cm at 82.8 years, 9.8 cm at 102.3 years, 10.2 cm at 106.9 years, 11.0 cm at 114.0 years

10 volcano sites in Mexico ≈3100 m above sea level (data from Dr. Víctor Peña, UNAM, doctoral thesis 2013)

235 Pinus montezumae trees: Radius: 9.5 cm to 43.8 cm, Annual radius increment: 0.04 to 0.86 cm

192 trees with annual growth from 1944 to 2004 (60 years):

FREE KNOT
Optimal knot age per tree: 4.2 to 53.9 years, interquartile range: 24.1 to 38.6 years
Overall (R^2) per tree: 0.345 to 0.997, interquartile range: 0.888 to 0.976

KNOT FIXED AT 30 YEARS
Overall (R^2) per tree: 0.285 to 0.997, interquartile range: 0.861 to 0.968

235 trees with annual growth from 1974 to 2004 (30 years):

(Subscript[Q, 1]) => 4.6 to 37.5 cm
a   => -15.7 to 95.0 (extremely wide range)
b   => -2.65 to 0.61 (cm^(-1)) ⇒ TPQ = 0.38 cm to over-exponential growth
(R^2) => 0.84 to 0.9998

PROBLEM: Correlations among regression coefficients close to 1

SOLUTION: Fix b at the the inverse of half the quantity range ⇒ If range 10 to 30 cm, b = -1/20 = -0.05

(Subscript[Q, 1]) => 8.3 to 37.5 cm
a   => -5.2 to -1.5
b   => fixed at -0.12 to -0.027 (cm^(-1)), depending on tree ⇒ TPQ = 8.3 to 37.3 cm
(R^2) => 0.956 to 0.9998
Correlation between (Subscript[Q, 1]) and a ⇒ -4·(10^(-12)) to 1·(10^(-11))

Generalized Multivariate Analysis of Variance (GMANOVA) or “classical growth curve model” established by Potthoff & Roy (1964) for short time series data (5 to 50 time point reasonable)

Y = T.A.G

A polynomial function of any degree is used, e.g., Y = (Subscript[c, 0]) + (Subscript[c, 1])·T + (Subscript[c, 2])·(T^2)

HERE: (Ei[-b·(Subscript[Q, i])] - Ei[-b·(Subscript[Q, 1])]) = A·((Subscript[t, i]) - (Subscript[t, 1]))  with  (e^a) = A = initial relative growth

corresponding to  Y = (Subscript[c, 1])·T

Y = T.A.G  with  n = 235 trees, k = 10 volcano sites, p = 31 time points

YMatrix=(({{-0.2274486713590229`, -0.208544469135713`, -0.20839648084115114`, -0.375822335255392`, -0.4292098515096592`...}, {-0.21816630061413167`, -0.17976177844228647`, -0.2032168612887304`, -0.3595535338166649`, -0.4165661036760462`...}, {-0.21309753849134316`, -0.1554890831061393`, -0.19382594005021359`, -0.34366052487155074`, -0.3934900786831539`...}, {-0.1974181091328493`, -0.13204968320632648`, -0.18997663228044903`, -0.32480103989648623`, -0.35980058327743314`...}, {-0.1817447911834975`, -0.12824911153619944`, -0.1778273173732663`, -0.2877249540627751`, -0.35285523390609397`...}, {-0.17477252501747342`, -0.12141249328523585`, -0.16485174636101974`, -0.25624112103460184`, -0.33256950575851696`...}, {-0.14276162155957328`, -0.11451931361911072`, -0.15041689458065144`, -0.2448495090034568`, -0.3075147714695534`...}, {-0.12650624032708668`, -0.10207128235658769`, -0.13304607878571373`, -0.21767083846752855`, -0.2944151619289037`...}, {-0.08762968361650314`, -0.09460723537206883`, -0.08691229967615977`, -0.15218032328429532`, -0.2725767083456325`...}, {-0.07509758113866272`, -0.08187589709778598`, -0.0684866789244507`, -0.10668583958057032`, -0.1998468593749536`...}, {-0.060661344026970765`, -0.06989041090885895`, -0.057173486384355376`, -0.10080816935460035`, -0.12173597685174986`...}, {-0.055214456760733555`, -0.05552077107469611`, -0.03570894835846827`, -0.06800111707756806`, -0.05783280723944917`...}, {-0.04693443980973422`, -0.04017924046387722`, -0.02677716179855838`, -0.05183309254755142`, 0.011218597654839568`...}, {-0.04063297657753395`, -0.018942912975550286`, -0.007242039009841994`, -0.045894574811154776`, 0.024435078613512573`...}, {-0.01730826780321504`, -0.010700426533420071`, 0.007578937847400713`, -0.037610446153710564`, 0.034457787205271595`...}, {0.03266681038714325`, 0.0009001486440087803`, 0.03205625665601208`, -0.02402227774140986`, 0.04836561002768636`...}, {0.05563240449040352`, 0.027098028280035047`, 0.045207629992837006`, -0.0042846928745012836`, 0.061679999754040926`...}, {0.0817565217209626`, 0.03980695721853822`, 0.053674623758378504`, 0.035806785929724416`, 0.06941890491051184`...}, {0.08901616347388575`, 0.05495991911227427`, 0.058224763854738404`, 0.0745355520822677`, 0.09652688519010466`...}, {0.09368523040187138`, 0.057727353901323664`, 0.06687359911912294`, 0.0967828230566814`, 0.11292242861670143`...}, {0.10326888036770332`, 0.06400683376819605`, 0.07330800180810937`, 0.11393852933607529`, 0.1315808880355056`...}, {0.11300077339101988`, 0.0707066540868424`, 0.08581654546215911`, 0.1432696709243546`, 0.17472694911426423`...}, {0.11750050076451513`, 0.0766815479661771`, 0.10378426078323799`, 0.20029254692005494`, 0.2174166142780538`...}, {0.12131577644240221`, 0.09554444603341361`, 0.12270696550942017`, 0.2674926471534713`, 0.25518813478550184`...}, {0.12654999717009763`, 0.11238777393213084`, 0.12812103952322107`, 0.2702967232775084`, 0.2768157965718414`...}, {0.13315467995442298`, 0.12919745630570922`, 0.13462165446624152`, 0.2752470186862894`, 0.29152722323849867`...}, {0.14670930333561016`, 0.13620933974182226`, 0.14119850304723092`, 0.2801981888570182`, 0.3054271353593534`...}, {0.14908298918815333`, 0.1471015979919863`, 0.1477762112291432`, 0.29176148298774485`, 0.32469794036489175`...}, {0.15923987296143083`, 0.15159295846722642`, 0.15809718255915217`, 0.3020961585626172`, 0.3507428177530545`...}, {0.16509034049547577`, 0.16550933731528028`, 0.16513004923943186`, 0.31247750844568034`, 0.37014584019072894`...}, {0.17772430287322183`, 0.1843837563489008`, 0.17968234085721124`, 0.33744875361434157`, 0.38111901035678697`...}}))               TMatrix=(({{-15.`}, {-14.`}, {-13.`}, {-12.`}, {-11.`}, {-10.`}, {-9.`}, {-8.`}, {-7.`}, {-6.`}, {-5.`}, {-4.`}, {-3.`}, {-2.`}, {-1.`}, {0.`}, {1.`}, {2.`}, {3.`}, {4.`}, {5.`}, {6.`}, {7.`}, {8.`}, {9.`}, {10.`}, {11.`}, {12.`}, {13.`}, {14.`}, {15.`}}))

GMatrix=(({{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...}}))



Solution vector of initial relative growth rates:

AMatrix = Inverse[Transpose[TMatrix].TMatrix].
Transpose[TMatrix].QMatrix.Transpose[GMatrix].Inverse[GMatrix.Transpose[GMatrix]]

{0.017367846364308794`,0.01738905276826874`,0.015280296805135025`,0.013304793119881614`,0.01598096460928944`,0.013870230834962301`,0.015358359230102871`,0.0523492816325794`,0.021664069978226987`,0.021831280443677437`}

Sum-of-squares matrix:
SMatrix = MatrixForm[QMatrix.(IdentityMatrix[n] - Transpose[GMatrix].Inverse[GMatrix.Transpose[GMatrix]].GMatrix).Transpose[QMatrix]]

(({{7.747840834394986`, 6.967070238983878`, 6.153372772872473`, 5.383933922819434`, 4.650539941028728`...}, {6.967070238983878`, 6.297677708334401`, 5.575192200048873`, 4.890014849956493`, 4.238498319785522`...}, {6.153372772872468`, 5.575192200048875`, 4.971179810400497`, 4.385697177610464`, 3.8102002470669376`...}, {5.383933922819438`, 4.8900148499564935`, 4.385697177610462`, 3.899056567565258`, 3.401215867170604`...}, {4.650539941028726`, 4.238498319785522`, 3.8102002470669394`, 3.401215867170603`, 2.999546639882121`...}, {4.150777219599045`, 3.788061963899026`, 3.4112583998440744`, 3.05202812985324`, 2.7028810293898538`...}, {3.4572913935289935`, 3.164231880567711`, 2.8565748347061692`, 2.5645340139455275`, 2.28118870702223`...}, {2.895727232271776`, 2.6561175265532926`, 2.4031282616947665`, 2.1625359505992527`, 1.9314691381935118`...}, {2.5046273693979857`, 2.2940739897958733`, 2.0766252926423614`, 1.8689342488758185`, 1.6718401634848075`...}, {2.004526235344442`, 1.8348355278587827`, 1.6545306494209409`, 1.4836823112020219`, 1.326755472481372`...}, {1.4946574500748269`, 1.370795532250152`, 1.2388132312905522`, 1.1138686598421637`, 1.0017508826533044`...}, {1.079363307819414`, 0.9920197095560038`, 0.8981479854341103`, 0.8099466395577661`, 0.7326210423585443`...}, {0.6892523073196045`, 0.6375224924098446`, 0.5810099298453545`, 0.5299803927905467`, 0.48619513167964495`...}, {0.3558336291360276`, 0.330010900579125`, 0.2961624641229762`, 0.26540716098201644`, 0.24872160474768346`...}, {-0.007728282766628337`, 0.0020132404001184934`, 0.0008729295972831497`, 0.0022767694043443164`, 0.01736759418791569`...}, {-0.43667564235891687`, -0.38611944585711366`, -0.34447616485350185`, -0.3019836207514729`, -0.24704441374951838`...}, {-0.7432840073555043`, -0.6700688673349405`, -0.5958051687310889`, -0.5240244243094222`, -0.4445564698647432`...}, {-1.0828423776800835`, -0.9799273875104886`, -0.8721164907431498`, -0.7697446035171912`, -0.6625086334373419`...}, {-1.6758623821835126`, -1.5211270568593718`, -1.3542752114190193`, -1.1951035304190023`, -1.0375455249424812`...}, {-2.11994069147913`, -1.9276088647901335`, -1.7147355334537826`, -1.5154658540707564`, -1.3240721349396005`...}, {-2.479705620883024`, -2.2548506550412593`, -2.006351056998529`, -1.7717970813189763`, -1.5513864389563292`...}, {-2.7649449918693803`, -2.5149957592999908`, -2.2426553705814554`, -1.984104286995833`, -1.741010680589526`...}, {-3.1164816163219244`, -2.8338804399675697`, -2.5292861150507218`, -2.2404662300547606`, -1.9696663540422832`...}, {-3.5011806903030274`, -3.1836555787841645`, -2.8456763826808618`, -2.524916148253854`, -2.2256441076635727`...}, {-3.729531253210182`, -3.394247215386329`, -3.0390817384441617`, -2.701758073922387`, -2.3862637711552157`...}, {-3.9545608379963304`, -3.598412570778048`, -3.2244168628427357`, -2.868196563231714`, -2.5319518607356564`...}, {-4.2342216301492215`, -3.850613934819402`, -3.452469862651523`, -3.0708866055019017`, -2.708622256798798`...}, {-4.486887521579205`, -4.08140670039364`, -3.659747782995569`, -3.2573681945342288`, -2.8742284330628096`...}, {-4.8070733309492555`, -4.376715865691087`, -3.9325193358000425`, -3.509003037281681`, -3.0995548045810724`...}, {-5.0408002303205475`, -4.591626083986798`, -4.121331007705368`, -3.6735705477473237`, -3.2461799857257625`...}, {-5.353092747223862`, -4.872879654512378`, -4.377822101675613`, -3.904723860289369`, -3.4505559109109325`...}}))

Standard errors:

DAMatrix = n*(n-k-1) /(n-k-p) /(n-k-p+1)*KroneckerProduct[Inverse[GMatrix.Transpose[GMatrix]], Inverse[Transpose[TMatrix].PseudoInverse[SMatrix].TMatrix]];

Diagonal[DAMatrix]

{0.0004291118884887058`,0.0005363898606108823`,0.00039732582267472757`,0.0004469915505090685`,0.0004469915505090685`,0.0005646209059061918`,0.0007151864808145098`,0.0004291118884887058`,0.000487627146009893`,0.0005363898606108823`}

Outliers: points in a curve that are more than 4 standard deviations off the regression line; the whole curve is marked and eliminated in that case:

Autocorrelation (lag = 1) of median standardized residuals of 31 time points: 0.761

Revelation of non-obvious outliers:

Overall significance test in GMANOVA model for significant differences among means: (Subscript[F, S])= 192.3  ⇒  P = 3·(10^(-30))

GMANOVA pairwise comparisons: 22 (49%) of 45 comparisons are significant

Statistical comparison intervals

Advantages of the approach of combining the PL model with the GMANOVA model:

- Flexible and accurate growth curves (although no direct modeling of negative or cycling growth modeling)
- Autocorrelation and non-constant variance taken into account
- Standard errors very low: here around 18 times lower than without GMANOVA approach

Future research:

- Taking several independent (explaining) variables into account
- Modeling more than one segment simultaneously

Thanks:

- to my collaborators Dietrich von Rosen (Swedish University of Agricultural Sciences in Uppsala) and Víctor Peña (UNAM in Mexico City)
- Stephen Wolfram and collaborators for having developed Mathematica!