USING MATHEMATICA IN A BIOMEDICAL
STATISTICS CLASSROOM:
BLOOD GLUCOSE LEVEL ANALYSIS

Iveta Bryjová and Stanislav Hledík

Institute of Physics, Silesian University in Opava, Czech Republic

Abstract

This contribution demonstrates an example of the possible use of Mathematica in statistics courses for biomedical students. The text, combined seamlessly with statistical computations using built-in features of Mathematica, provides the students a step-by-step guide through pitfalls of statistical analysis using a particular case study: blood glucose measurement records of a patient suffering from type 1 (juvenile) diabetes mellitus. Both descriptive and inferential statistics are included.  Attention is particularly paid to the ability of students to select the correct hypothesis test.

Objective

Many people fear of dislike mathematics, on the other hand, statistical methods are essential for graduates, postgraduates, and researchers in virtually all sciences. Students of biomedical sciences are usually made to take some compulsory undergraduate course on statistics. However, even if they manage to pass the exams successfully, they are often unsure about applying statistical methods and have difficulty knowing which statistical test to be (and which not to be) used. Some students even consider statistical courses boring and hard to understand.

A lot of textbooks on statistics take profound mathematical approach and often go into details which usually dim the basic concepts of statistics if their audience lack a strong mathematical background. Although there are good statistical textbooks designed for life scientists (e.g., McKillup, 2006), we feel that students (in our case the biomedical students) need examples providing a step-by-step guide through pitfalls of statistical analysis. Mathematica appears to be an ideal tool to support such a goal.

This contribution is one of our teaching materials proposed to achieve the intention outlined above. It is based on subject matter close to every student of biomedicine, which serves as a basis for introducing and practicing essential concepts of statistics and their correct practical usage.

Preliminaries

A brief medical introduction

Blood glucose and insulin

Glucose is the primary source of energy for cells. It is transported from the intestines or liver to body cells via the bloodstream, and is made available for cell absorption via the hormone insulin, produced by the body primarily in the pancreas.

The blood glucose level is defined as the amount of glucose present in the blood. Normally the mammals maintain the reference range of blood glucose level in mammals between about 3.6 and 5.8 millimoles/liter (65 and 104 milligrams/deciliter), with small surges after the meals. The human body naturally tightly regulates blood glucose levels as a part of metabolic homeostasis.

Glucose level fluctuation

The blood glucose level fluctuates throughout the day—usually lowest in the morning, before the first meal of the day (termed “the fasting level”), and rises after meals for an hour or two. Standard course of blood glucose fluctuation of a healthy person during the day with three dishes (breakfast, lunch and dinner) looks like the following (plot based on paper by Daly, 1998):

Diabetes mellitus

Diabetes mellitus (often simply referred to as diabetes) is a group of metabolic diseases in which a person has persistent hyperglycemia from any of several causes (Ganong, 2005; Kahn et al., 2004), most frequently

either because the body does not produce enough insulin (presently requires the person to inject insulin), referred to as Type 1 diabetes (also IDDM = insulin-dependent diabetes mellitus, or juvenile diabetes),

or because cells fail to use insulin properly (sometimes combined with an absolute insulin deficiency), referred to as Type 2 diabetes (also NIDDM = non-insulin-dependent diabetes mellitus, or adult-onset diabetes),

or when pregnant women, who have never had diabetes before, have a high blood glucose level during pregnancy (gestational diabetes). This type may precede development of Type 2 DM.

Diabetes have been treatable since insulin became available in 1921. While gestational diabetes usually resolves after delivery, Type 1 and 2 are chronic conditions that usually cannot be cured—the ill person is permanently dependent on external source of insulin. Diabetes is a serious illness with consequences which influence the patients for all their rest of life.  Without proper treatments it can cause many complications, both acute (including hypoglycemia, diabetic ketoacidosis, or nonketotic hyperosmolar coma) and long-term ones (including cardiovascular disease, chronic renal failure, retinal damage).

Diabetic person measures blood glucose level by means of a glucometer before each meal. Based on the result of the measurement, the person applies necessary quantity of insulin. The course of daily blood glucose is much more oscillated—an ideal model is to get as near as possible to the physiological values shown in the figure of standard course of blood glucose fluctuation.

Data acquisition and experimental design

Subject of investigation: 5-to-6-year-old (at the time of recording) boy suffering from Type 1 diabetes mellitus since his age of 4.

Monitoring of blood glucose level: by patient’s mother systematically in two periods using the FreeStyle Mini glucometer by Abbott Laboratories.

Data arrangement: the recordings are grouped into two contiguous periods (denotes as GL1 and GL2 hereinafter) separated by a time gap. Experimental units are the days in which glucose level were recorded. Each period is divided into four samples: the glucose levels are read four times a day prior to breakfast (fasting level), prior to lunch, prior to dinner, and prior to sleep (PB, PL, PD, and PS hereinafter). As a result, we have eight samples to be compared.

Worth mentioning: data recorded are rather time series sampling continuous dependence of glucose level on time for a single human instead of monitoring more persons. This inevitably leads to the following two notes:

The results in any case cannot be generalized outside the subject of investigation. The purpose of the presented work is educational above all. However, as a by-product, it could serve as a guideline for the treatment discipline and the proper care.

The date/time dependence of glucose level is not treated at all. Instead of this, the statistical comparison of glucose levels at different parts of a day during periods of order several months are done. The time series analysis is beyond the scope of the present work and will be treated in another notebook. However, the comparison of statistics between the periods can be viewed as a simple example of time evolution analysis.

A few useful pieces of code and settings

The following cell gathers global settings for keeping the Mathematica codes brief but still intelligible:

