Research Statistics Central

 Problem:   You are running a standard quadratic (polynomial) regression analysis, and are specifically interested in the X and Y values at the peak.  If you use standard regression software, typically there will be no option that allows the peak to be estimated, with standard errors. Example:   You are studying Growth as a function of Age.  Of particular interest is when maximum Growth occurs, and at what Age. SAS code to generate artificial data, and run the analysis is: data one; do Age=1 to 20; Growth=95 + 2.7*Age - .3*Age*Age + 5*rannor(22); end; proc nlin plots=fit; parms int=2 lin=1 quad=1; model Growth = int + lin*Age + quad*Age*Age; estimate 'Age at peak' -lin/(2*quad); estimate 'Growth at peak' int + lin*(-lin/(2*quad)) + quad*(-lin/(2*quad))*(-lin/(2*quad)); run; The standard quadratic regression model with intercept, linear and quadratic slopes, is coded into Proc NLIN which has the ability to estimate any fun...

Factorial treatment designs are popular, due to advantages of research on multiple treatment factors and how they interact.  But if the design includes a control treatment that is not part of the factorial, problems occur in estimation of least squares means.  A typical example is shown here, with 2 fertilizer and 3 irrigation treatments, giving 6 factorial treatment combinations, plus a control that is defined by a 3rd level of fertilizer, and a 4th level of irrigation: Fert1:Irrig2               Fert2:Irrig1             Fert1:Irrig1 Fert2:Irrig3               Fert1:Irrig3             Fert2:Irrig2           Control Other situations might have the control sharing a level of one of the factors, for example the control might be defined as Fert2:Irrig4.  But this still causes problems with estim...

Objective :  We are fitting nonlinear regression lines to data, but have multiple groups (treatments), each with its own line.  Since groups are a factor of interest, particularly in how they change the lines (parameters of the model), we want to compare parameter estimates among the groups. First approach is to fit the nonlinear model to each group separately, then compare the parameter estimates using t-tests.  The code below generates a random example dataset, with 8 replicates for each of 5 treatments, all measured over 12 days.  Then Proc Nlmixed is used to fit the model explaining change in prate with water changes over the days, "by treat", and parameter estimates are output to data ppp.  This ppp dataset is processed to collect the estimates and standard errors, and t-tests are calculated for all 5*(5-1)/2 comparisons.  Code will need to be customized for new data, including number and values of treatments, degrees of freedom, and parameter names...

No guarantees, but roughly if means differ by 3*SE then they are statistically significant.   This is based on the Least Significant Difference, which is 2*sqrt(2)*SE.  Often people use non-overlapping confidence intervals as a decision rule, but this is equivalent to 4*SE, which is a bit conservative. Things that make 3*SE fail: 1)   Actually statistical differences depend on the standard error of difference, SED, not SE.   Anything in the model that makes these differ will make the rule fail, such as covariates and blocking factors. 2)   In general, mixed models with random effects will make the rule fail, because random variance is included in SE, but not in SED.   But this will make 3*SE rule conservative, 3*SED will be even smaller.   If 3*SE suggests a statistical difference, difference most likely exists. Also take a look at Error Bars paper .

Objective and Background :  You have text data that is UTF encoded and need SAS/R to read and write datasets with that encoding.  If you have ever printed or viewed text information, and seen something like Giuffr?Ÿ’e?ƒe?Ÿƒ?ÿ?›ƒ?ªƒ?›?Ÿ’e›ƒ?ª­?Ÿƒeee, then you are running into this encoding issue.  Computers store text using numbers, with each number assigned to a particular character.  See  https://en.wikipedia.org/wiki/ASCII  to find that the character & is stored as 38 when using the ASCII encoding.  Unicode is popular internationally because it encodes special characters such as accented letters, and UTF-8 is a widely used version ( https://en.wikipedia.org/wiki/UTF-8 ).  In UTF-8 the & character is stored as 26, and you can imagine how the jumbled example above arises from the confusion of what letters are being stored. Solution 1 :  Use options to request that individual datasets be read and written in a particular encodin...

Objective :  Transformed data, for example log(y), is analyzed to correct normality or equal variance requirements.  But we want to report means and standard errors in the original units. SAS example : data one;  do treat=1 to 3;  do rep=1 to 5;    y=10 + treat+ exp(rannor(111));    logy=log(y);    output;  end;end; run; proc mixed plots=all;   class treat;   model y=treat;   lsmeans treat/pdiff; run; proc mixed plots=all;   class treat;   model logy=treat;   lsmeans treat/pdiff; run; The original data, variable y, might have units of pounds.  If a transformation is needed, we simply calculate a new variable by applying a mathematical function known to improve normality or equal variance, and run the same analysis on the new variable.  Commonly used choices are listed in the second table below. However, looking at the results for both analyses we see treat Mean Y S...