(**** Initialization cell #1 ****)
Clear[binWd,boxWhiskerGlobal,cen,data,hor,homSc,k,lis,lisb,mean,rdTime,rdTimV,pat,per,pat,pos,pt,sampNo,sf,showR,showL,sigLev,tmpa,tmpb,tyData,ver,violOf]
pdfS="PDF";cdfS="CDF";pb="PB";pl="PL";pd="PD";ps="PS";labelGL="blood glucose [mmol/L]";perStr="period";rdtStr="read time";bwdStr="bin width";pvStr="p-value";tsStr="test";tabHd={{"Period 1","Period 2"},{"start-end","days"}};
ltRtMenu={0->"DistChart",1->"BWChart",2->"HistoPDF",3->"HistoCDF",4->"InvHisPDF",5->"InvHisCDF",6->"Q-Q plot",7->"P-P plot"};
imSizA=370;imSizB=750;imSizC=360;imSizD=725;frMarg=10;
periodList[data_]:=Row[{Row[First[First[data]],"/"],Row[First[Last[data]],"/"]},"-"](* periods listing formatting *)
showRecords[tyData_]:=tyData->Manipulate[Short[Style[per[[rdTime]],stdColor[[rdTime]]],6],Row[{Control[{{per,ToExpression[tyData<>"1"],perStr},{ToExpression[tyData<>"1"]->"1",ToExpression[tyData<>"2"]->"2"}}],Control[{{rdTime,1,rdtStr},{1->pb,2->pl,3->pd,4->ps}}]},Spacer[20]],Deployed->False]
coeffOfVariation[lis_]:=If[(tmpa=Mean[lis])==0,0,StandardDeviation[lis]/tmpa]
minQuartMax[lis_]:=Flatten@{Min[lis],Quartiles[lis],Max[lis]}(* min,lqrt,med,uqrt,max *)
safeSkewness[lis_]:=Sqrt[Length[lis]/6]Skewness[lis](* relative to its stderr (normality assumed) *)
kurtosis[lis_]:=Kurtosis[lis]-3(* kurtosis of the normal distribution equal to zero (exceed kurtosis) *)
safeKurtosis[lis_]:=Sqrt[Length[lis]/24]kurtosis[lis](* relative to its stderr (normality assumed) *)
numStats={{Length,Mean,TrimmedMean,Commonest,minQuartMax,InterquartileRange},{Min,Max,Variance,StandardDeviation,MeanDeviation,MedianDeviation},{coeffOfVariation,Skewness,kurtosis,safeSkewness,safeKurtosis}};
stdRed=RGBColor[.816,.071,.071];stdGreen=RGBColor[.071,.816,.071];stdBlue=RGBColor[.071,.071,.816];stdBlack=Black;stdColor={stdGreen,stdRed,stdBlue,stdBlack};
brightRed=Red;brightGreen=Green;brightBlue=Blue;brightBlack=Black;brightColor={brightGreen,brightRed,brightBlue,brightBlack};
paleRed=RGBColor[1,.5,.5];paleGreen=RGBColor[.5,1,.5];paleBlue=RGBColor[.5,.5,1];paleBlack=Gray;paleColor={paleGreen,paleRed,paleBlue,paleBlack};
fadedRed=RGBColor[1,.9,.9];fadedGreen=RGBColor[.9,1,.9];fadedBlue=RGBColor[.9,.9,1];fadedBlack=RGBColor[.9,.9,.9];fadedColor={fadedGreen,fadedRed,fadedBlue,fadedBlack};
timeOfDay[pt_]:=Which[pt==1,Style[pb,15,stdColor[[pt]]],pt==2,Style[pl,15,stdColor[[pt]]],pt==3,Style[pd,15,stdColor[[pt]]],pt==4,Style[ps,15,stdColor[[pt]]]]
showStats[sf_,data_,pt_]:=TableForm[{Through[sf[data]]},TableDirections->Column,TableHeadings->{{timeOfDay[pt]},ToString/@sf},TableSpacing->{2,2}]
cdfRange=PlotRange->{{0,28},{0,1.02}};cdfInvRange=PlotRange->{{0,1.02},{0,28}};cdfLabel=FrameLabel->{Style[labelGL,12],Style[cdfS,12]};cdfInvLabel=FrameLabel->{None,Style[labelGL,12]};
pdfRange=PlotRange->{{0,28},{0,.13}};pdfInvRange=PlotRange->{{0,.13},{0,28}};pdfLabel=FrameLabel->{Style[labelGL,12],Style[pdfS,12]};pdfInvLabel=FrameLabel->{None,Style[labelGL,12]};
binTypes={Automatic->"Auto","Sturges","Scott","FreedmanDiaconis"->"F-D","Knuth","Wand",{.5}->"  0.5",{1}->"  1",{1.5}->"  1.5",{2}->"  2",{3}->"  3",{4}->"  4",{5}->"  5",{6}->"  6",{10}->"10",{15}->"15"};
homscTests={"BrownForsythe","Conover","Levene"};homscShTsts=homscTests;
locEqTests={"CompleteBlockF","FriedmanRank","KruskalWallis","KSampleT"};locEqShTests={"KruskalWallis","KSampleT"};
histoAltTicks=FrameTicks->{{Automatic,None},{Automatic,None}};histoAltAR=AspectRatio->.8;
vRange=First[pdfRange]->First[Last[pdfRange]];vvRange=First[pdfRange]->{First[Last[pdfRange]],First[Last[pdfRange]]};
commonChartLabels=ChartLabels->{{Style[pb,16,stdColor[[1]]],Style[pl,16,stdColor[[2]]],Style[pd,16,stdColor[[3]]],Style[ps,16,stdColor[[4]]]},{Style["1",16],Style["2",16]}};
posOfPat[lis_,pat_]:=First@First@Position[lis,pat]
dateListLabel[lis_,pat_]:=Style[timeOfDay[(pos=posOfPat[lis,pat])],15,stdColor[[pos]]];
histoLabel[pt_]:=Style[timeOfDay[pt],15,stdRed];
bxWhStyle={{"Fences",Thick},{"MeanDiamond",LightGray},{"MeanMarker","+",Black},{"MedianNotch",.25,Gray},{"MedianMarker",.1,White},{"Outliers",stdBlue},{"FarOutliers",stdRed},{"Whiskers",Thin}};
pBLDScw=Table[{fadedColor[[k]],Disk[{Cos[(1/10) (5+2 k) π],Sin[(1/10) (5+2 k) π]},.25],Text[timeOfDay[k],{Cos[(1/10) (5+2 k) π],Sin[(1/10) (5+2 k) π]}]},{k,4}];
pBLDSccw=Table[{fadedColor[[k]],Disk[{Cos[(1/10) (5-2 k) π],Sin[(1/10) (5-2 k) π]},.25],Text[timeOfDay[k],{Cos[(1/10) (5-2 k) π],Sin[(1/10) (5-2 k) π]}]},{k,4}];
pBLDStxt={Text[Style["Period 1",16],{-.6,1.2}],Text[Style["Period 2",16],{.6,1.2}]};pBLDShead=Arrowheads[{-.06,.06}];
testArrowV[hor_,ver_,cen_]:=Arrow[BezierCurve[{{hor,ver},{cen,0},{hor,-ver}}]]
testArrowH[hor_,ver_,cen_]:=Arrow[BezierCurve[{{hor,ver},{0,cen},{-hor,ver}}]]
sigLevList={.01->"1",.05->"5",.1->"10"};sigLevInit=.05;sigLevStr="sig.lev.%";
SetOptions[Plot,Axes->False,Frame->False];
SetOptions[ListPlot,PlotRange->{{0,28},{0,28}},PlotStyle->{stdBlack,PointSize[.01]},Axes->False,Frame->True,FrameLabel->{Style[labelGL,12],Style[labelGL,12]},ImageSize->imSizC,AspectRatio->.8];
SetOptions[DateListPlot,vRange,ImageSize->imSizA,Joined->False,Filling->Bottom,AspectRatio->.5,DateTicksFormat->{"Year","/","Month"},FrameLabel->{Style["calendar date",12],Style[labelGL,12]}];
SetOptions[DateListLogPlot,vRange,ImageSize->imSizA,Joined->False,Filling->Bottom,AspectRatio->.5,DateTicksFormat->{"Year","/","Month"},FrameLabel->{Style["calendar date",12],Style[labelGL,12]}];
SetOptions[BoxWhiskerChart,ChartLabels->{Style[pb,12,stdGreen],Style[pl,12,stdRed],Style[pd,12,stdBlue],Style[ps,12,Black]},ChartStyle->{RGBColor[1,.8,0]},BarSpacing->1,FrameLabel->{None,Style[labelGL,12]},GridLines->Automatic,GridLinesStyle->Directive[fadedBlack,Dashed],ImageSize->imSizA,AspectRatio->.8];
SetOptions[DistributionChart,ChartLabels->{Style[pb,12,stdGreen],Style[pl,12,stdRed],Style[pd,12,stdBlue],Style[ps,12,Black]},ChartStyle->RGBColor[1,.8,0],BarSpacing->1,ChartBaseStyle->EdgeForm[],ChartElementFunction->"PointDensity",FrameLabel->{None,Style[labelGL,12]},ImageSize->imSizA,AspectRatio->.8];
SetOptions[Histogram,ChartElementFunction->"FadingRectangle",Axes->False,Frame->True,ImageSize->imSizC,AspectRatio->1/GoldenRatio];
SetOptions[QuantilePlot,vvRange,PlotStyle->stdRed,Joined->False,Filling->None,FrameLabel->{None,Style["empirical quantiles",12]},ImageSize->imSizC,AspectRatio->1/GoldenRatio];
SetOptions[ProbabilityPlot,PlotStyle->stdRed,Joined->False,Filling->None,FrameLabel->{None,Style["empirical probabilities",12]},ImageSize->imSizC,AspectRatio->1/GoldenRatio];
SetOptions[BarChart,Axes->False,Frame->True,PlotRange->All,FrameLabel->{Style[pvStr,12],Style[tsStr,12]},FrameTicks->{{None,None},{Automatic,None}},BarOrigin->Left,BarSpacing->.5,ChartElementFunction->"Rectangle",ChartBaseStyle->EdgeForm[None],ImageSize->imSizC,AspectRatio->1/GoldenRatio];
SetOptions[Manipulate,ControlPlacement->Top,AppearanceElements->"ResetButton",FrameMargins->Automatic,ContinuousAction->True,Deployed->True,SaveDefinitions->True];
Quiet@Needs["MultivariateStatistics`"]
pValCorrelation[lis_,lisb_]:=Module[{len,cor,tst},len=Length[lis];cor=Correlation[lis,lisb];If[lis==lisb,0,(tst=cor (Sqrt[((len-2)/(1-(cor^2)))]);1-CDF[StudentTDistribution[len-2],tst])]]
pValSpearmanRankCorrelation[lis_,lisb_]:=Module[{len,cor,tst},len=Length[lis];cor=N@SpearmanRankCorrelation[lis,lisb];If[lis==lisb,0,(tst=cor (Sqrt[((len-2)/(1-(cor^2)))]);1-CDF[StudentTDistribution[len-2],tst])]]
pValKendallRankCorrelation[lis_,lisb_]:=Module[{len,cor,sdv},len=Length[lis];cor=N@KendallRankCorrelation[lis,lisb];If[lis==lisb,0,(sdv=(Sqrt[((4 len+10)/(9 len (len-1)))]);1-CDF[NormalDistribution[0,sdv],cor])]]

Data import and preprocessing

The data was first exported from the .xlsx to .csv format using OpenOffice Calc.

Data import into Mathematica (Mangano, 2010):

(**** Initialization cell #2 ****)
GlucoLev1=Import[FileNameJoin[{NotebookDirectory[],"GlucoLev_1.csv"}],"DateStringFormat"->{"Year"," ","Month"," ","Day"}];
GlucoLev2=Import[FileNameJoin[{NotebookDirectory[],"GlucoLev_2.csv"}],"DateStringFormat"->{"Year"," ","Month"," ","Day"}];
GlucoLev=Join[GlucoLev1,GlucoLev2];
Panel[TraditionalForm@Style[TableForm[Transpose[{periodList/@{GlucoLev1,GlucoLev2},(glLen=Length/@{GlucoLev1,GlucoLev2})}],TableHeadings->tabHd],"Text",FontSize->18],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD]

Splitting data into samples suitable for further treatment. Indices 1, 2, 3, 4 correspond to PB (fasting level), PL, PD, and PS, respectively :

(**** Initialization cell #3 ****)
GL1=Table[GlucoLev1[[All,i]],{i,2,5}];(* gl.lev.: 1-col. fmt for each part of a day, per. 1 *)
GL2=Table[GlucoLev2[[All,i]],{i,2,5}];(* same for period 2 *)
GL=Table[GlucoLev[[All,i]],{i,2,5}];(* amalgamation of both periods *)
DateVsGL1=Table[GlucoLev1[[All,{1,i}]],{i,2,5}];(* date-vs-level: 2-column format, Period 1 *)
DateVsGL2=Table[GlucoLev2[[All,{1,i}]],{i,2,5}];(* same for Period 2 *)
DateVsGL=Table[GlucoLev[[All,{1,i}]],{i,2,5}];(* amalgamation of both periods *)

What hypotheses are subject to testing?

The crucial point of interest are two four-sample null hypotheses that the daytime means/medians within each period are equal (McKillup, 2006),

(Subscript[H, 0]):(Subscript[μ, PB]) = (Subscript[μ, PL]) = (Subscript[μ, PD]) = (Subscript[μ, PS])            (Subscript[H, 0]):(Subscript[ν, PB]) = (Subscript[ν, PL]) = (Subscript[ν, PD]) = (Subscript[ν, PS])

where (Subscript[μ, PB]), (Subscript[μ, PL]), (Subscript[μ, PD]), (Subscript[μ, PS]) ((Subscript[ν, PB]), (Subscript[ν, PL]),(Subscript[ν, PD]), (Subscript[ν, PS])) denote corresponding means/medians in Period 1 (Period 2), respectively, against the alternative hypotheses, (Subscript[H, A])=(Subscript[(Overscript[H, –]), 0]). Test like this will be carried out for each period.

Panel[Graphics[{pBLDScw,pBLDSccw,pBLDStxt,pBLDShead,paleBlack,Thick,testArrowV[-.4,.66,.2],testArrowV[-.71,.26,.5],testArrowV[.4,.66,-.2],testArrowV[.71,.26,-.5],Text[Style["two four-sample tests",16,stdBlack],{0,0},{0,0},{0,1}]}],ImageSize->imSizD,Alignment->{Center,Center},Background->White]

Second, we shall deal with four two-sample null hypotheses that the means/medians across the Periods 1 and 2 are equal,

(Subscript[H, 0]):(Subscript[μ, PB]) = (Subscript[ν, PB])            (Subscript[H, 0]):(Subscript[μ, PL]) = (Subscript[ν, PL])            (Subscript[H, 0]):(Subscript[μ, PD]) = (Subscript[ν, PD])            (Subscript[H, 0]):(Subscript[μ, PS]) = (Subscript[ν, PS])

where ((Subscript[μ, PB]), (Subscript[ν, PB])), ((Subscript[μ, PL]), (Subscript[ν, PL])), ((Subscript[μ, PD]), (Subscript[ν, PD])), and ((Subscript[μ, PS]), (Subscript[ν, PS])) denote corresponding pairs of means/medians in Periods 1 and 2, against the alternative hypotheses, (Subscript[H, A])=(Subscript[(Overscript[H, –]), 0]). Test like this will be carried out for each measuring time PB, PL, PD, and PS (“cross-period” comparison).

Panel[Graphics[{pBLDScw,pBLDSccw,pBLDStxt,pBLDShead,paleBlack,Thick,testArrowH[-.4,.66,.1],testArrowH[-.4,-.66,-.1],testArrowH[-.71,.3,.05],testArrowH[-.71,-.3,-.05],Text[Style["four two-sample tests",16,stdBlack],{0,0},{0,0},{1,0}]}],ImageSize->imSizD,Alignment->{Center,Center},Background->White]

We intentionally write “means/medians” since we don’t know yet which method (ANOVA, nonparametric Kruskal–Wallis test for independent samples or nonparametric Friedman test for dependent samples in case of four-sample hypotheses, and Student t-test or nonparametric Mann–Whitney test in case of two-sample hypotheses) will eventually be used. Some test are more sensitive, and some are less sensitive to  significant disruption of independence, normality and/or homoscedasticity of the sample distributions. Hence, tests for these properties must be carried out prior to tests for equal means/medians.

Descriptive statistics

In this section, basic descriptive statistics with both numerical and graphical output shall be presented. This is crucial for visual assessment of the data, and for later selection of proper tools of inferential statistics.

Displaying numerical values

A look at the raw data is not particularly useful except that one sees the structure of the arrays GLperiod[[rdTime]] and DateVsGLperiod[[rdTime]] defined above is as follows. Remember that period ∈ {1,2} and rdTime ∈ {PB,PL,PD,PS}:

TabView[showRecords/@{"GL","DateVsGL"}]

Time dependence of data

In the above subsection Data acquisition and experimental design we intentionally relinquished carrying out any time dependence analysis. However, the plots of the temporal glucose level evolution (including the gap between the periods) are worth displaying for an overall idea.

Manipulate[If[rdTime==5,Grid@Partition[(* all *)
Table[DateListPlot[data,FillingStyle->paleColor[[posOfPat[DateVsGL,data]]],
Epilog:>Text[dateListLabel[DateVsGL,data],{{2008,7,1},25}]],
{data,DateVsGL}],2],
DateListPlot[DateVsGL[[rdTime]],(* per one *)
FillingStyle->paleColor[[rdTime]],ImageSize->imSizB]],
{{rdTime,5,rdtStr},{1->pb,2->pl,3->pd,4->ps,5->"all"}}]

Charting and data visualization

Mathematica makes possible a flexible and convenient approach to charting and information visualization, providing high levels of automation. She lets you immediately take data and produce compelling dynamic visualizations. Statistical visualization is used to understand how data is distributed and how they compare to other datasets or distributions.

Box-and-whisker and distribution charts compare a number of data distributions to each other.

Histograms effectively estimate the various distribution functions through binning the data.

Quantile-quantile and probability-probability plots compare data to a reference distribution.

The following additions to the Mathematica built-ins were adopted:

The coefficient of variation is defined as the ratio of the standard deviation to the mean.

The built-in Kurtosis provided by Mathematica is actually the fourth standardized moment (Subscript[β, 2]). The redefined kurtosis (sometimes called exceed kurtosis) (Subscript[γ, 2])≡(Subscript[β, 2])-3 is a correction to make the kurtosis value of the normal distribution equal to zero (Nelson, 2004).

The standard error of skewness and (exceed) kurtosis, as an estimator of these statistics of an underlying normal distribution, is (((6/N))^(1/2)) and (((24/N))^(1/2)), respectively. For that reason we introduce safeSkewness and safeKurtosis defined as Skewness and kurtosis divided by the above standard error estimators. It is good practice to rely on skewness and kurtosis only when the safe variants are at least of several units of, better, much greater then 1 (Press, 2007).

Manipulate[𝒩=DistributionFitTest[per[[rdTime]],
Automatic,"HypothesisTestData",Method->Automatic];
GraphicsGrid[{{Which[
showL==0,DistributionChart[Thread[Style[per,paleColor]]],
showL==1,BoxWhiskerChart[Thread[Style[per,paleColor]],bxWhStyle],
showL==2,Show[Histogram[per[[rdTime]],binWd,pdfS,pdfRange,pdfLabel,histoAltAR,ChartStyle->paleColor[[rdTime]]],Plot[PDF[𝒩["FittedDistribution"],x],{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showL==3,Show[Histogram[per[[rdTime]],binWd,cdfS,cdfRange,cdfLabel,histoAltAR,ChartStyle->paleColor[[rdTime]]],Plot[CDF[𝒩["FittedDistribution"],x],
{x,First@cdfRange[[2,1]],Last@cdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showL==4,Show[Histogram[per[[rdTime]],binWd,pdfS,BarOrigin->Left,pdfInvRange,ChartStyle->paleColor[[rdTime]],histoAltTicks,pdfInvLabel,histoAltAR,Epilog->{Text[histoLabel[rdTime],{.11,25}],Text[Style[pdfS,12],{.065,1.5}]}],ParametricPlot[{PDF[𝒩["FittedDistribution"],x],x},{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showL==5,Show[Histogram[per[[rdTime]],binWd,cdfS,BarOrigin->Left,cdfInvRange,ChartStyle->paleColor[[rdTime]],histoAltTicks,cdfInvLabel,histoAltAR,Epilog->{Text[histoLabel[rdTime],{.863,3}],Text[Style[cdfS,12],{.51,1.5}]}],ParametricPlot[{CDF[𝒩["FittedDistribution"],x],x},
{x,First@cdfRange[[2,1]],Last@cdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showL==6,QuantilePlot[per[[rdTime]],PlotStyle->stdColor[[rdTime]],histoAltAR,Epilog->{Text[histoLabel[rdTime],{23.7,3}],Text[Style["theoretical quartiles",12],{14,1.5}]}],
showL==7,ProbabilityPlot[per[[rdTime]],PlotStyle->stdColor[[rdTime]],histoAltAR,Epilog->{Text[histoLabel[rdTime],{.846,.107}],Text[Style["theoretical probabilities",12],{.5,.054}]}]],(* </Which> *)
Dynamic[Which[
showR==0,DistributionChart[Thread[Style[per,paleColor]]],
showR==1,BoxWhiskerChart[Thread[Style[per,paleColor]],bxWhStyle],
showR==2,Show[Histogram[per[[rdTime]],binWd,pdfS,pdfRange,pdfLabel,histoAltAR,ChartStyle->paleColor[[rdTime]]],Plot[PDF[𝒩["FittedDistribution"],x],{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showR==3,Show[Histogram[per[[rdTime]],binWd,cdfS,cdfRange,cdfLabel,histoAltAR,ChartStyle->paleColor[[rdTime]]],Plot[CDF[𝒩["FittedDistribution"],x],
{x,First@cdfRange[[2,1]],Last@cdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showR==4,Show[Histogram[per[[rdTime]],binWd,pdfS,BarOrigin->Left,pdfInvRange,ChartStyle->paleColor[[rdTime]],histoAltTicks,pdfInvLabel,histoAltAR,Epilog->{Text[histoLabel[rdTime],{.11,25}],Text[Style[pdfS,12],{.065,1.5}]}],ParametricPlot[{PDF[𝒩["FittedDistribution"],x],x},{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showR==5,Show[Histogram[per[[rdTime]],binWd,cdfS,BarOrigin->Left,cdfInvRange,ChartStyle->paleColor[[rdTime]],histoAltTicks,cdfInvLabel,histoAltAR,Epilog->{Text[histoLabel[rdTime],{.863,3}],Text[Style[cdfS,12],{.51,1.5}]}],ParametricPlot[{CDF[𝒩["FittedDistribution"],x],x},
{x,First@cdfRange[[2,1]],Last@cdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],
showR==6,QuantilePlot[per[[rdTime]],PlotStyle->stdColor[[rdTime]],histoAltAR,Epilog->{Text[histoLabel[rdTime],{23.7,3}],Text[Style["theoretical quartiles",12],{14,1.5}]}],
showR==7,ProbabilityPlot[per[[rdTime]],PlotStyle->stdColor[[rdTime]],histoAltAR,Epilog->{Text[histoLabel[rdTime],{.846,.107}],Text[Style["theoretical probabilities",12],{.5,.054}]}]
]](* </Which,Dynamic> *)},{Dynamic[Panel[Column[Table[Module[{sf},sf=numStats[[k]];
showStats[sf,per[[rdTime]],rdTime]],{k,3}],Frame->All,FrameStyle->{stdColor[[rdTime]],Thick},Background->fadedColor[[rdTime]],Spacings->2],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Center,Center},ImageSize->700]],SpanFromLeft}},
Spacings->{20,0},ImageSize->imSizB],(* </GraphicsRow> *)
Row[{Control[{{per,GL1,perStr},{GL1->"1",GL2->"2"}}],
Control[{{rdTime,1,rdtStr},{1->pb,2->pl,3->pd,4->ps}}],
Control[{{showL,1,"left"},ltRtMenu}],
Control[{{showR,4,"right"},ltRtMenu}],
Control[{{binWd,{2},bwdStr},binTypes}]},Spacer[6]]]

boxWhiskerGlobal=BoxWhiskerChart[Table[{Style[GL1[[i]],paleColor[[i]]],Style[GL2[[i]],paleColor[[i]]]},{i,4}],bxWhStyle,commonChartLabels,LabelStyle->16,FrameLabel->{None,Style[labelGL,16]},ImageSize->imSizB,AspectRatio->.5]

Observations

The probability distribution functions (PDF, left column) seem to eyes to be distinct from normal distribution (dashed lines) whatever its parameters may be (those drawn in histograms use sample mean and standard deviation from the numerical statistics above), on the other hand, the cumulative distribution functions (CDF, right column) seem to be in reasonable accord with normal distribution.

The PDF histograms and Q-Q plots indicate (or, at least they seem to eyes) that the sample distributions are different from normal distribution, maybe with the exceptions of

PB in both periods,

PS in Period 2.

On the other hand, the CDF histograms seem in any case to be in reasonable accord with normal distribution, which can be misleading when assessing normality at a rough guess.

The numerical statistics as well as the box-and-whisker chart show considerable variability of glucose levels in both periods and measuring times:

Minimum and maximum values differ almost (in some cases more than) by one order.

Coefficient of variation ranges approximately from 1/3 to 1/2, usually exceeding the latter value.

Interquartile range’s (IQR) extent is approximately 50–100% of the mean.

With one exception for each period, no outliers appear.

The box-and-whisker chart  indicates that the means (or medians) might differ significantly with the following exceptions:

PB between Periods 1 and 2,

PD between Periods 1 and 2,

between PL and PS for Period 1,

between PB, PL and PS for Period 2.

These results can be summarized into “suspicion” that the glucose level PD mean/median significantly differs from the remaining (PB, PL and PS) ones; this applies for Period1 and particularly for Period 2.

The safeSkewness defined above is mostly too small to be trustworthy; the exceptions are

PL, PD and PS in Period 1,

PD in Period 2,

all having markedly positive safe skewness within approximate interval (2.6,4.5).

The safeKurtosis defined above is mostly negative within approximate interval (-1.3,-2.4), i.e., the samples are slightly platykurtic, but such values are too small to be trustworthy. The exceptions from this rule are:

PD and PS in Period 1,

PD in Period 2,

all having negligible kurtosis. We can conclude that with the exception of the PL sample in Period 1 (which is platykurtic and positively skewed), all samples have either significant positive skewness (see 4. above) and negligible kurtosis, or they are symmetric but slightly platykurtic.

The result of the built-in Commonest must be treated cautiously. It simply returns a list of the elements that are the most common in a given list. But this is not the same as the meaningful notion of the mode, which can be demonstrated by choosing smaller bin widths and observing the resulting appearance. A typical example is the value of the mode about 15 as determined “by a bit of common sense” from the PDF histogram with a reasonable bin width (2, according to, e.g., Sturges or Freedman–Diaconis rule), while the Commonest provides half the value.

These vague and subjective statements based on visual judgment will be quantitatively reviewed later in section on Inferential statistics.

Guidelines for test selection

For exploration of hypothesis (1) in each period, the test candidates to be used are multiple-sample tests like one-way ANOVA, Kruskal–Wallis test and Friedman test (Litschmannová, 2011). However, it is still unclear which one., e.g., whether ANOVA is excluded for possible violation of independence, normality or homoscedasticity. Kruskal–Wallis test and Friedman test are good candidates also from the point of view of potential violation of these properties.

For exploration of hypothesis (2) in each of four daily measuring times, the test candidates are either Student t-test or Mann–Whitney test (Litschmannová, 2011). The reason is similar as in the previous paragraph: still unclear violation normality and homoscedasticity.

Aside from these two guidelines for subsequent inferential analysis, the descriptive statistics provide an intuitive picture of the evolution of diabetes compensation between two reported periods.

Conclusions of descriptive statistics

Period 1: While the PB glucose level maximum is about 15 mmol/L and skewness is insignificantly negative, the remaining PL, PD and PS maxima (at least the main ones) are close to physiological glucose level 5 mmol/L, and skewness is markedly positive. This facts can be interpreted as follows: while for PL, PD and PS the diabetes is more or less compensated, the PB level (the fasting level) is considerably shifted towards greater, uncompensated value. The reason might be morning insulin insufficiency, after longer time without insulin dose than during daytime.

Period 2: While for slightly bimodal PB glucose level distribution is almost symmetrical and the mode fluctuates within 13 mmol/L and 17 mmol/L, for markedly bimodal PL distribution even within 3 mmol/L and 15 mmol/L. Only distinctly unimodal PD distribution with significant positive skewness succeeds in keeping the mode at compensated level of about 5 mmol/L. The PS distribution’s mode is even about 17 mmol/L.

These observations can be interpreted as an indication that the behavior of the glucose level during the Period 2 may be evaluated markedly worse than during Period 1.

Inferential statistical analysis

In this section, the glucose level samples will be subjected to suitably chosen statistical tests, following the schemes of four-sample and two-sample tests (Ruskeepää, 2009). First, we shall test violation of the assumptions of parametric tests like ANOVA and Student t-test: independence, normality and homoscedasticity of samples.

Testing assumptions: Independence of samples

Question: are the PB, PL, PD, and PS glucose level samples within each period independent or not?

Answer: it is not trivial, since the readings form a time series (four subsequent records a day chained day by day). Each day, an unknown factor may emerge, which can influence all four readings in that day, e.g., an illness, inaccessibility of insulin in the worst case). Based on this, one can rather be inclined to conclude that the samples are not independent.

On the other hand, characteristic timescale of glucose level variability is small compared to 24 hours scale, so a short-term event during a day is able to alter the level significantly only for one particular monitored part of the day (e.g., stress, sports, exhaustion, consumption of an inappropriate meal, etc.). This  considerations make us incline to think of the samples as truly independent. Since none of the daylong factors is supposed to take place, we might conclude that the PB, PL, PD, and PS samples are independent.

So it is useful to carry out a quantitative test that there is no association between the samples. We use primarily the Kendall’s rank correlation test for its robustness and because the samples mostly are not normally distributed.

Manipulate[ℛ=LinearModelFit[Transpose[{per[[rdTime]],per[[rdTimV]]}],x,x];
Grid[{{Show[
ListPlot[Transpose[{per[[rdTime]],per[[rdTimV]]}],FrameStyle->{{Directive[Thick,stdColor[[rdTimV]]],Automatic},{Directive[Thick,stdColor[[rdTime]]],Automatic}},
Epilog->{Text[histoLabel[rdTime],{26.2,1.6}],Text[histoLabel[rdTimV],{1.6,26}]}],
Plot[ℛ["BestFit"],{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->{Thickness[.01],paleBlack,Opacity[.5]}]],
Column[{Panel[TraditionalForm@
Style[TableForm[{{N[Correlation[per[[rdTime]],per[[rdTimV]]]],N[pValCorrelation[per[[rdTime]],per[[rdTimV]]]]},
{N[SpearmanRankCorrelation[per[[rdTime]],per[[rdTimV]]]],N[pValSpearmanRankCorrelation[per[[rdTime]],per[[rdTimV]]]]},{N[KendallRankCorrelation[per[[rdTime]],per[[rdTimV]]]],
(tmpb=N[pValKendallRankCorrelation[per[[rdTime]],per[[rdTimV]]]])}},
TableHeadings->{{"Pearson's r","Spearman's ρ","Kendall's τ"},{"corr.coeff.","p-value"}}],"Text",FontSize->12],
FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizC],
Show[BarChart[Reverse[{pValCorrelation[per[[rdTime]],per[[rdTimV]]],pValSpearmanRankCorrelation[per[[rdTime]],per[[rdTimV]]],pValKendallRankCorrelation[per[[rdTime]],per[[rdTimV]]]}],ChartLabels->Placed[Reverse[{"Pearson","Spearman","Kendall"}],Left],ChartStyle->fadedBlack,AspectRatio->.46,ImageSize->imSizC],Graphics[{Dashed,stdBlack,Line[{{sigLev,0},{sigLev,4}}]}]]}]
},{Panel[TraditionalForm@Style["The null hypothesis that there is no association between the data is "<>If[tmpb>=sigLev,"not ",""]<>"rejected at the "<>ToString[100sigLev]<>" percent level based on the Kendall's correlation test.",(* bottom row *)
14],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD],SpanFromLeft}
},Alignment->{{Right,Left},{Top,Top}}], (* </Grid> *)
Row[{Control[{{per,GL1,perStr},{GL1->"1",GL2->"2"}}],
Control[{{rdTime,1,"h-axis"},{1->pb,2->pl,3->pd,4->ps}}],
Control[{{rdTimV,2,"v-axis"},{1->pb,2->pl,3->pd,4->ps}}],
Control[{{sigLev,sigLevInit,sigLevStr},sigLevList}]
},Spacer[20]]]

Conclusion

We do not reject the null hypothesis that there is no association between the following pairs of glucose level samples at the 5% significance level, i.e., these pairs of samples can be considered as independent in Period 1:

PB vs. PD,

PB vs. PS.

All remaining 4 pairs of Period 1, i.e., PB vs. PL, PL vs. PD, PL vs. PS, and PD vs. PS, are significantly dependent.

Similarly, the following pairs of samples can be considered as independent for Period 2:

PB vs. PL,

PB vs. PD,

PB vs. PS,

PL vs. PD.

The remaining 2 pairs of Period 2, i.e., PL vs. PS and PD vs. PS, are significantly dependent.

Testing assumptions: Normality

Normality of data is one of the basic assumptions if one wants to use parametric tests. Although normality is not a critical assumption with sample sizes of order (10^2) as in this case, it is worth testing it anyway.

Manipulate[𝒩=DistributionFitTest[per[[rdTime]],
Automatic,"HypothesisTestData",Method->Automatic,SignificanceLevel->sigLev];
Grid[{If[showR≠0,
{QuantilePlot[per[[rdTime]],PlotStyle->stdColor[[rdTime]],(* upper row: QQ/PP *)Epilog->{Text[histoLabel[rdTime],{23.7,3}],Text[Style["theoretical quartiles",12],{14,1.5}]}],
ProbabilityPlot[per[[rdTime]],PlotStyle->stdColor[[rdTime]],Epilog->{Text[histoLabel[rdTime],{.846,.107}],Text[Style["theoretical probabilities",12],{.5,.054}]}]},
{Show[Histogram[per[[rdTime]],binWd,pdfS,(* upper row: histograms *)
pdfRange,pdfLabel,ChartStyle->paleColor[[rdTime]]],Plot[PDF[𝒩["FittedDistribution"],x],{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]],Show[Histogram[per[[rdTime]],binWd,cdfS,cdfRange,cdfLabel,ChartStyle->paleColor[[rdTime]]],
Plot[CDF[𝒩["FittedDistribution"],x],
{x,First@cdfRange[[2,1]],Last@cdfRange[[2,1]]},PlotStyle->{Thickness[.01],brightColor[[rdTime]]}]]}],
{Framed[Quiet@𝒩["TestDataTable",All],FrameMargins->frMarg],(* middle row *)
Show[BarChart[Reverse[𝒩["PValue",All]],ChartLabels->Placed[Reverse[𝒩["AllTests"]],Left],ChartStyle->paleColor[[rdTime]],AspectRatio->.4,ImageSize->360],Graphics[{Dashed,stdBlack,Line[{{sigLev,0},{sigLev,6}}]}]]},
{Panel[TraditionalForm@Style[𝒩["TestConclusion",Automatic],(* bottom row *)
14],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD],SpanFromLeft}
},Alignment->{{Right,Left},{Bottom,Top}}],(* </Grid> *)
Row[{Control[{{per,GL1,perStr},{GL1->"1",GL2->"2"}}],
Control[{{rdTime,1,rdtStr},{1->pb,2->pl,3->pd,4->ps}}],
Control[{{binWd,{2},bwdStr},binTypes}],
Control[{{sigLev,sigLevInit,sigLevStr},sigLevList}],
Control[{{showR,0,"plot"},{0->"Hist",1->"Qq/Pp"}}]},Spacer[6]]]

Conclusion

We reject the null hypothesis that glucose level samples are drawn from a normal distribution at the 5% significance level for all samples with the exceptions of:

PB in Period 1,

PB and PS in Period 2,

which can be considered as normally distributed.

Testing assumptions: Homoscedasticity

The built-in suite of variance equivalence test performs a hypothesis test on the data with null hypothesis that the true population variances are equal, and alternative hypothesis that at least one is different.

Since normality is violated for at least two read times in each period and for at least one sample at each read time (except of PB),  we exclude the Bartlett’s and Fisher’s tests which are sensitive to a normality assumption.

Manipulate[𝒩=Table[DistributionFitTest[per[[rdTime]],Automatic,"HypothesisTestData",Method->Automatic],{per,{GL1,GL2}},{rdTime,4}];
(Subscript[ℋ, 4])=VarianceEquivalenceTest[ToExpression["GL"<>ToString[per]],"HypothesisTestData",Method->Automatic,SignificanceLevel->sigLev];
(Subscript[ℋ, 2])=VarianceEquivalenceTest[{GL1[[rdTime]],GL2[[rdTime]]},"HypothesisTestData",Method->Automatic,SignificanceLevel->sigLev];
Grid[{If[homSc==4,
{Plot[{PDF[𝒩[[per,1]]["FittedDistribution"],x],
PDF[𝒩[[per,2]]["FittedDistribution"],x],
PDF[𝒩[[per,3]]["FittedDistribution"],x],
PDF[𝒩[[per,4]]["FittedDistribution"],x]},
{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->Table[{Thickness[.01],stdColor[[k]]},{k,4}],Frame->True,PlotRange->Last@pdfRange,FrameLabel->Last@pdfLabel,ImageSize->imSizC],
Plot[{CDF[𝒩[[per,1]]["FittedDistribution"],x],
CDF[𝒩[[per,2]]["FittedDistribution"],x],
CDF[𝒩[[per,3]]["FittedDistribution"],x],
CDF[𝒩[[per,4]]["FittedDistribution"],x]},{x,First@cdfRange[[2,1]],Last@cdfRange[[2,1]]},PlotStyle->Table[{Thickness[.01],stdColor[[k]]},{k,4}],Frame->True,PlotRange->Last@cdfRange,FrameLabel->Last@cdfLabel,ImageSize->imSizC]},
{Plot[{PDF[𝒩[[1,rdTime]]["FittedDistribution"],x],
PDF[𝒩[[2,rdTime]]["FittedDistribution"],x]},
{x,First@pdfRange[[2,1]],Last@pdfRange[[2,1]]},PlotStyle->Table[{Thickness[.01],stdColor[[rdTime]]},{2}],Frame->True,PlotRange->Last@pdfRange,FrameLabel->Last@pdfLabel,ImageSize->imSizC],
Plot[{CDF[𝒩[[1,rdTime]]["FittedDistribution"],x],CDF[𝒩[[2,rdTime]]["FittedDistribution"],x]},{x,First@cdfRange[[2,1]],Last@cdfRange[[2,1]]},PlotStyle->Table[{Thickness[.01],stdColor[[rdTime]]},{2}],Frame->True,PlotRange->Last@cdfRange,FrameLabel->Last@cdfLabel,ImageSize->imSizC]
}] ,(* </If> *)
If[homSc==4,
{Framed[Quiet@(Subscript[ℋ, 4])["TestDataTable",
homscTests],FrameMargins->frMarg],Quiet@Show[BarChart[
Reverse[(Subscript[ℋ, 4])["PValue",homscShTsts]],ChartLabels->Placed[Reverse[homscShTsts],Left],
ChartStyle->fadedColor[[4]],AspectRatio->.26,ImageSize->360],Graphics[{Dashed,stdBlack,Line[{{sigLev,0},{sigLev,Length[homscTests]+1}}]}]]},
{Framed[Quiet@(Subscript[ℋ, 2])["TestDataTable",
homscTests],FrameMargins->frMarg],Quiet@Show[BarChart[Reverse[Quiet@(Subscript[ℋ, 2])["PValue",
homscTests]],ChartLabels->Placed[Reverse[homscTests],Left],ChartStyle->paleColor[[rdTime]],AspectRatio->.26,ImageSize->360],Graphics[{Dashed,stdBlack,Line[{{sigLev,0},{sigLev,Length[homscTests]+1}}]}]]}],
If[homSc==4,
{Panel[TraditionalForm@Style[(Subscript[ℋ, 4])["TestConclusion",Automatic],14],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD],SpanFromLeft},
{Panel[TraditionalForm@Style[(Subscript[ℋ, 2])["TestConclusion",Automatic],
14],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD],SpanFromLeft}]},
Alignment->{{Right,Left},{Bottom,Top}}],
Row[{Control[{{homSc,4,"across"},{4->"read times",2->"periods"}}],
Control[{{per,1,perStr},{1->"1",2->"2"}}],
Control[{{rdTime,1,rdtStr},{1->pb,2->pl,3->pd,4->ps}}],
Control[{{sigLev,sigLevInit,sigLevStr},sigLevList}]},Spacer[6]]]

Conclusion

For each period we reject the null hypothesis that the variances of PB, PL, PD and PS are equal at the 5% significance level. (Homoscedasticity of each period is disproved.)

For each read time PB and PL we reject the null hypothesis that their variances in both periods are equal at the 5% significance level.

For each read time PD and PS we do not reject the null hypothesis that their variances in both periods are equal at the 5% significance level. The PD and PS samples can be considered as homoscedastic across periods.

Testing assumptions: Summary

Manipulate[Graphics[{pBLDScw,pBLDSccw,Text[Style["Period 1",16],{-.6,1.2}],Text[Style["Period 2",16],{.6,1.2}],
If[violOf==1,
{pBLDShead,paleBlack,Thick,testArrowV[-.75,.18,-.6],
Arrow[BezierCurve[{{-.55,.58},{-.45,.33},{-.7,.3}}]],
Arrow[BezierCurve[{{-.55,-.58},{-.45,-.33},{-.7,-.3}}]],
If[sigLev≠.01,Arrow[BezierCurve[{{-.72,.25},{-.2,0},{-.46,-.6}}]]],
Arrow[BezierCurve[{{.55,-.58},{.45,-.33},{.7,-.3}}]],
Arrow[BezierCurve[{{.72,.25},{.2,0},{.46,-.6}}]]}],
If[violOf==2,{If[sigLev≠.1,Table[{paleBlack,Thickness[.02],Circle[{Cos[(((5+2k) π)/10)],Sin[(((5+2k) π)/10)]},.25]},{k,2,4}],Table[{paleBlack,Thickness[.02],Circle[{Cos[(((5+2k) π)/10)],Sin[(((5+2k) π)/10)]},.25]},{k,1,4}]],
If[sigLev≠.1,Table[{paleBlack,Thickness[.02],Circle[{Cos[(((5-2k) π)/10)],Sin[(((5-2k) π)/10)]},.25]},{k,2,3}],Table[{paleBlack,Thickness[.02],Circle[{Cos[(((5-2k) π)/10)],Sin[(((5-2k) π)/10)]},.25]},{k,2,4}]]}],
If[violOf==3,
{pBLDShead,paleBlack,Thick,testArrowV[-.4,.66,.2],testArrowV[-.71,.26,.5],testArrowV[.4,.66,-.2],testArrowV[.71,.26,-.5],
If[sigLev≠.01,{testArrowH[-.35,.8,.7],testArrowH[-.71,.36,.45]}]}]}],
Row[{Control[{{violOf,1,"violation of"},
{1->"independence",2->"normality",3->"homoscedasticity"}}],
Control[{{sigLev,sigLevInit,sigLevStr},sigLevList}]},Spacer[10]],Alignment->Center]

Tests for equal means/medians

This is the key point of our effort.

The four-sample hypotheses

The above results indicate that the use of ANOVA for testing null hypotheses (1) is inappropriate due to

normality violation (not the main reason since ANOVA is not too much sensitive to a normality assumption, particularly when the samples size is greater than 30),

significant disruption of  independence and homoscedasticity.

Hence, we test the null hypothesis (1) using two nonparametric methods:

Friedman test for testing zero median difference of two or more non-independent samples (see also discussion in the subsection Independence of samples).

Kruskal–Wallis test for testing zero median difference (aka nonparametric ANOVA, in fact an extension of nonparametric Mann–Whitney test to more than two samples).

The post-hoc analysis is carried out using the nonparametric Kruskal–Wallis test or parametric K-Sample T test (whichever is autodetected by Mathematica as more appropriate) with Bonferroni correction.

The two-sample hypotheses

Null hypotheses (2) that the daytime means/medians across each period are equal will also be done using

Kruskal–Wallis median test (Mann–Whitney test in case of median comparison of two samples),

K-Sample T mean test,

(whichever is autodetected by Mathematica as more appropriate).

Manipulate[
(Subscript[ℒ, 4])=LocationEquivalenceTest[ToExpression["GL"<>ToString[per]],"HypothesisTestData",Method->Automatic,SignificanceLevel->sigLev];
(Subscript[ℒ, 2])=LocationEquivalenceTest[{GL1[[rdTime]],GL2[[rdTime]]},"HypothesisTestData",Method->Automatic,SignificanceLevel->sigLev];
Grid[{
If[sampNo==4,
{BoxWhiskerChart[Thread[Style[ToExpression["GL"<>ToString[per]],paleColor]],bxWhStyle,BarSpacing->2.5,ImageSize->imSizD,AspectRatio->.3],SpanFromLeft},
{BoxWhiskerChart[{Style[GL1[[rdTime]],paleColor[[rdTime]]],Style[GL2[[rdTime]],paleColor[[rdTime]]]},bxWhStyle,ChartLabels->{Style[First[tabHd][[1]],12,stdColor[[rdTime]]],Style[First[tabHd][[2]],12,stdColor[[rdTime]]]},BarSpacing->3,ImageSize->imSizD,AspectRatio->.3],SpanFromLeft}
], (* </If> *)
If[sampNo==4,
{Framed[Quiet@(Subscript[ℒ, 4])["TestDataTable",
locEqTests],FrameMargins->frMarg],Quiet@Show[BarChart[
Reverse[(Subscript[ℒ, 4])["PValue",locEqTests]],ChartLabels->Placed[Reverse[locEqTests],Left],
ChartStyle->fadedColor[[4]],AspectRatio->.25,ImageSize->360],Graphics[{Dashed,stdBlack,Line[{{sigLev,0},{sigLev,Length[locEqTests]+1}}]}]]},
{Framed[Quiet@(Subscript[ℒ, 2])["TestDataTable",
locEqShTests],FrameMargins->2.2frMarg],Quiet@Show[BarChart[Reverse[Quiet@(Subscript[ℒ, 2])["PValue",
locEqShTests]],ChartLabels->Placed[Reverse[locEqShTests],Left],ChartStyle->paleColor[[rdTime]],AspectRatio->.18,ImageSize->360],Graphics[{Dashed,stdBlack,Line[{{sigLev,0},{sigLev,Length[locEqShTests]+1}}]}]]}
],(* </If> *)
If[sampNo==4,
{Panel[TraditionalForm@Style[(Subscript[ℒ, 4])["TestConclusion",Automatic],14],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD],SpanFromLeft},
{Panel[TraditionalForm@Style[(Subscript[ℒ, 2])["TestConclusion",Automatic],
14],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD],SpanFromLeft}
],(* </If> *)
If[sampNo==4,
{Text[Style["Post-hoc analysis using a Bonferroni correction:",15]],SpanFromLeft},{Text[Style["No post-hoc analysis required.",15]],SpanFromLeft}],
If[sampNo==4,
{Panel[TraditionalForm@Style[TableForm@Transpose@{Subsets[Table[k,{k,4}],{2}],
Quiet@Thread[LocationEquivalenceTest[Subsets[Table[ToExpression["GL"<>ToString[per]][[k]],{k,4}],{2}],"TestConclusion",SignificanceLevel->sigLev/Binomial[Length[ToExpression["GL"<>ToString[per]]],2]]]}/.{1->Style[PB,stdColor[[1]]],2->Style[PL,stdColor[[2]]],3->Style[PD,stdColor[[3]]],4->Style[PS,stdColor[[4]]]},"Text",FontSize->10],FrameMargins->frMarg,ImageMargins->Automatic,BaseStyle->stdBlack,Alignment->{Left,Center},ImageSize->imSizD],SpanFromLeft
},{"",SpanFromLeft}]},
Alignment->{{Left,Left},{Bottom,Top}}],
Row[{Control[{{sampNo,4,"me(di)ans across"},{4->"read times",2->"periods"}}],
Control[{{per,1,perStr},{1->"1",2->"2"}}],
Control[{{rdTime,1,rdtStr},{1->pb,2->pl,3->pd,4->ps}}],
Control[{{sigLev,sigLevInit,sigLevStr},sigLevList}]},Spacer[6]]]

Note

Observe that the Student’s t-test (variant for unknown but same population variances) provides essetially the same results as above:

TTest/@Table[{GL1[[k]],GL2[[k]]},{k,4}]

Conclusions

The results of this statistical research can be summarized as follows:

The four-sample hypotheses

The null hypothesis that medians of PB, PL, PD, and PS glucose levels are equal at the 1%, 5% and 10% significance level for both Period 1 and Period 2 is rejected. The post-hoc analysis shows that:

In the Period 1, there is a difference between medians of all pairs with the exception of the PL–PD (at the 1% and 5% level) and PL–PS pairs (at the 1%, 5% and 10% level), i.e., the glucose level decreases from breakfast to dinner at the 10% significance level.

In the Period 2, there is a difference between medians of three pairs containing the PD sample at the 1%, 5% and 10% level, for the remaining three pairs, the pairwise null hypothesis cannot be rejected at any significance level used. In other words, particularly the PD glucose level is significantly less than the others, and is responsible for rejection of four-sample null hypothesis (1).

The two-sample hypotheses

The null hypothesis that mean and median of glucose level between periods is rejected (at 5% and 10% level) for PL and PS samples (at 1%, only for the PS samples). Using the results of descriptive statistical analysis, these findings can be interpreted as follows: from Period 1 to Period 2, the mean and median of the PL and PS samples significantly increased, while the remaining did not.

Medical guidelines

Interpretation of the statistical analysis is left only to doctors. From the long-term medical point of view, the level of HbA1C (glycated hemoglobin) is the decisive quantity for diabetes compensation. The present contribution is mentioned for educational purposes only.

These observations can only be interpreted as an indication that the behavior of the glucose level during the Period 2 may be evaluated statistically markedly worse than during Period 1.

Perspectives

We expect to extend the research by adding new data and other statistical tools, particularly regression analysis, time series analysis, and spectral analysis, which could potentially reveal some unknown facts in time dependence of glucose level.

Acknowledgment

The first author (I. B.) thanks to the parents of the child patient (her nephew) for providing his glucose level records.

Thank you for attention

References

Daly, M. E. (1998), Acute effects on insulin sensitivity and diurnal metabolic profiles of a high-sucrose compared with a high starch diet, Am. J. Clin. Nutr., 67, pages 1186–1196.
http://www.ajcn.org/content/67/6/1186.full.pdf

Ganong, W. F. (2005), Přehled lékařské fyziologie, 20. vydání, Galén, Praha, 890 pages, ISBN 80-7262-311-7. (In Czech)
http://www.lekarskeknihy.cz/?id=80-7262-311-7

Kahn, C. R., Weir, G. C., King, G. L., Moses, A. C., Smith, R. J. & Jacobson, A. M. (2004), Joslin’s Diabetes Mellitus, 14th edition, Lippincott Williams & Wilkins, ISBN 978-0781727969.
https://www.joslin.org/jstore/joslin_diabetes_mellitus_14th_edition.html

Litschmannová, M. (2011), Úvod do statistiky, VŠB-TU Ostrava, Fakulta elektrotechniky a informatiky, Ostrava. (In Czech)
http://mi21.vsb.cz/modul/uvod-do-statistiky

Mangano, S. (2010), Mathematica Cookbook. Building Blocks for Science, Engineering, Finance, Music, and More, 1st edition, O’Reilly, Sebastopol, xxiv+801 pages, ISBN 978-0-596-52099-1.
http://www.oreilly.com/

McKillup, S. (2006), Statistics Explained. An Introductory Guide for Life Scientists, Cambridge University Press, Cambridge, ISBN 0-521-54316-9.
http://www.cambridge.org/9780521835503

Nelson, D. (2004), The Penguin Dictionary of Statistics, Penguin Reference Series, Penguin Books, London, ISBN 0-140-51496-1.
http://www.penguin.com/

Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007), Numerical Recipes: The Art of Scientific Computing, 3rd edition, Cambridge University Press, Cambridge, 1256 pages, ISBN 0-521-88068-8.
http://www.nr.com/

Ruskeepää, H (2009), Mathematica Navigator. Mathematics, Statistics, and Graphics, 3rd edition, Academic Press (an imprint of Elsevier), Amsterdam, 1111 pages, ISBN 978-0-12-374164-6.
http://elsevier.com/