Objective : I am running a statistical analysis in SAS, and the default ODS graphics look good, but I need them to be publication quality. SAS can automatically create some nice graphs, and has greatly increased the availability of graphs within procedures.  If you like what you see, you might copy graphs directly from the SAS output window, or possibly you save graphs and output to a pdf or other external file format.  But this output will be low quality, generally 75 dpi.  Instead, add the following statements to write graphics directly to files, allowing control of format and quality. ods graphics on /       width=7in       imagefmt=tiff       imagemap=off       imagename="MyPlot"       border=off; ods listing file="Body.rtf" style=journal gpath="."  dpi=600; Once these statements have been submitted, all graphs created by subsequent procedures will be written  to files na...

Objective :  You are running ANOVAs or regressions in SAS, and wish there was a way to avoid writing the dozens of commands needed to conduct the analysis and generate recommended diagnostics and summary of results, not to mention the hundreds of possible options that might be needed to access recommended methods.  A possible solution is to download a copy of danda.sas below, and use this macro collection to run the dozens of commands with one statement.  We will also have future posts covering various uses of danda.sas, giving examples as always. danda.sas is under continued development, check this page for updates. Date                       Version               Link 2021/03/15             2.12.030          danda.sas 2021/03/15       ...

Objective :  You are running a linear model, for example ANOVA or regression, and are using the "ANOVA table" to decide which terms in the model are influencing the dependent variable.  You check the numerator and denominator degrees of freedom, as recommended , to guard against modeling errors and use of wrong error terms.  Reported values disagree with what you expected, so now what? Make sure your expected numbers are calculated correctly : An example model is (last term is the residual error) [Model 1]       y = u + block + treat + block*treat + rep(block*treat) If numbers of levels are b=2 for blocks, t=2 for treats, and r=5 for reps, then we expect degrees of freedom (DF) to be DF[block] = b-1 = 1 DF[treat] = t-1 = 1 DF[block*treat] = (b-1)*(t-1) = 1*1 = 1 DF[rep(block*treat)] = (r-1)*(b)*(t) = 4*2*2 = 16 Number of observations is b*t*r = 20, and DF add to 19, which is 20 minus the one DF for the intercept, as expected. DF rules are: ...

Objective :  Individuals are measured for several characteristics, and we want to group the individuals based on similarities across all characteristics. Statistical options are cluster analysis, principle components (PCA), and biplots.  PCA has the advantage of combining correlated variables together, reducing the complexity of explaining why individuals cluster together.  And biplots adds some nice features to PCA.  This post compares these choices. Example :  10 farms measured for nutrients in grass, and the primary question is to see if/which farms have similar nutrient profiles. Create a random dataset with 4 nutrients measured on 4 pastures in each of 10 farms. Run the SAS code for producing biplots.  Here we restrict the number of PC to 2 (n=2), in general you would use PCA to decide how many components are needed.  Prinqual requires all variables to be processed by a Transform statement, here we use the identity transformation so the variab...

Objective : We are comparing means using ANOVA, and our treatment levels are amounts of something.  Thus regression hypotheses may shed light on how the treatments differ, for example is there an overall linear trend for the response variable to increase or decrease with treatment level.  This is addressed by adding orthogonal polynomial contrasts to our ANOVA, which may require that we add contrast coefficients. Example :  Treatments are amounts of corn in the diet, specifically 62%, 65%, 68%, 71% and 74%. SAS :  IML product has an orthogonal polynomial calculator.  Additional code here attempts to make the coefficients whole numbers by dividing by the smallest non-zero number.  Note IML may not be available, depending on your license. proc iml; trtlevels={0.62, 0.65,0.68,0.71,0.74}; **this is only user input; ntrt=nrow(trtlevels); coeff=orpol(trtlevels); coeff = coeff[,2:ntrt]; div=abs(coeff); zerloc=loc(div<1e-14); if n...

Objective :  We need to know how many observations to collect so our estimate of the mean has a useful precision.  For example, how many animals should be measured in order to have an 80% chance that the 95% confidence interval for weight will be no wider than 20 kg?  In addition to those 3 numbers, we also need an estimate of the std. deviation.  Suppose the best situation expected has SD=20kg, but we also want to see what changes if SD=40kg, SAS:   Run this code proc power;    onesamplemeans ci=t       alpha = 0.05       halfwidth = 10       stddev = 20 40       probwidth = 0.80       ntotal = .; run; data adjust;  samplesize=22;  population=1200;  adjsamplesize=ceil(samplesize/(1 + ((samplesize-1)/population))); run; proc print; run; The 95% confidence interval is requested by setting alpha=0.05. The 80% chance that our experiment will satisf...

 In our role as statisticians, we see many identical questions from researchers on how to perform common statistical methods.  Rather than repeatedly give the same answers, we plan to post the current recommended practice for each type of problem, and update as needed.  Our focus will be on R and SAS software, widely used and what we have the most experience with.  Like most statistical consultants, we are educators, hoping to teach you how to think about and do statistics.    Statistics is a broad topic, reaching into almost every area of science.  Since the different sciences have different types of data, different sources of variation, statistical methodology can be area specific.  The questions we see mainly come from agriculture, so there will be a bias towards methods most widely used there. We will strive to place sufficient key words in the posts to enable searches to uncover appropriate material.    Time permitting, we will p...