| Title: | 'SAS' Linear Model |
|---|---|
| Description: | This is a core implementation of 'SAS' procedures for linear models - GLM, REG, ANOVA, TTEST, FREQ, and UNIVARIATE. Some R packages provide Type II and Type III SS. However, the results of nested and complex designs are often different from those of 'SAS'. Different results do not necessarily mean incorrectness. However, many want the same results as 'SAS'. This package aims to achieve that. Reference: Littell RC, Stroup WW, Freund RJ (2002, ISBN:0-471-22174-0). |
| Authors: | Kyun-Seop Bae [aut, cre] |
| Maintainer: | Kyun-Seop Bae <[email protected]> |
| License: | GPL-3 |
| Version: | 1.0.1 |
| Built: | 2026-07-14 17:53:26 UTC |
| Source: | https://github.com/cran/sasLM |
This is a core implementation of 'SAS' procedures for linear models - GLM, REG, and ANOVA. Some packages provide type II and type III SS. However, the results of nested and complex designs are often different from those of 'SAS'. A different result does not necessarily mean incorrectness. However, many want the same results as 'SAS'. This package aims to achieve that. Reference: Littell RC, Stroup WW, Freund RJ (2002, ISBN:0-471-22174-0).
This will serve those who want SAS PROC GLM, REG, and ANOVA in R.
Kyun-Seop Bae [email protected]
## SAS PROC GLM Script for Typical Bioequivalence Data # PROC GLM DATA=BEdata; # CLASS SEQ SUBJ PRD TRT; # MODEL LNCMAX = SEQ SUBJ(SEQ) PRD TRT; # RANDOM SUBJ(SEQ)/TEST; # LSMEANS TRT / DIFF=CONTROL("R") CL ALPHA=0.1; # ODS OUTPUT LSMeanDiffCL=LSMD; # DATA LSMD; SET LSMD; # PE = EXP(DIFFERENCE); # LL = EXP(LowerCL); # UL = EXP(UpperCL); # PROC PRINT DATA=LSMD; RUN; ## ## SAS PROC GLM equivalent BEdata = af(BEdata, c("SEQ", "SUBJ", "PRD", "TRT")) # Columns as factor formula1 = log(CMAX) ~ SEQ/SUBJ + PRD + TRT # Model GLM(formula1, BEdata) # ANOVA tables of Type I, II, III SS RanTest(formula1, BEdata, Random="SUBJ") # Hypothesis test with SUBJ as random ci0 = CIest(formula1, BEdata, "TRT", c(-1, 1), 0.90) # 90% CI exp(ci0[, c("Estimate", "Lower CL", "Upper CL")]) # 90% CI of GMR ## 'nlme' or SAS PROC MIXED is preferred for an unbalanced case ## SAS PROC MIXED equivalent # require(nlme) # Result = lme(log(CMAX) ~ SEQ + PRD + TRT, random=~1|SUBJ, data=BEdata) # summary(Result) # VarCorr(Result) # ci = intervals(Result, 0.90) ; ci # exp(ci$fixed["TRTT",]) #### SAS PROC GLM Script for Typical Bioequivalence Data # PROC GLM DATA=BEdata; # CLASS SEQ SUBJ PRD TRT; # MODEL LNCMAX = SEQ SUBJ(SEQ) PRD TRT; # RANDOM SUBJ(SEQ)/TEST; # LSMEANS TRT / DIFF=CONTROL("R") CL ALPHA=0.1; # ODS OUTPUT LSMeanDiffCL=LSMD; # DATA LSMD; SET LSMD; # PE = EXP(DIFFERENCE); # LL = EXP(LowerCL); # UL = EXP(UpperCL); # PROC PRINT DATA=LSMD; RUN; ## ## SAS PROC GLM equivalent BEdata = af(BEdata, c("SEQ", "SUBJ", "PRD", "TRT")) # Columns as factor formula1 = log(CMAX) ~ SEQ/SUBJ + PRD + TRT # Model GLM(formula1, BEdata) # ANOVA tables of Type I, II, III SS RanTest(formula1, BEdata, Random="SUBJ") # Hypothesis test with SUBJ as random ci0 = CIest(formula1, BEdata, "TRT", c(-1, 1), 0.90) # 90% CI exp(ci0[, c("Estimate", "Lower CL", "Upper CL")]) # 90% CI of GMR ## 'nlme' or SAS PROC MIXED is preferred for an unbalanced case ## SAS PROC MIXED equivalent # require(nlme) # Result = lme(log(CMAX) ~ SEQ + PRD + TRT, random=~1|SUBJ, data=BEdata) # summary(Result) # VarCorr(Result) # ci = intervals(Result, 0.90) ; ci # exp(ci$fixed["TRTT",]) ##
Insert empty digital-signature form fields (AcroForm /Sig fields) into a PDF
so it can be signed in Adobe Acrobat Reader (the free reader) with one click: open
the PDF, click a field, and sign with a Digital ID. The field is added with a
pure base-R PDF incremental update, with no external tools or packages, so it works
wherever R does. It is designed for the simple PDFs produced by IQLM
and OQLM.
Acrobat Reader can also sign a PDF that has no field at all (“All tools” > “Use a certificate” > “Digitally sign”, then drag a rectangle). This function only makes the workflow click-to-sign by pre-placing labelled fields.
addSigFieldLM(pdf, out = pdf, page = 1L, fieldNames = c("Performed_by", "Reviewed_by"), rects = NULL)addSigFieldLM(pdf, out = pdf, page = 1L, fieldNames = c("Performed_by", "Reviewed_by"), rects = NULL)
pdf |
path to the input PDF. |
out |
path to write the result; defaults to overwriting |
page |
1-based page number to place the field(s) on (the signature page is
page 1 of |
fieldNames |
character vector of field names; one signature field is added per name. The names appear in Acrobat's Signature panel. |
rects |
optional list of numeric length-4 rectangles |
The function appends an incremental-update section that adds an /AcroForm to
the document catalog (with /SigFlags 3), one /Widget signature
annotation per field, and the field references to the page's /Annots. The
original content is left untouched. The resulting fields are unsigned; the actual
cryptographic signature is applied by Acrobat Reader using the signer's Digital ID.
Only classic cross-reference-table PDFs are supported (those written by R's
pdf device); the function stops on cross-reference-stream
PDFs.
Invisibly, the output path. A PDF with signature fields is written as a side effect.
Kyun-Seop Bae <[email protected]>
#IQLM("sasLM-IQ-Report.pdf", performedBy = "Kyun-Seop Bae") #addSigFieldLM("sasLM-IQ-Report.pdf") # adds Performed_by / Reviewed_by fields ## or in one step: #IQLM("sasLM-IQ-Report.pdf", performedBy = "Kyun-Seop Bae", sigField = TRUE) ## then open in Acrobat Reader and click each field to sign with your Digital ID.#IQLM("sasLM-IQ-Report.pdf", performedBy = "Kyun-Seop Bae") #addSigFieldLM("sasLM-IQ-Report.pdf") # adds Performed_by / Reviewed_by fields ## or in one step: #IQLM("sasLM-IQ-Report.pdf", performedBy = "Kyun-Seop Bae", sigField = TRUE) ## then open in Acrobat Reader and click each field to sign with your Digital ID.
Conveniently convert some columns of a data.frame into factors.
af(DataFrame, Cols)af(DataFrame, Cols)
DataFrame |
a |
Cols |
column names or indices to be converted |
It performs conversion of some columns in a data.frame into factors conveniently.
Returns a data.frame with converted columns.
Kyun-Seop Bae [email protected]
ANOVA with Type I SS.
aov1(Formula, Data, BETA=FALSE, Resid=FALSE)aov1(Formula, Data, BETA=FALSE, Resid=FALSE)
Formula |
a conventional formula for a linear model. |
Data |
a |
BETA |
if |
Resid |
if |
It performs the core function of SAS PROC GLM, and returns Type I SS. This also accepts continuous independent variables.
The result table is comparable to that of SAS PROC ANOVA.
Df |
degrees of freedom |
Sum Sq |
sum of squares for the set of contrasts |
Mean Sq |
mean square |
F value |
F value for the F distribution |
Pr(>F) |
probability of a larger F value |
The next returns are optional.
Parameter |
Parameter table with standard error, t value, p value. |
Fitted |
Fitted values or y hat. This is returned only with the |
Residual |
Weighted residuals. This is returned only with the |
Kyun-Seop Bae [email protected]
aov1(uptake ~ Plant + Type + Treatment + conc, CO2) aov1(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE) aov1(uptake ~ Plant + Type + Treatment + conc, CO2, Resid=TRUE) aov1(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE, Resid=TRUE)aov1(uptake ~ Plant + Type + Treatment + conc, CO2) aov1(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE) aov1(uptake ~ Plant + Type + Treatment + conc, CO2, Resid=TRUE) aov1(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE, Resid=TRUE)
ANOVA with Type II SS.
aov2(Formula, Data, BETA=FALSE, Resid=FALSE)aov2(Formula, Data, BETA=FALSE, Resid=FALSE)
Formula |
a conventional formula for a linear model. |
Data |
a |
BETA |
if |
Resid |
if |
It performs the core function of SAS PROC GLM, and returns Type II SS. This also accepts continuous independent variables.
The result table is comparable to that of SAS PROC ANOVA.
Df |
degrees of freedom |
Sum Sq |
sum of squares for the set of contrasts |
Mean Sq |
mean square |
F value |
F value for the F distribution |
Pr(>F) |
probability of a larger F value |
The next returns are optional.
Parameter |
Parameter table with standard error, t value, p value. |
Fitted |
Fitted values or y hat. This is returned only with the |
Residual |
Weighted residuals. This is returned only with the |
Kyun-Seop Bae [email protected]
aov2(uptake ~ Plant + Type + Treatment + conc, CO2) aov2(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE) aov2(uptake ~ Plant + Type + Treatment + conc, CO2, Resid=TRUE) aov2(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE, Resid=TRUE) aov2(uptake ~ Type, CO2) aov2(uptake ~ Type - 1, CO2)aov2(uptake ~ Plant + Type + Treatment + conc, CO2) aov2(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE) aov2(uptake ~ Plant + Type + Treatment + conc, CO2, Resid=TRUE) aov2(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE, Resid=TRUE) aov2(uptake ~ Type, CO2) aov2(uptake ~ Type - 1, CO2)
ANOVA with Type III SS.
aov3(Formula, Data, BETA=FALSE, Resid=FALSE)aov3(Formula, Data, BETA=FALSE, Resid=FALSE)
Formula |
a conventional formula for a linear model. |
Data |
a |
BETA |
if |
Resid |
if |
It performs the core function of SAS PROC GLM, and returns Type III SS. This also accepts continuous independent variables.
The result table is comparable to that of SAS PROC ANOVA.
Df |
degrees of freedom |
Sum Sq |
sum of squares for the set of contrasts |
Mean Sq |
mean square |
F value |
F value for the F distribution |
Pr(>F) |
probability of a larger F value |
The next returns are optional.
Parameter |
Parameter table with standard error, t value, p value. |
Fitted |
Fitted values or y hat. This is returned only with the |
Residual |
Weighted residuals. This is returned only with the |
Kyun-Seop Bae [email protected]
aov3(uptake ~ Plant + Type + Treatment + conc, CO2) aov3(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE) aov3(uptake ~ Plant + Type + Treatment + conc, CO2, Resid=TRUE) aov3(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE, Resid=TRUE)aov3(uptake ~ Plant + Type + Treatment + conc, CO2) aov3(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE) aov3(uptake ~ Plant + Type + Treatment + conc, CO2, Resid=TRUE) aov3(uptake ~ Plant + Type + Treatment + conc, CO2, BETA=TRUE, Resid=TRUE)
The data is from 'Canner PL. An overview of six clinical trials of aspirin in coronary heart disease. Stat Med. 1987'
aspirinCHDaspirinCHD
A data frame with 6 rows.
y1death event count of aspirin group
n1total subjects of the aspirin group
y2death event count of placebo group
n2total subjects of the placebo group
This data is for educational purposes.
Canner PL. An overview of six clinical trials of aspirin in coronary heart disease. Stat Med. 1987;6:255-263.
Contains Cmax data from a real bioequivalence study.
BEdataBEdata
A data frame with 91 observations on the following 6 variables.
ADMAdmission or Hospitalization Group Code: 1, 2, or 3
SEQGroup or Sequence character code: 'RT' or 'TR'
PRDPeriod numeric value: 1 or 2
TRTTreatment or Drug code: 'R' or 'T'
SUBJSubject ID
CMAXCmax values
This contains real data from a 2x2 bioequivalence study, which has three different hospitalization groups. See Bae KS, Kang SH. Bioequivalence data analysis for the case of separate hospitalization. Transl Clin Pharmacol. 2017;25(2):93-100. doi.org/10.12793/tcp.2017.25.2.93
Trailing zeros after integers are somewhat annoying. This removes them from the vector of strings.
bk(ktab, rpltag=c("n", "N"), dig=10)bk(ktab, rpltag=c("n", "N"), dig=10)
ktab |
an output of |
rpltag |
tag string of replacement rows. This is usually "n", which means the sample count. |
dig |
maximum digits of decimals in the |
This is convenient if used with tsum0, tsum1, tsum2, or tsum3. This requires knitr::kable.
A new processed vector of strings. The class is still knitr_kable.
Kyun-Seop Bae [email protected]
## OUTPUT example # t0 = tsum0(CO2, "uptake", c("mean", "median", "sd", "length", "min", "max")) # bk(kable(t0)) # requires knitr package # # | | x| # |:------|--------:| # |mean | 27.21310| # |median | 28.30000| # |sd | 10.81441| # |n | 84 | # |min | 7.70000| # |max | 45.50000| # t1 = tsum(uptake ~ Treatment, CO2, # e=c("mean", "median", "sd", "min", "max", "length"), # ou=c("chilled", "nonchilled"), # repl=list(c("median", "length"), c("med", "N"))) # # bk(kable(t1, digits=3)) # requires knitr package # # | | chilled| nonchilled| Combined| # |:----|-------:|----------:|--------:| # |mean | 23.783| 30.643| 27.213| # |med | 19.700| 31.300| 28.300| # |sd | 10.884| 9.705| 10.814| # |min | 7.700| 10.600| 7.700| # |max | 42.400| 45.500| 45.500| # |N | 42 | 42 | 84 |## OUTPUT example # t0 = tsum0(CO2, "uptake", c("mean", "median", "sd", "length", "min", "max")) # bk(kable(t0)) # requires knitr package # # | | x| # |:------|--------:| # |mean | 27.21310| # |median | 28.30000| # |sd | 10.81441| # |n | 84 | # |min | 7.70000| # |max | 45.50000| # t1 = tsum(uptake ~ Treatment, CO2, # e=c("mean", "median", "sd", "min", "max", "length"), # ou=c("chilled", "nonchilled"), # repl=list(c("median", "length"), c("med", "N"))) # # bk(kable(t1, digits=3)) # requires knitr package # # | | chilled| nonchilled| Combined| # |:----|-------:|----------:|--------:| # |mean | 23.783| 30.643| 27.213| # |med | 19.700| 31.300| 28.300| # |sd | 10.884| 9.705| 10.814| # |min | 7.700| 10.600| 7.700| # |max | 42.400| 45.500| 45.500| # |N | 42 | 42 | 84 |
Functions such as GLM, REG, and aov1 can be run by levels of a variable.
BY(FUN, Formula, Data, By, ...)BY(FUN, Formula, Data, By, ...)
FUN |
Function name to be called, such as GLM or REG |
Formula |
a conventional formula for a linear model. |
Data |
a |
By |
a variable name in the |
... |
arguments to be passed to |
This mimics the BY clause of SAS procedures.
a list of FUN function outputs. The names of the list are the levels of the By variable.
Kyun-Seop Bae [email protected]
BY(GLM, uptake ~ Treatment + as.factor(conc), CO2, By="Type") BY(REG, uptake ~ conc, CO2, By="Type")BY(GLM, uptake ~ Treatment + as.factor(conc), CO2, By="Type") BY(REG, uptake ~ conc, CO2, By="Type")
Get the point estimate and its confidence interval with a given contrast and alpha value using the t distribution.
CIest(Formula, Data, Term, Contrast, conf.level=0.95)CIest(Formula, Data, Term, Contrast, conf.level=0.95)
Formula |
a conventional formula for a linear model |
Data |
a |
Term |
a factor name to be estimated |
Contrast |
a level vector. Levels are alphabetically ordered by default. |
conf.level |
confidence level of the confidence interval |
Get the point estimate and its confidence interval with a given contrast and alpha value using the t distribution.
Estimate |
point estimate of the input linear contrast |
Lower CL |
lower confidence limit |
Upper CL |
upper confidence limit |
Std. Error |
standard error of the point estimate |
t value |
value for the t distribution |
Df |
degrees of freedom |
Pr(>|t|) |
probability of a larger absolute t value from the t distribution with the residual degrees of freedom |
Kyun-Seop Bae [email protected]
CIest(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, "TRT", c(-1, 1), 0.90) # 90% CICIest(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, "TRT", c(-1, 1), 0.90) # 90% CI
Collinearity diagnostics with tolerance, VIF, eigenvalue, condition index, and variance proportions
Coll(Formula, Data)Coll(Formula, Data)
Formula |
formula of the model |
Data |
input data as a matrix or a |
Sometimes collinearity diagnostics after multiple linear regression are necessary.
Tol |
tolerance of independent variables |
VIF |
variance inflation factor of independent variables |
Eigenvalue |
eigenvalue of Z'Z (crossproduct) of standardized independent variables |
Cond. Index |
condition index |
Proportions of variances |
under the names of coefficients |
Kyun-Seop Bae [email protected]
Coll(mpg ~ disp + hp + drat + wt + qsec, mtcars)Coll(mpg ~ disp + hp + drat + wt + qsec, mtcars)
Do an F test with a given set of contrasts.
CONTR(L, Formula, Data, mu=0)CONTR(L, Formula, Data, mu=0)
L |
contrast matrix. Each row is a contrast. |
Formula |
a conventional formula for a linear model |
Data |
a |
mu |
a vector of mu for the hypothesis L. The length should be equal to the row count of L. |
It performs an F test with a given set of contrasts (a matrix). It is similar to the CONTRAST clause of SAS PROC GLM. This can test the hypothesis that the linear combination (function)'s mean vector is mu.
Returns the sum of squares and its F value and p-value.
Df |
degrees of freedom |
Sum Sq |
sum of squares for the set of contrasts |
Mean Sq |
mean square |
F value |
F value for the F distribution |
Pr(>F) |
probability of a larger F value |
Kyun-Seop Bae [email protected]
CONTR(t(c(0, -1, 1)), uptake ~ Type, CO2) # sum of square GLM(uptake ~ Type, CO2) # compare with the aboveCONTR(t(c(0, -1, 1)), uptake ~ Type, CO2) # sum of square GLM(uptake ~ Type, CO2) # compare with the above
Testing correlation between numeric columns of data with the Pearson method.
Cor.test(Data, conf.level=0.95)Cor.test(Data, conf.level=0.95)
Data |
a matrix or a |
conf.level |
confidence level |
It uses all numeric columns of the input data. It uses "pairwise.complete.obs" rows.
Row names show which columns are used for the test.
Estimate |
point estimate of correlation |
Lower CL |
lower confidence limit |
Upper CL |
upper confidence limit |
t value |
t value of the t distribution |
Df |
degrees of freedom |
Pr(>|t|) |
probability with the t distribution |
Kyun-Seop Bae [email protected]
Cor.test(mtcars)Cor.test(mtcars)
Testing correlation between two numeric vectors by Fisher's Z transformation.
corFisher(x, y, conf.level=0.95, rho=0)corFisher(x, y, conf.level=0.95, rho=0)
x |
the first input numeric vector |
y |
the second input numeric vector |
conf.level |
confidence level |
rho |
population correlation rho under the null hypothesis |
This accepts only two numeric vectors.
N |
sample size, length of input vectors |
r |
sample correlation |
Fisher.z |
Fisher's z |
bias |
bias to correct |
rho.hat |
point estimate of population rho |
conf.level |
confidence level for the confidence interval |
lower |
lower limit of confidence interval |
upper |
upper limit of confidence interval |
rho0 |
population correlation rho under the null hypothesis |
p.value |
p value under the null hypothesis |
Kyun-Seop Bae [email protected]
Fisher RA. Statistical Methods for Research Workers. 14e. 1973
corFisher(mtcars$disp, mtcars$hp, rho=0.6)corFisher(mtcars$disp, mtcars$hp, rho=0.6)
Calculates the sum of squares of a contrast from an lfit result.
cSS(K, rx, mu=0, eps=1e-8)cSS(K, rx, mu=0, eps=1e-8)
K |
contrast matrix. Each row is a contrast. |
rx |
a result of the |
mu |
a vector of mu for the hypothesis K. The length should be equal to the row count of K. |
eps |
Values less than this are considered zero. |
It calculates the sum of squares with a given contrast matrix and an lfit result. It corresponds to SAS PROC GLM CONTRAST. This can test the hypothesis that the linear combination (function)'s mean vector is mu.
Returns the sum of squares and its F value and p-value.
Df |
degrees of freedom |
Sum Sq |
sum of squares for the set of contrasts |
Mean Sq |
mean square |
F value |
F value for the F distribution |
Pr(>F) |
probability of a larger F value |
Kyun-Seop Bae [email protected]
rx = REG(uptake ~ Type, CO2, summarize=FALSE) cSS(t(c(0, -1, 1)), rx) # sum of square GLM(uptake ~ Type, CO2) # compare with the aboverx = REG(uptake ~ Type, CO2, summarize=FALSE) cSS(t(c(0, -1, 1)), rx) # sum of square GLM(uptake ~ Type, CO2) # compare with the above
Cumulative alpha values for repeated hypothesis tests with changing bound z-values and times of test (ti).
CumAlpha(z, side=2, ti=NULL, c0=NULL, Seed=5)CumAlpha(z, side=2, ti=NULL, c0=NULL, Seed=5)
z |
vector of upper z-value bounds for the repeated hypothesis test |
side |
1=one-sided test, 2=two-sided test |
ti |
vector of times (or information amount) of test. All values should be in [0, 1] and sorted. If not specified, equal intervals are assumed. |
c0 |
correlation matrix. If not specified, Brownian motion is assumed. |
Seed |
seed value for the |
It calculates cumulative alpha-values for the repeated hypothesis test with a vector of upper bound z-values. If the times of test are not specified, linear (proportional) increase of information amount and Brownian motion of z-values are assumed, i.e. the correlation is sqrt(t_i/t_j).
The result is a matrix.
ti |
time of test |
cum.alpha |
cumulative alpha values |
Kyun-Seop Bae [email protected]
Reboussin DM, DeMets DL, Kim K, Lan KKG. Computations for group sequential boundaries using the Lan-DeMets function method. Controlled Clinical Trials. 2000;21:190-207.
CumAlpha(z=rep(qnorm(1 - 0.05/2), 10)) # two-side Z-test with alpha=0.05 for ten timesCumAlpha(z=rep(qnorm(1 - 0.05/2), 10)) # two-side Z-test with alpha=0.05 for ten times
Coefficient of variation in percentage.
CV(y)CV(y)
y |
a numeric vector |
It removes NA.
Coefficient of variation in percentage.
Kyun-Seop Bae [email protected]
CV(mtcars$mpg)CV(mtcars$mpg)
Plot pairwise differences by a common method.
Diffogram(Formula, Data, Term, conf.level=0.95, adj="lsd", ...)Diffogram(Formula, Data, Term, conf.level=0.95, adj="lsd", ...)
Formula |
a conventional formula for a linear model |
Data |
a |
Term |
a factor name to be estimated |
conf.level |
confidence level of the confidence interval |
adj |
"lsd", "tukey", "scheffe", "bon", or "duncan" to adjust p-value and confidence limit |
... |
arguments to be passed to |
This usually shows the shortest interval. It corresponds to the PDIFF option of SAS PROC GLM. For the adjustment method "dunnett", see the PDIFF function.
no return value, but a plot on the current device
Kyun-Seop Bae [email protected]
Diffogram(uptake ~ Type*Treatment + as.factor(conc), CO2, "as.factor(conc)")Diffogram(uptake ~ Type*Treatment + as.factor(conc), CO2, "as.factor(conc)")
Calculate the drift value with given upper bounds (z-values), times of test, and power.
Drift(bi, ti=NULL, Power=0.9)Drift(bi, ti=NULL, Power=0.9)
bi |
upper bound z-values |
ti |
times of test. These should be in the range of [0, 1]. If omitted, equal intervals are assumed. |
Power |
target power at the final test |
It calculates the drift value with given upper bound z-values, times of test, and power. If the times of test are not given, equal intervals are assumed. mvtnorm::pmvt (with noncentrality) is better than pmvnorm in calculating power and sample size. However, Lan-DeMets used the multivariate normal rather than the multivariate noncentral t distribution. This function follows Lan-DeMets for consistency with previous results.
Drift value for the given condition
Kyun-Seop Bae [email protected]
Reboussin DM, DeMets DL, Kim K, Lan KKG. Computations for group sequential boundaries using the Lan-DeMets function method. Controlled Clinical Trials. 2000;21:190-207.
Drift(seqBound(ti=(1:5)/5)[, "up.bound"])Drift(seqBound(ti=(1:5)/5)[, "up.bound"])
Makes a contrast matrix for Type I SS using the forward Doolittle method.
e1(XpX, eps=1e-8)e1(XpX, eps=1e-8)
XpX |
the crossproduct of a design or model matrix. This should have appropriate column names. |
eps |
A value less than this is considered zero. |
It makes a contrast matrix for Type I SS. If zapsmall is used, the result becomes less accurate.
A contrast matrix for Type I SS.
Kyun-Seop Bae [email protected]
x = ModelMatrix(uptake ~ Plant + Type + Treatment + conc, CO2) round(e1(crossprod(x$X)), 12)x = ModelMatrix(uptake ~ Plant + Type + Treatment + conc, CO2) round(e1(crossprod(x$X)), 12)
Makes a contrast matrix for Type II SS.
e2(x, eps=1e-8)e2(x, eps=1e-8)
x |
an output of |
eps |
A value less than this is considered zero. |
It makes a contrast matrix for Type II SS. If zapsmall is used, the result becomes less accurate.
A contrast matrix for Type II SS.
Kyun-Seop Bae [email protected]
round(e2(ModelMatrix(uptake ~ Plant + Type + Treatment + conc, CO2)), 12) round(e2(ModelMatrix(uptake ~ Type, CO2)), 12) round(e2(ModelMatrix(uptake ~ Type - 1, CO2)), 12)round(e2(ModelMatrix(uptake ~ Plant + Type + Treatment + conc, CO2)), 12) round(e2(ModelMatrix(uptake ~ Type, CO2)), 12) round(e2(ModelMatrix(uptake ~ Type - 1, CO2)), 12)
Makes a contrast matrix for Type III SS.
e3(x, eps=1e-8)e3(x, eps=1e-8)
x |
an output of |
eps |
A value less than this is considered zero. |
It makes a contrast matrix for Type III SS. If zapsmall is used, the result becomes less accurate.
A contrast matrix for Type III SS.
Kyun-Seop Bae [email protected]
round(e3(ModelMatrix(uptake ~ Plant + Type + Treatment + conc, CO2)), 12)round(e3(ModelMatrix(uptake ~ Plant + Type + Treatment + conc, CO2)), 12)
Calculates a formula table for the expected mean square of the given contrast. The default is for Type III SS.
EMS(Formula, Data, Type=3, eps=1e-8)EMS(Formula, Data, Type=3, eps=1e-8)
Formula |
a conventional formula for a linear model |
Data |
a |
Type |
type of sum of squares. The default is 3. Type 4 is not supported yet. |
eps |
A value less than this is considered zero. |
This is necessary for further hypothesis tests of nesting factors.
A coefficient matrix for Type III expected mean square
Kyun-Seop Bae [email protected]
f1 = log(CMAX) ~ SEQ/SUBJ + PRD + TRT EMS(f1, BEdata) EMS(f1, BEdata, Type=1) EMS(f1, BEdata, Type=2)f1 = log(CMAX) ~ SEQ/SUBJ + PRD + TRT EMS(f1, BEdata) EMS(f1, BEdata, Type=1) EMS(f1, BEdata, Type=2)
Estimates Linear Functions with a given GLM result.
est(L, X, rx, conf.level=0.95, adj="lsd", paired=FALSE)est(L, X, rx, conf.level=0.95, adj="lsd", paired=FALSE)
L |
a matrix of linear contrast rows to be tested |
X |
a model (design) matrix from |
rx |
a result of the |
conf.level |
confidence level of the confidence limit |
adj |
adjustment method for grouping. This supports "tukey", "bon", "scheffe", "duncan", and "dunnett". This only affects grouping, not the confidence interval. |
paired |
If this is |
It tests rows of linear functions. A linear function means a linear combination of estimated coefficients. It corresponds to the ESTIMATE statement of SAS PROC GLM. The same sample size per group is assumed for the Tukey adjustment.
Estimate |
point estimate of the input linear contrast |
Lower CL |
lower confidence limit by the "lsd" method |
Upper CL |
upper confidence limit by the "lsd" method |
Std. Error |
standard error of the point estimate |
t value |
value for the t distribution, for methods other than "scheffe" |
F value |
value for the F distribution, for the "scheffe" method only |
Df |
degrees of freedom of the residuals |
Pr(>|t|) |
probability of a larger absolute t value from the t distribution with the residual degrees of freedom, for methods other than "scheffe" |
Pr(>F) |
probability of a larger F value from the F distribution with the residual degrees of freedom, for the "scheffe" method only |
Kyun-Seop Bae [email protected]
x = ModelMatrix(uptake ~ Type, CO2) rx = REG(uptake ~ Type, CO2, summarize=FALSE) est(t(c(0, -1, 1)), x$X, rx) # Quebec - Mississippi t.test(uptake ~ Type, CO2) # compare with the abovex = ModelMatrix(uptake ~ Type, CO2) rx = REG(uptake ~ Type, CO2, summarize=FALSE) est(t(c(0, -1, 1)), x$X, rx) # Quebec - Mississippi t.test(uptake ~ Type, CO2) # compare with the above
Estimates Linear Function with a formula and a dataset.
ESTM(L, Formula, Data, conf.level=0.95)ESTM(L, Formula, Data, conf.level=0.95)
L |
a matrix of linear function rows to be tested |
Formula |
a conventional formula for a linear model |
Data |
a |
conf.level |
confidence level of the confidence limit |
It tests rows of linear functions. A linear function means a linear combination of estimated coefficients. It is similar to the ESTIMATE statement of SAS PROC GLM. This is a convenient version of the est function.
Estimate |
point estimate of the input linear contrast |
Lower CL |
lower confidence limit |
Upper CL |
upper confidence limit |
Std. Error |
standard error of the point estimate |
t value |
value for the t distribution |
Df |
degrees of freedom |
Pr(>|t|) |
probability of a larger absolute t value from the t distribution with the residual degrees of freedom |
Kyun-Seop Bae [email protected]
ESTM(t(c(0, -1, 1)), uptake ~ Type, CO2) # Quevec - MississippiESTM(t(c(0, -1, 1)), uptake ~ Type, CO2) # Quevec - Mississippi
Checks the estimability of row vectors of coefficients.
estmb(L, X, g2, eps=1e-8)estmb(L, X, g2, eps=1e-8)
L |
row vectors of coefficients |
X |
a model (design) matrix from |
g2 |
g2 generalized inverse of |
eps |
An absolute value less than this is considered to be zero. |
It checks the estimability of L, row vectors of coefficients. This corresponds to the ESTIMATE statement of SAS PROC GLM. See <Kennedy Jr. WJ, Gentle JE. Statistical Computing. 1980> p361 or <Golub GH, Styan GP. Numerical Computations for Univariate Linear Models. 1971>.
a vector of logical values indicating which rows are estimable (as TRUE)
Kyun-Seop Bae [email protected]
Exit probabilities with the given drift, upper bounds, and times of test.
ExitP(Theta, bi, ti=NULL)ExitP(Theta, bi, ti=NULL)
Theta |
drift value defined by Lan-DeMets. See the reference. |
bi |
upper bound z-values |
ti |
times of test. These should be in the range of [0, 1]. If omitted, even intervals are assumed. |
It calculates exit probabilities and cumulative exit probabilities with the given drift, upper z-bounds, and times of test. If the times of test are not given, even intervals are assumed. mvtnorm::pmvt (with noncentrality) is better than pmvnorm in calculating power and sample size. However, Lan-DeMets used the multivariate normal rather than the multivariate noncentral t distribution. This function follows Lan-DeMets for consistency with previous results.
The result is a matrix.
ti |
time of test |
bi |
upper z-bound |
cum.alpha |
cumulative alpha-value |
Kyun-Seop Bae [email protected]
Reboussin DM, DeMets DL, Kim K, Lan KKG. Computations for group sequential boundaries using the Lan-DeMets function method. Controlled Clinical Trials. 2000;21:190-207.
b0 = seqBound(ti=(1:5)/5)[, "up.bound"] ExitP(Theta = Drift(b0), bi = b0)b0 = seqBound(ti=(1:5)/5)[, "up.bound"] ExitP(Theta = Drift(b0), bi = b0)
A generalized inverse is usually not unique. Some programs use this algorithm to get a unique generalized inverse matrix. This uses the SWEEP operator and works for non-square matrices also.
g2inv(A, eps=1e-08)g2inv(A, eps=1e-08)
A |
a matrix to be inverted |
eps |
A value less than this is considered zero. |
See 'SAS Technical Report R106, The Sweep Operator: Its Importance in Statistical Computing' by J. H. Goodnight for details.
g2 inverse
Kyun-Seop Bae [email protected]
Searle SR, Khuri AI. Matrix Algebra Useful for Statistics. 2e. John Wiley and Sons Inc. 2017.
A = matrix(c(1, 2, 4, 3, 3, -1, 2, -2, 5, -4, 0, -7), byrow=TRUE, ncol=4) ; A g2inv(A)A = matrix(c(1, 2, 4, 3, 3, -1, 2, -2, 5, -4, 0, -7), byrow=TRUE, ncol=4) ; A g2inv(A)
A generalized inverse is usually not unique. Some programs use this algorithm to get a unique generalized inverse matrix.
G2SWEEP(A, Augmented=FALSE, eps=1e-08)G2SWEEP(A, Augmented=FALSE, eps=1e-08)
A |
a matrix to be inverted. If |
Augmented |
If this is |
eps |
A value less than this is considered zero. |
The generalized inverse of g2-type is used by some software to do linear regression. See 'SAS Technical Report R106, The Sweep Operator: Its Importance in Statistical Computing' by J. H. Goodnight for details.
when Augmented=FALSE |
ordinary g2 inverse |
when Augmented=TRUE |
g2 inverse and beta hats in the last column and the last row, and the sum of squares error (SSE) in the last cell |
attribute "rank" |
the rank of the input matrix |
Kyun-Seop Bae [email protected]
f1 = uptake ~ Type + Treatment # formula x = ModelMatrix(f1, CO2) # Model matrix and relevant information y = model.frame(f1, CO2)[, 1] # observation vector nc = ncol(x$X) # number of columns of model matrix XpY = crossprod(x$X, y) aXpX = rbind(cbind(crossprod(x$X), XpY), cbind(t(XpY), crossprod(y))) ag2 = G2SWEEP(aXpX, Augmented=TRUE) b = ag2[1:nc, (nc + 1)] ; b # Beta hat iXpX = ag2[1:nc, 1:nc] ; iXpX # g2 inverse of X'X SSE = ag2[(nc + 1), (nc + 1)] ; SSE # Sum of Square Error DFr = nrow(x$X) - attr(ag2, "rank") ; DFr # Degree of freedom for the residual # Compare the below with the above REG(f1, CO2) aov1(f1, CO2)f1 = uptake ~ Type + Treatment # formula x = ModelMatrix(f1, CO2) # Model matrix and relevant information y = model.frame(f1, CO2)[, 1] # observation vector nc = ncol(x$X) # number of columns of model matrix XpY = crossprod(x$X, y) aXpX = rbind(cbind(crossprod(x$X), XpY), cbind(t(XpY), crossprod(y))) ag2 = G2SWEEP(aXpX, Augmented=TRUE) b = ag2[1:nc, (nc + 1)] ; b # Beta hat iXpX = ag2[1:nc, 1:nc] ; iXpX # g2 inverse of X'X SSE = ag2[(nc + 1), (nc + 1)] ; SSE # Sum of Square Error DFr = nrow(x$X) - attr(ag2, "rank") ; DFr # Degree of freedom for the residual # Compare the below with the above REG(f1, CO2) aov1(f1, CO2)
Geometric coefficient of variation in percentage.
geoCV(y)geoCV(y)
y |
a numeric vector |
It removes NA. This is sqrt(exp(var(log(y))) - 1)*100.
Geometric coefficient of variation in percentage.
Kyun-Seop Bae [email protected]
geoCV(mtcars$mpg)geoCV(mtcars$mpg)
Geometric mean without NA values.
geoMean(y)geoMean(y)
y |
a vector of numerics |
It removes NA in the input vector.
geometric mean value
Kyun-Seop Bae [email protected]
geoMean(mtcars$mpg)geoMean(mtcars$mpg)
GLM is the main function of this package.
GLM(Formula, Data, BETA=FALSE, EMEAN=FALSE, Resid=FALSE, conf.level=0.95, Weights=1)GLM(Formula, Data, BETA=FALSE, EMEAN=FALSE, Resid=FALSE, conf.level=0.95, Weights=1)
Formula |
a conventional formula for a linear model. |
Data |
a |
BETA |
if |
EMEAN |
if |
Resid |
if |
conf.level |
confidence level for the confidence limit of the least square mean |
Weights |
weights for the weighted least squares. This should be a scalar or a vector of the same length as the number of rows of |
It performs the core function of SAS PROC GLM. Least square means for the interaction term of three variables are not supported yet.
The result is comparable to that of SAS PROC GLM.
ANOVA |
ANOVA table for the model |
Fitness |
Some measures of goodness of fit such as R-square and CV |
Type I |
Type I sum of squares table |
Type II |
Type II sum of squares table |
Type III |
Type III sum of squares table |
Parameter |
Parameter table with standard error, t value, p value. |
Expected Mean |
Least square (or expected) mean table with confidence limits. This is returned only with the EMEAN=TRUE option. |
Fitted |
Fitted values or y hat in the original scale, as SAS OUTPUT P= does, even with Weights. This is returned only with the Resid=TRUE option. |
Residual |
Residuals in the original scale, as SAS OUTPUT R= does, even with Weights. This is returned only with the Resid=TRUE option. |
Kyun-Seop Bae [email protected]
REG, aov1, aov2, aov3, LSM, PDIFF
GLM(uptake ~ Type*Treatment + conc, CO2[-1,]) # Making data unbalanced GLM(uptake ~ Type*Treatment + conc, CO2[-1,], BETA=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], EMEAN=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], Resid=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], BETA=TRUE, EMEAN=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], BETA=TRUE, EMEAN=TRUE, Resid=TRUE)GLM(uptake ~ Type*Treatment + conc, CO2[-1,]) # Making data unbalanced GLM(uptake ~ Type*Treatment + conc, CO2[-1,], BETA=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], EMEAN=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], Resid=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], BETA=TRUE, EMEAN=TRUE) GLM(uptake ~ Type*Treatment + conc, CO2[-1,], BETA=TRUE, EMEAN=TRUE, Resid=TRUE)
Generate a self-contained PDF report documenting whether the sasLM package is correctly installed, intact, loadable, and operational in the user's own R environment. It is intended as Installation Qualification (IQ) evidence, in the spirit of commercial statistical-software validation. The report uses only base R and the package's own pdf helpers, so it requires no LaTeX, pandoc, or other external tools. The signature (approval) page is placed first and the report uses 1 inch margins on every side, with Letter paper in a United States locale and A4 elsewhere.
IQLM(fileName = "sasLM-IQ-Report.pdf", pkgs = "sasLM", functional = TRUE, performedBy = "", paper = "auto", sigField = FALSE)IQLM(fileName = "sasLM-IQ-Report.pdf", pkgs = "sasLM", functional = TRUE, performedBy = "", paper = "auto", sigField = FALSE)
fileName |
file name to save the PDF report. |
pkgs |
character vector of package names to qualify; defaults to
|
functional |
if |
performedBy |
name of the person performing the qualification, printed on the signature page. Defaults to the login name. |
paper |
paper size: |
sigField |
if |
The report contains: the test environment (R version, platform, OS, locale, library
paths); the installed package version, location and declared-dependency satisfaction;
file integrity via checkMD5sums (PASS / FAIL / WARN when no
manifest is present - see writeMD5LM); namespace load and core exports;
an optional functional verification; an overall QUALIFIED / NOT QUALIFIED verdict; a
sessionInfo appendix; and a per-file md5 checksum appendix.
Invisibly, a list with fileName, qualified, paper, checks
(a data frame of every check), and the counts nPass, nFail,
nWarn. A PDF file is written to fileName.
Kyun-Seop Bae <[email protected]>
#IQLM() #res <- IQLM(performedBy = "Jane Doe"); res$qualified#IQLM() #res <- IQLM(performedBy = "Jane Doe"); res$qualified
Tests if the input matrix is a correlation matrix or not.
is.cor(m, eps=1e-16)is.cor(m, eps=1e-16)
m |
a presumed correlation matrix |
eps |
epsilon value. An absolute value less than this is considered zero. |
A diagonal component does not need to be exactly 1, but it should be close to 1.
TRUE or FALSE
Kyun-Seop Bae [email protected]
Kurtosis with a conventional formula.
Kurtosis(y)Kurtosis(y)
y |
a vector of numerics |
It removes NA in the input vector.
Estimate of kurtosis
Kyun-Seop Bae [email protected]
Standard error of the estimated kurtosis with a conventional formula.
KurtosisSE(y)KurtosisSE(y)
y |
a vector of numerics |
It removes NA in the input vector.
Standard error of the estimated kurtosis
Kyun-Seop Bae [email protected]
The estimate of the lower bound of the confidence limit using the t-distribution
LCL(y, conf.level=0.95)LCL(y, conf.level=0.95)
y |
a vector of numerics |
conf.level |
confidence level |
It removes NA in the input vector.
The estimate of the lower bound of the confidence limit using the t-distribution
Kyun-Seop Bae [email protected]
Fits a least squares linear model.
lfit(x, y, eps=1e-8)lfit(x, y, eps=1e-8)
x |
a result of |
y |
a column vector of the response (dependent) variable |
eps |
A value less than this is considered zero. |
A minimal version of the least squares fit of a linear model
coefficients |
beta coefficients |
g2 |
g2 inverse |
rank |
rank of the model matrix |
DFr |
degrees of freedom for the residual |
SSE |
sum of squares error |
SST |
sum of squares total |
DFr2 |
degrees of freedom of the residual for the beta coefficient |
Kyun-Seop Bae [email protected]
f1 = uptake ~ Type*Treatment + conc x = ModelMatrix(f1, CO2) y = model.frame(f1, CO2)[,1] lfit(x, y)f1 = uptake ~ Type*Treatment + conc x = ModelMatrix(f1, CO2) y = model.frame(f1, CO2)[,1] lfit(x, y)
Coefficients are calculated with the g2 inverse. The output is similar to summary(lm()).
lr(Formula, Data, eps=1e-8)lr(Formula, Data, eps=1e-8)
Formula |
a conventional formula for a linear model |
Data |
a |
eps |
A value less than this is considered zero. |
It uses G2SWEEP to get the g2 inverse. The result is similar to summary(lm()) without options.
The result is comparable to that of SAS PROC REG.
Estimate |
point estimate of parameters, coefficients |
Std. Error |
standard error of the point estimate |
t value |
value for the t distribution |
Pr(>|t|) |
probability of a larger absolute t value from the t distribution with the residual degrees of freedom |
Kyun-Seop Bae [email protected]
lr(uptake ~ Plant + Type + Treatment + conc, CO2) lr(uptake ~ Plant + Type + Treatment + conc - 1, CO2) lr(uptake ~ Type, CO2) lr(uptake ~ Type - 1, CO2)lr(uptake ~ Plant + Type + Treatment + conc, CO2) lr(uptake ~ Plant + Type + Treatment + conc - 1, CO2) lr(uptake ~ Type, CO2) lr(uptake ~ Type - 1, CO2)
Usually, the first step in multiple linear regression is to perform simple linear regressions with each single independent variable.
lr0(Formula, Data)lr0(Formula, Data)
Formula |
a conventional formula for a linear model. The intercept will always be added. |
Data |
a |
It performs simple linear regression for each independent variable.
Each row means one simple linear regression with that row name as the only independent variable.
Intercept |
estimate of the intercept |
SE(Intercept) |
standard error of the intercept |
Slope |
estimate of the slope |
SE(Slope) |
standard error of the slope |
Rsq |
R-squared for the simple linear model |
Pr(>F) |
p-value of the slope or the model |
Kyun-Seop Bae [email protected]
lr0(uptake ~ Plant + Type + Treatment + conc, CO2) lr0(mpg ~ ., mtcars)lr0(uptake ~ Plant + Type + Treatment + conc, CO2) lr0(mpg ~ ., mtcars)
Estimates least squares means using the g2 inverse.
LSM(Formula, Data, Term, conf.level=0.95, adj="lsd", hideNonEst=TRUE, PLOT=FALSE, descend=FALSE, ...)LSM(Formula, Data, Term, conf.level=0.95, adj="lsd", hideNonEst=TRUE, PLOT=FALSE, descend=FALSE, ...)
Formula |
a conventional formula for a linear model |
Data |
a |
Term |
a term name to be returned. If there is only one independent variable, this can be omitted. |
conf.level |
confidence level for the confidence limit |
adj |
adjustment method for grouping; "lsd" (default), "tukey", "bon", "duncan", and "scheffe" are available. This does not affect the SE, Lower CL, and Upper CL of the output table. |
hideNonEst |
logical. whether to hide non-estimable values |
PLOT |
logical. whether to plot LSMs and their confidence intervals |
descend |
logical. This specifies whether the plotting order is ascending or descending. |
... |
arguments to be passed to |
It is equivalent to the LSMEANS statement of SAS PROC GLM. The result of the second example below may be different from emmeans. This is because SAS and this function calculate the mean of the transformed continuous variable. However, emmeans calculates the average before the transformation. An interaction of three variables is not supported yet. For the "dunnett" adjustment method, see the PDIFF function.
Returns a table of expectations, t values and p-values.
Group |
group character. This appears when the model is one-way ANOVA or when the |
LSmean |
point estimate of the least squares mean |
LowerCL |
lower confidence limit at the given confidence level by the "lsd" method |
UpperCL |
upper confidence limit at the given confidence level by the "lsd" method |
SE |
standard error of the point estimate |
Df |
degrees of freedom of the point estimate |
Kyun-Seop Bae [email protected]
LSM(uptake ~ Type, CO2[-1,]) LSM(uptake ~ Type - 1, CO2[-1,]) LSM(uptake ~ Type*Treatment + conc, CO2[-1,]) LSM(uptake ~ Type*Treatment + conc - 1, CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + log(conc), CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + log(conc) - 1, CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + as.factor(conc), CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + as.factor(conc) - 1, CO2[-1,]) LSM(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata) LSM(log(CMAX) ~ SEQ/SUBJ + PRD + TRT - 1, BEdata)LSM(uptake ~ Type, CO2[-1,]) LSM(uptake ~ Type - 1, CO2[-1,]) LSM(uptake ~ Type*Treatment + conc, CO2[-1,]) LSM(uptake ~ Type*Treatment + conc - 1, CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + log(conc), CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + log(conc) - 1, CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + as.factor(conc), CO2[-1,]) LSM(log(uptake) ~ Type*Treatment + as.factor(conc) - 1, CO2[-1,]) LSM(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata) LSM(log(CMAX) ~ SEQ/SUBJ + PRD + TRT - 1, BEdata)
Maximum without NA values.
Max(y)Max(y)
y |
a vector of numerics |
It removes NA values from the input vector.
maximum value
Kyun-Seop Bae [email protected]
Mean without NA values.
Mean(y)Mean(y)
y |
a vector of numerics |
It removes NA values from the input vector.
mean value
Kyun-Seop Bae [email protected]
Median without NA values.
Median(y)Median(y)
y |
a vector of numerics |
It removes NA values from the input vector.
median value
Kyun-Seop Bae [email protected]
Minimum without NA values.
Min(y)Min(y)
y |
a vector of numerics |
It removes NA values from the input vector.
minimum value
Kyun-Seop Bae [email protected]
This model matrix is similar to model.matrix, but it does not omit unnecessary columns.
ModelMatrix(Formula, Data, KeepOrder=FALSE, XpX=FALSE)ModelMatrix(Formula, Data, KeepOrder=FALSE, XpX=FALSE)
Formula |
a conventional formula for a linear model |
Data |
a |
KeepOrder |
If |
XpX |
If |
It makes the model (design) matrix for GLM.
Model matrix and attributes similar to the output of model.matrix.
X |
design matrix, i.e. model matrix |
XpX |
cross-product of the design matrix, X'X |
terms |
detailed information about terms such as formula and labels |
termsIndices |
term indices |
assign |
assignment of columns for each term in order, a different way of expressing term indices |
Kyun-Seop Bae [email protected]
This is comparable to SAS PROC TTEST, except that it uses summarized input (sufficient statistics).
mtest(m1, s1, n1, m0, s0, n0, conf.level=0.95)mtest(m1, s1, n1, m0, s0, n0, conf.level=0.95)
m1 |
mean of the first (test, active, experimental) group |
s1 |
sample standard deviation of the first group |
n1 |
sample size of the first group |
m0 |
mean of the second (reference, control, placebo) group |
s0 |
sample standard deviation of the second group |
n0 |
sample size of the second group |
conf.level |
confidence level |
This uses summarized input. This also produces confidence intervals of means and variances by group.
The output format is comparable to that of SAS PROC TTEST.
Kyun-Seop Bae [email protected]
mtest(5.4, 10.5, 3529, 5.1, 8.9, 5190) # NEJM 388;15 p1386mtest(5.4, 10.5, 3529, 5.1, 8.9, 5190) # NEJM 388;15 p1386
Number of observations excluding NA values.
N(y)N(y)
y |
a vector of numerics |
It removes NA values from the input vector.
Count of the observations
Kyun-Seop Bae [email protected]
Sequential O'Brien-Fleming upper bounds for the cumulative Z-test on accumulating data. Z values are correlated. This is usually used for group sequential design.
OBFBound(K, alpha=0.05, side=2, ti=NULL, c0=NULL)OBFBound(K, alpha=0.05, side=2, ti=NULL, c0=NULL)
K |
count of tests, including the final one |
alpha |
goal alpha value for the last test at time 0. |
side |
1=one-sided test, 2=two-sided test |
ti |
times for test. These should be in [0, 1]. If not specified, equal intervals are assumed. |
c0 |
correlation matrix. If not specified, Brownian motion is assumed. |
It calculates O'Brien-Fleming upper z-bounds and cumulative alpha-values for the repeated test in group sequential design.
The result is a matrix.
ti |
time of test |
z |
O'Brien-Fleming upper z-bound |
cum.alpha |
cumulative alpha-value |
Kyun-Seop Bae [email protected]
OBFBound(K=2) OBFBound(K=3) OBFBound(K=4) OBFBound(K=5)OBFBound(K=2) OBFBound(K=3) OBFBound(K=4) OBFBound(K=5)
Generate a self-contained PDF report that documents whether the sasLM linear-model engine reproduces pre-specified reference results on the user's own machine, within a pre-specified relative tolerance. It is intended as Operational Qualification (OQ) evidence. For each scenario the report lists every computed value - the ANOVA table, the fit statistics, and the Type I, II and III sums of squares - next to its reference value, with the absolute and relative difference and a pass flag, then an overall verdict. It uses only base R and the package's pdf helpers (no LaTeX). The signature page is first; 1 inch margins; Letter paper in a United States locale and A4 elsewhere.
OQLM(fileName = "sasLM-OQ-Report.pdf", cfg = NULL, tol = 0.001, refDir = system.file("OQ", package = "sasLM"), performedBy = "", paper = "auto", sigField = FALSE)OQLM(fileName = "sasLM-OQ-Report.pdf", cfg = NULL, tol = 0.001, refDir = system.file("OQ", package = "sasLM"), performedBy = "", paper = "auto", sigField = FALSE)
fileName |
file name to save the PDF report. |
cfg |
a data frame of scenario definitions. If |
tol |
acceptance tolerance: a value passes when the symmetric relative
difference |
refDir |
directory holding ‘OQcfg.csv’, the input data, and the reference
csv files. Defaults to the |
performedBy |
name of the person performing the qualification, printed on the signature page. Defaults to the login name. |
paper |
paper size: |
sigField |
if |
For each scenario the function reads the (frozen) input data, makes the named columns
factors with af, fits the model with GLM, flattens the
result, and compares every reference value to the corresponding computed value. A
scenario passes only when all of its values pass.
The first release covers GLM with scenarios from Harvey (ARS20-8) and
Goodnight (1976). The bundled reference values are the SAS PROC GLM / textbook
results documented in the package validation report (Validation of 'sasLM'
Package); they establish concordance with SAS, the de-facto reference. A new
GLM scenario can be added with no code change, by adding an input csv, a
reference csv, and one row to ‘OQcfg.csv’ in refDir. Adding a different
procedure (such as REG, aov3, TTEST or
UNIV) additionally requires extending OQLM itself (procedure
dispatch and result flattening).
Invisibly, a list with fileName, qualified, tol,
nScenarios, nFailScenarios, nValues, nFailValues,
paper, and results (a per-scenario list whose comp element is the
full table of every comparison). A PDF is written to fileName.
Kyun-Seop Bae <[email protected]>
#OQLM() #res <- OQLM(); res$qualified; res$results[[1]]$comp#OQLM() #res <- OQLM(); res$qualified; res$results[[1]]$comp
Odds ratio between two groups
OR(y1, n1, y2, n2, conf.level=0.95)OR(y1, n1, y2, n2, conf.level=0.95)
y1 |
positive event count of the test (the first) group |
n1 |
total count of the test (the first) group |
y2 |
positive event count of the control (the second) group |
n2 |
total count of the control (the second) group |
conf.level |
confidence level |
It calculates the odds ratio of two groups. No continuity correction is done here. If you need the percent scale, multiply the output by 100.
The result is a data.frame.
odd1 |
odds from the first group, y1/(n1 - y1) |
odd2 |
odds from the second group, y2/(n2 - y2) |
OR |
odds ratio, odd1/odd2 |
SElog |
standard error of log(OR) |
lower |
lower confidence limit of OR |
upper |
upper confidence limit of OR |
Kyun-Seop Bae [email protected]
RD, RR, RDmn1, RRmn1, ORmn1, RDmn, RRmn, ORmn
OR(104, 11037, 189, 11034) # no continuity correctionOR(104, 11037, 189, 11034) # no continuity correction
Odds ratio and its confidence interval of two groups with stratification by the Cochran-Mantel-Haenszel method
ORcmh(d0, conf.level=0.95)ORcmh(d0, conf.level=0.95)
d0 |
A |
conf.level |
confidence level |
It calculates the odds ratio and its confidence interval of two groups. This can be used for meta-analysis also.
The following output will be returned for each stratum and the common value.
odd1 |
odds from the first group, y1/(n1 - y1) |
odd2 |
odds from the second group, y2/(n2 - y2) |
OR |
odds ratio, odd1/odd2. The point estimate of the common OR is calculated with the MH weights. |
SElog |
standard error of log(OR) |
lower |
lower confidence limit of OR |
upper |
upper confidence limit of OR |
Kyun-Seop Bae [email protected]
RDmn1, RRmn1, ORmn1, RDmn, RRmn, ORmn, RDinv, RRinv, ORinv
d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") ORcmh(d1)d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") ORcmh(d1)
Odds ratio and its confidence interval of two groups with stratification by the inverse variance method
ORinv(d0, conf.level=0.95)ORinv(d0, conf.level=0.95)
d0 |
A |
conf.level |
confidence level |
It calculates the odds ratio and its confidence interval of two groups by the inverse variance method. This supports stratification. This can be used for meta-analysis also.
The following output will be returned for each stratum and the common value.
odd1 |
odds from the first group, y1/(n1 - y1) |
odd2 |
odds from the second group, y2/(n2 - y2) |
OR |
odds ratio, odd1/odd2. The point estimate of the common OR is calculated with the inverse variance weights. |
SElog |
standard error of log(OR) |
lower |
lower confidence limit of OR |
upper |
upper confidence limit of OR |
Kyun-Seop Bae [email protected]
RDmn1, RRmn1, ORmn1, RDmn, RRmn, ORmn, RDinv, RRinv, ORcmh
d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") ORinv(d1)d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") ORinv(d1)
Odds ratio and its score confidence interval of two groups with stratification by the Miettinen and Nurminen method
ORmn(d0, conf.level=0.95, eps=1e-8)ORmn(d0, conf.level=0.95, eps=1e-8)
d0 |
A |
conf.level |
confidence level |
eps |
absolute value less than eps is regarded as negligible |
It calculates the common odds ratio and its score confidence interval of two groups with stratification. The confidence interval is asymmetric, and there is no standard error in the output. For the stratified case, the inverse variance weighted score statistic with the bias correction is used, following Laud. The common odds ratio point estimate is the zero of the weighted score statistic, and the confidence bounds are found with the uniroot function. The result agrees with ratesci::scoreci(contrast="OR", stratified=TRUE, skew=FALSE) to at least 7 significant digits. For a single stratum, it returns the classical Miettinen-Nurminen interval of ORmn1. This can be used for meta-analysis also.
The following output will be returned for each stratum and common value. There is no standard error.
odd1 |
odds from the first group, y1/(n1 - y1). For the common value, it is calculated with the weights at the point estimate. |
odd2 |
odds from the second group, y2/(n2 - y2). For the common value, it is calculated with the weights at the point estimate. |
OR |
odds ratio of the stratum. The common OR is the zero of the inverse variance weighted score statistic. |
lower |
lower confidence limit of OR |
upper |
upper confidence limit of OR |
Kyun-Seop Bae [email protected]
Miettinen O, Nurminen M. Comparative analysis of two rates. Stat Med 1985;4:213-26
Laud PJ. Equal-tailed confidence intervals for comparison of rates. Pharmaceutical Statistics 2017;16:334-348
RDmn1, RRmn1, ORmn1, RDmn, RRmn, RDinv, RRinv, ORinv, ORcmh
d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") ORmn(d1) d2 = data.frame(y1=c(4, 2, 10), n1=c(20, 20, 20), y2=c(8, 11, 2), n2=c(20, 20, 20)) ORmn(d2)d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") ORmn(d1) d2 = data.frame(y1=c(4, 2, 10), n1=c(20, 20, 20), y2=c(8, 11, 2), n2=c(20, 20, 20)) ORmn(d2)
Odds ratio and its score confidence interval of two groups without stratification
ORmn1(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)ORmn1(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)
y1 |
positive event count of the test (the first) group |
n1 |
total count of the test (the first) group |
y2 |
positive event count of the control (the second) group |
n2 |
total count of the control (the second) group |
conf.level |
confidence level |
eps |
absolute value less than eps is regarded as negligible |
It calculates the odds ratio and its score confidence interval of the two groups. The confidence interval is asymmetric, and there is no standard error in the output. This does not support stratification. This implementation uses the uniroot function, which usually gives at least 5 significant digits. In contrast, the PropCIs::orscoreci function uses an incremental or decremental search by a factor of 1.001, which gives less than 3 significant digits.
There is no standard error.
odd1 |
odds from the first group, y1/(n1 - y1) |
odd2 |
odds from the second group, y2/(n2 - y2) |
OR |
odds ratio, odd1/odd2 |
lower |
lower confidence limit of OR |
upper |
upper confidence limit of OR |
Kyun-Seop Bae [email protected]
Miettinen O, Nurminen M. Comparative analysis of two rates. Stat Med 1985;4:213-26
RDmn1, RRmn1, RDmn, RRmn, ORmn
ORmn1(104, 11037, 189, 11034)ORmn1(104, 11037, 189, 11034)
It plots bands of the confidence interval and prediction interval for simple linear regression.
pB(Formula, Data, Resol=300, conf.level=0.95, lx, ly, ...)pB(Formula, Data, Resol=300, conf.level=0.95, lx, ly, ...)
Formula |
a formula |
Data |
a |
Resol |
resolution for the output |
conf.level |
confidence level |
lx |
x position of the legend |
ly |
y position of the legend |
... |
arguments to be passed to |
It plots. Discard the return values. If lx or ly is missing, the legend position is calculated automatically.
Ignore the return values.
Kyun-Seop Bae [email protected]
pB(hp ~ disp, mtcars) pB(mpg ~ disp, mtcars)pB(hp ~ disp, mtcars) pB(mpg ~ disp, mtcars)
Testing partial correlation between many columns of data with the Pearson method.
Pcor.test(Data, x, y)Pcor.test(Data, x, y)
Data |
a numeric matrix or |
x |
names of columns to be tested |
y |
names of control columns |
It performs multiple partial correlation tests. It uses "complete.obs" rows of the x and y columns.
Row names show which columns are used for the test.
Estimate |
point estimate of correlation |
Df |
degrees of freedom |
t value |
t value of the t distribution |
Pr(>|t|) |
probability with the t distribution |
Kyun-Seop Bae [email protected]
Pcor.test(mtcars, c("mpg", "hp", "qsec"), c("drat", "wt"))Pcor.test(mtcars, c("mpg", "hp", "qsec"), c("drat", "wt"))
Four standard diagnostic plots for regression.
pD(rx, Title=NULL)pD(rx, Title=NULL)
rx |
a result of |
Title |
title to be printed on the plot |
The most frequently used diagnostic plots are 'observed vs. fitted', 'standardized residual vs. fitted', 'distribution plot of standardized residuals', and 'Q-Q plot of standardized residuals'.
Four diagnostic plots on a page.
Kyun-Seop Bae [email protected]
pD(lm(uptake ~ Plant + Type + Treatment + conc, CO2), "Diagnostic Plot")pD(lm(uptake ~ Plant + Type + Treatment + conc, CO2), "Diagnostic Plot")
Estimates pairwise differences by a common method.
PDIFF(Formula, Data, Term, conf.level=0.95, adj="lsd", ref, PLOT=FALSE, reverse=FALSE, ...)PDIFF(Formula, Data, Term, conf.level=0.95, adj="lsd", ref, PLOT=FALSE, reverse=FALSE, ...)
Formula |
a conventional formula for a linear model |
Data |
a |
Term |
a factor name to be estimated |
conf.level |
confidence level of the confidence interval |
adj |
"lsd", "tukey", "scheffe", "bon", "duncan", or "dunnett" to adjust the p-value and confidence limit |
ref |
reference or control level for the Dunnett test. If missing, the first level is used, as SAS does. |
PLOT |
whether to plot the diffogram |
reverse |
reverse A - B to B - A |
... |
arguments to be passed to |
It corresponds to the PDIFF option of SAS PROC GLM.
Returns a table of expectations, t values and p-values. Output columns may vary according to the adjustment option.
Estimate |
point estimate of the input linear contrast |
Lower CL |
lower confidence limit |
Upper CL |
upper confidence limit |
Std. Error |
standard error of the point estimate |
t value |
value for the t distribution |
Df |
degrees of freedom |
Pr(>|t|) |
probability of a larger absolute t value from the t distribution with the residual's degrees of freedom |
Kyun-Seop Bae [email protected]
PDIFF(uptake ~ Type*Treatment + as.factor(conc), CO2, "as.factor(conc)") PDIFF(uptake ~ Type*Treatment + as.factor(conc), CO2, "as.factor(conc)", adj="tukey")PDIFF(uptake ~ Type*Treatment + as.factor(conc), CO2, "as.factor(conc)") PDIFF(uptake ~ Type*Treatment + as.factor(conc), CO2, "as.factor(conc)", adj="tukey")
Cumulative alpha values for the cumulative hypothesis test with a fixed upper bound z-value in group sequential design.
PocockBound(K=2, alpha=0.05, side=2)PocockBound(K=2, alpha=0.05, side=2)
K |
total number of tests |
alpha |
alpha value at the final test |
side |
1=one-sided test, 2=two-sided test |
Pocock suggested a fixed upper bound z-value for the cumulative hypothesis test in group sequential designs.
a fixed upper bound z-value for the K times repeated hypothesis test with a final alpha-value. Attributes are:
ti |
time of test. Equal intervals are assumed. |
cum.alpha |
cumulative alpha value |
Kyun-Seop Bae [email protected]
Reboussin DM, DeMets DL, Kim K, Lan KKG. Computations for group sequential boundaries using the Lan-DeMets function method. Controlled Clinical Trials. 2000;21:190-207.
PocockBound(K=2) # Z-value of upper bound for the two-stage design PocockBound(K=3) # Z-value of upper bound for the two-stage design PocockBound(K=4) # Z-value of upper bound for the two-stage design PocockBound(K=5) # Z-value of upper bound for the two-stage designPocockBound(K=2) # Z-value of upper bound for the two-stage design PocockBound(K=3) # Z-value of upper bound for the two-stage design PocockBound(K=4) # Z-value of upper bound for the two-stage design PocockBound(K=5) # Z-value of upper bound for the two-stage design
Nine residual diagnostic plots.
pResD(rx, Title=NULL)pResD(rx, Title=NULL)
rx |
a result of |
Title |
title to be printed on the plot |
SAS-style residual diagnostic plots.
Nine residual diagnostic plots on a page.
Kyun-Seop Bae [email protected]
pResD(lm(uptake ~ Plant + Type + Treatment + conc, CO2), "Residual Diagnostic Plot")pResD(lm(uptake ~ Plant + Type + Treatment + conc, CO2), "Residual Diagnostic Plot")
Interquartile range (Q3 - Q1) with a conventional formula.
QuartileRange(y, Type=2)QuartileRange(y, Type=2)
y |
a vector of numerics |
Type |
a type specifier to be passed to the |
It removes NA in the input vector. Type 2 is the SAS default, while Type 6 is the SPSS default.
The value of the interquartile range
Kyun-Seop Bae [email protected]
The range, maximum - minimum, as a scalar value.
Range(y)Range(y)
y |
a vector of numerics |
It removes NA in the input vector.
A scalar value of the range
Kyun-Seop Bae [email protected]
Hypothesis test with a specified type of SS using random effects as error terms. This corresponds to SAS PROC GLM's RANDOM /TEST statement.
RanTest(Formula, Data, Random="", Type=3, eps=1e-8)RanTest(Formula, Data, Random="", Type=3, eps=1e-8)
Formula |
a conventional formula for a linear model |
Data |
a |
Random |
a vector of random effects. All should be specified as primary terms, not as interaction terms. All interaction terms with a random factor are regarded as random effects. |
Type |
Sum of squares type to be used as contrast |
eps |
A value less than this is considered as zero. |
Type can be from 1 to 3. All interaction terms with a random factor are regarded as random effects. Here the error term should not be the MSE.
Returns ANOVA and E(MS) tables with the specified type of SS.
Kyun-Seop Bae [email protected]
RanTest(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, Random="SUBJ") fBE = log(CMAX) ~ ADM/SEQ/SUBJ + PRD + TRT RanTest(fBE, BEdata, Random=c("ADM", "SUBJ")) RanTest(fBE, BEdata, Random=c("ADM", "SUBJ"), Type=2) RanTest(fBE, BEdata, Random=c("ADM", "SUBJ"), Type=1)RanTest(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, Random="SUBJ") fBE = log(CMAX) ~ ADM/SEQ/SUBJ + PRD + TRT RanTest(fBE, BEdata, Random=c("ADM", "SUBJ")) RanTest(fBE, BEdata, Random=c("ADM", "SUBJ"), Type=2) RanTest(fBE, BEdata, Random=c("ADM", "SUBJ"), Type=1)
Risk (proportion) difference between two groups
RD(y1, n1, y2, n2, conf.level=0.95)RD(y1, n1, y2, n2, conf.level=0.95)
y1 |
positive event count of the test (the first) group |
n1 |
total count of the test (the first) group |
y2 |
positive event count of the control (the second) group |
n2 |
total count of the control (the second) group |
conf.level |
confidence level |
It calculates the risk difference between the two groups. No continuity correction here. If you need the percent scale, multiply the output by 100.
The result is a data.frame.
p1 |
proportion from the first group |
p2 |
proportion from the second group |
RD |
risk difference, p1 - p2 |
SE |
standard error of RD |
lower |
lower confidence limit of RD |
upper |
upper confidence limit of RD |
Kyun-Seop Bae [email protected]
RR, OR, RDmn1, RRmn1, ORmn1, RDmn, RRmn, ORmn
RD(104, 11037, 189, 11034) # no continuity correctionRD(104, 11037, 189, 11034) # no continuity correction
Risk difference and its confidence interval between two groups with stratification by the inverse variance method
RDinv(d0, conf.level=0.95)RDinv(d0, conf.level=0.95)
d0 |
A |
conf.level |
confidence level |
It calculates the risk difference and its confidence interval between two groups by the inverse variance method. The common risk difference is given by both the fixed effect model and the DerSimonian-Laird random effects model, with Cochran's Q test for heterogeneity. If you need the percent scale, multiply the output by 100. This can be used for meta-analysis also.
RDs |
risk difference and its confidence interval of each stratum |
Heterogeneity |
Cochran's Q statistic for heterogeneity across the strata and its p-value |
tau2 |
between-strata variance estimated by the method of moments |
Fixed |
common risk difference, its standard error, and its confidence interval by the fixed effect model with the inverse variance weights |
Random |
common risk difference, its standard error, and its confidence interval by the DerSimonian-Laird random effects model |
Kyun-Seop Bae [email protected]
RDmn1, RRmn1, ORmn1, RDmn, RRmn, ORmn, RRinv, ORinv, ORcmh
d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RDinv(d1)d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RDinv(d1)
Risk difference and its score confidence interval between two groups with stratification by the Miettinen and Nurminen method
RDmn(d0, conf.level=0.95, eps=1e-8)RDmn(d0, conf.level=0.95, eps=1e-8)
d0 |
A |
conf.level |
confidence level |
eps |
an absolute value less than eps is regarded as negligible |
It calculates the risk difference and its score confidence interval between the two groups. The confidence interval is asymmetric, and there is no standard error in the output. If you need the percent scale, multiply the output by 100. This supports stratification. This implementation uses the uniroot function, which usually gives at least 5 significant digits. This can be used for meta-analysis also.
The following output will be returned for each stratum and the common value. There is no standard error.
p1 |
proportion from the first group, y1/n1 |
p2 |
proportion from the second group, y2/n2 |
RD |
risk difference, p1 - p2. The point estimate of the common RD is calculated with the MN weight. |
lower |
lower confidence limit of RD |
upper |
upper confidence limit of RD |
Kyun-Seop Bae [email protected]
Miettinen O, Nurminen M. Comparative analysis of two rates. Stat Med 1985;4:213-26
RDmn1, RRmn1, ORmn1, RRmn, ORmn, RDinv, RRinv, ORinv, ORcmh
d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RDmn(d1) d2 = data.frame(y1=c(4, 2, 10), n1=c(20, 20, 20), y2=c(8, 11, 2), n2=c(20, 20, 20)) RDmn(d2)d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RDmn(d1) d2 = data.frame(y1=c(4, 2, 10), n1=c(20, 20, 20), y2=c(8, 11, 2), n2=c(20, 20, 20)) RDmn(d2)
Risk difference and its score confidence interval between two groups without stratification
RDmn1(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)RDmn1(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)
y1 |
positive event count of the test (the first) group |
n1 |
total count of the test (the first) group. The maximum allowable value is 1e8. |
y2 |
positive event count of the control (the second) group |
n2 |
total count of the control (the second) group. The maximum allowable value is 1e8. |
conf.level |
confidence level |
eps |
an absolute value less than eps is regarded as negligible |
It calculates the risk difference and its score confidence interval between the two groups. The confidence interval is asymmetric, and there is no standard error in the output. If you need the percent scale, multiply the output by 100. This does not support stratification. This implementation uses the uniroot function, which usually gives at least 5 significant digits.
There is no standard error.
p1 |
proportion from the first group, y1/n1 |
p2 |
proportion from the second group, y2/n2 |
RD |
risk difference, p1 - p2 |
lower |
lower confidence limit of RD |
upper |
upper confidence limit of RD |
Kyun-Seop Bae [email protected]
Miettinen O, Nurminen M. Comparative analysis of two rates. Stat Med 1985;4:213-26
RRmn1, ORmn1, RDmn, RRmn, ORmn
RDmn1(104, 11037, 189, 11034)RDmn1(104, 11037, 189, 11034)
REG is similar to SAS PROC REG.
REG(Formula, Data, conf.level=0.95, HC=FALSE, Resid=FALSE, Weights=1, summarize=TRUE)REG(Formula, Data, conf.level=0.95, HC=FALSE, Resid=FALSE, Weights=1, summarize=TRUE)
Formula |
a conventional formula for a linear model |
Data |
a |
conf.level |
confidence level for the confidence limit |
HC |
heteroscedasticity-related output is required, such as HC0, HC3, and White's first and second moment specification test |
Resid |
if |
Weights |
weights for each observation, usually the inverse of each variance. This should be a scalar or a vector of the same length as the number of rows of |
summarize |
If this is |
It performs the core function of SAS PROC REG.
The result is comparable to that of SAS PROC REG.
The first part is the ANOVA table.
The second part is measures of fitness.
The third part is the estimates of coefficients.
Estimate |
point estimate of parameters, coefficients |
Estimable |
estimability: 1=TRUE, 0=FALSE. This appears only when at least one inestimability occurs. |
Std. Error |
standard error of the point estimate |
Lower CL |
lower confidence limit with conf.level |
Upper CL |
upper confidence limit with conf.level |
Df |
degrees of freedom |
t value |
value for the t distribution |
Pr(>|t|) |
probability of a larger absolute t value from the t distribution with the residual's degrees of freedom |
The above result is repeated using HC0 and HC3, followed by White's first and second moment specification test, if the HC option is specified. The t values and their p values with HC1 and HC2 are between those of HC0 and HC3.
Fitted |
Fitted value or y hat in the original scale as SAS OUTPUT P= does, even with Weights. This is returned only with the Resid=TRUE option. |
Residual |
Residuals in the original scale as SAS OUTPUT R= does, even with Weights. This is returned only with the Resid=TRUE option. |
If summarize=FALSE, REG returns;
coefficients |
beta coefficients |
g2 |
g2 inverse |
rank |
rank of the model matrix |
DFr |
degrees of freedom for the residual |
SSE |
sum of squared errors |
Kyun-Seop Bae [email protected]
REG(uptake ~ Plant + Type + Treatment + conc, CO2) REG(uptake ~ conc, CO2, HC=TRUE) REG(uptake ~ conc, CO2, Resid=TRUE) REG(uptake ~ conc, CO2, HC=TRUE, Resid=TRUE) REG(uptake ~ conc, CO2, summarize=FALSE)REG(uptake ~ Plant + Type + Treatment + conc, CO2) REG(uptake ~ conc, CO2, HC=TRUE) REG(uptake ~ conc, CO2, Resid=TRUE) REG(uptake ~ conc, CO2, HC=TRUE, Resid=TRUE) REG(uptake ~ conc, CO2, summarize=FALSE)
regD provides rich diagnostics such as student residual, leverage (hat), Cook's D, studentized deleted residual, DFFITS, and DFBETAS.
regD(Formula, Data)regD(Formula, Data)
Formula |
a conventional formula for a linear model |
Data |
a |
It performs the conventional regression analysis. This does not use the g2 inverse; therefore, it cannot handle a singular matrix. If the model (design) matrix is not of full rank, use REG or fewer parameters.
Coefficients |
conventional coefficients summary with Wald statistics |
Diagnostics |
Diagnostics table for detecting outliers or influential/leverage points. This includes the fitted value (Predicted), residual (Residual), standard error of the residual (SE_Resid), studentized residual (Student_Res), hat (Leverage), Cook's D, studentized deleted residual (RStudent), DFFITS, and COVRATIO. |
DFBETAS |
Column names are the names of coefficients. Each row shows how much each coefficient is affected by deleting the corresponding observation. |
Kyun-Seop Bae [email protected]
regD(uptake ~ conc, CO2)regD(uptake ~ conc, CO2)
Relative Risk between the two groups
RR(y1, n1, y2, n2, conf.level=0.95)RR(y1, n1, y2, n2, conf.level=0.95)
y1 |
positive event count of the test (the first) group |
n1 |
total count of the test (the first) group |
y2 |
positive event count of the control (the second) group |
n2 |
total count of the control (the second) group |
conf.level |
confidence level |
It calculates the relative risk of the two groups. No continuity correction here. If you need the percent scale, multiply the output by 100.
The result is a data.frame.
p1 |
proportion from the first group |
p2 |
proportion from the second group |
RR |
relative risk, p1/p2 |
SElog |
standard error of log(RR) |
lower |
lower confidence limit of RR |
upper |
upper confidence limit of RR |
Kyun-Seop Bae [email protected]
RD, OR, RDmn1, RRmn1, ORmn1, RDmn, RRmn, ORmn
RR(104, 11037, 189, 11034) # no continuity correctionRR(104, 11037, 189, 11034) # no continuity correction
Relative risk and its confidence interval of two groups with stratification by the inverse variance method
RRinv(d0, conf.level=0.95)RRinv(d0, conf.level=0.95)
d0 |
A |
conf.level |
confidence level |
It calculates the relative risk and its confidence interval of two groups by the inverse variance method on the log scale. The common relative risk is given by both the fixed effect model and the DerSimonian-Laird random effects model, with Cochran's Q test for heterogeneity. This can be used for meta-analysis also.
RRs |
relative risk, its confidence interval, and the percent weights (pwi for the fixed effect model, pwsi for the random effects model) of each stratum |
Heterogeneity |
Cochran's Q statistic for heterogeneity across the strata and its p-value |
tau2 |
between-strata variance estimated by the method of moments |
Fixed |
common relative risk and its confidence interval by the fixed effect model with the inverse variance weights |
Random |
common relative risk and its confidence interval by the DerSimonian-Laird random effects model |
Kyun-Seop Bae [email protected]
RDmn1, RRmn1, ORmn1, RDmn, RRmn, ORmn, RDinv, ORinv, ORcmh
d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RRinv(d1)d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RRinv(d1)
Relative risk and its score confidence interval of two groups with stratification by the Miettinen and Nurminen method
RRmn(d0, conf.level=0.95, eps=1e-8)RRmn(d0, conf.level=0.95, eps=1e-8)
d0 |
A |
conf.level |
confidence level |
eps |
an absolute value less than eps is regarded as negligible |
It calculates the relative risk and its score confidence interval of the two groups. The confidence interval is asymmetric, and there is no standard error in the output. This supports stratification. This implementation uses the uniroot function, which usually gives at least 5 significant digits, whereas the PropCIs::riskscoreci function uses a cubic equation approximation which gives only about 2 significant digits. This can be used for meta-analysis also.
The following output will be returned for each stratum and the common value. There is no standard error.
p1 |
proportion from the first group, y1/n1 |
p2 |
proportion from the second group, y2/n2 |
RR |
relative risk, p1/p2. The point estimate of the common RR is calculated with the MN weight. |
lower |
lower confidence limit of RR |
upper |
upper confidence limit of RR |
Kyun-Seop Bae [email protected]
Miettinen O, Nurminen M. Comparative analysis of two rates. Stat Med 1985;4:213-26
RDmn1, RRmn1, ORmn1, RDmn, ORmn, RDinv, RRinv, ORinv, ORcmh
d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RRmn(d1) d2 = data.frame(y1=c(4, 2, 10), n1=c(20, 20, 20), y2=c(8, 11, 2), n2=c(20, 20, 20)) RRmn(d2)d1 = matrix(c(25, 339, 28, 335, 23, 370, 40, 364), nrow=2, byrow=TRUE) colnames(d1) = c("y1", "n1", "y2", "n2") RRmn(d1) d2 = data.frame(y1=c(4, 2, 10), n1=c(20, 20, 20), y2=c(8, 11, 2), n2=c(20, 20, 20)) RRmn(d2)
Relative risk and its score confidence interval of the two groups without stratification
RRmn1(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)RRmn1(y1, n1, y2, n2, conf.level=0.95, eps=1e-8)
y1 |
positive event count of the test (the first) group |
n1 |
total count of the test (the first) group |
y2 |
positive event count of the control (the second) group |
n2 |
total count of the control (the second) group |
conf.level |
confidence level |
eps |
an absolute value less than eps is regarded as negligible |
It calculates the relative risk and its score confidence interval of the two groups. The confidence interval is asymmetric, and there is no standard error in the output. This does not support stratification. This implementation uses the uniroot function, which usually gives at least 5 significant digits, whereas the PropCIs::riskscoreci function uses a cubic equation approximation which gives only about 2 significant digits.
There is no standard error.
p1 |
proportion from the first group, y1/n1 |
p2 |
proportion from the second group, y2/n2 |
RR |
relative risk, p1/p2 |
lower |
lower confidence limit of RR |
upper |
upper confidence limit of RR |
Kyun-Seop Bae [email protected]
Miettinen O, Nurminen M. Comparative analysis of two rates. Stat Med 1985;4:213-26
RDmn1, ORmn1, RDmn, RRmn, ORmn
RRmn1(104, 11037, 189, 11034)RRmn1(104, 11037, 189, 11034)
Calculates the pooled variance and degrees of freedom using the Satterthwaite equation.
satt(vars, dfs, ws=c(1, 1))satt(vars, dfs, ws=c(1, 1))
vars |
a vector of variances |
dfs |
a vector of degrees of freedom |
ws |
a vector of weights |
The input can contain more than two variances.
Variance |
approximated variance |
Df |
degrees of freedom |
Kyun-Seop Bae [email protected]
Score confidence interval of a proportion in one group
ScoreCI(y, n, conf.level=0.95)ScoreCI(y, n, conf.level=0.95)
y |
positive event count of a group |
n |
total count of a group |
conf.level |
confidence level |
It calculates the score confidence interval of a proportion in one group. The confidence interval is asymmetric, and there is no standard error in the output. If you need the percent scale, multiply the output by 100.
The result is a data.frame. There is no standard error.
PE |
point estimate of the proportion |
Lower |
lower confidence limit of the proportion |
Upper |
upper confidence limit of the proportion |
Kyun-Seop Bae [email protected]
ScoreCI(104, 11037)ScoreCI(104, 11037)
Standard deviation of a sample.
SD(y)SD(y)
y |
a vector of numerics |
It removes NA in the input vector. The length of the vector should be larger than 1.
Sample standard deviation
Kyun-Seop Bae [email protected]
The estimate of the standard error of the sample mean
SEM(y)SEM(y)
y |
a vector of numerics |
It removes NA in the input vector.
The estimate of the standard error of the sample mean
Kyun-Seop Bae [email protected]
Sequential upper bounds for cumulative Z-test on accumulating data. Z values are correlated. This is usually used for a group sequential design.
seqBound(ti, alpha = 0.05, side = 2, t2 = NULL, asf = 1)seqBound(ti, alpha = 0.05, side = 2, t2 = NULL, asf = 1)
ti |
times for the tests. These should be in [0, 1]. |
alpha |
goal alpha value for the last test at time 1. |
side |
1=one-sided test, 2=two-sided test |
t2 |
fractions of the information amount. These should be in [0, 1]. If not available, ti will be used instead. |
asf |
alpha spending function. 1=O'Brien-Fleming type (approximate, not exact), 2=Pocock type (approximate, not exact), 3=alpha*ti, 4=alpha*ti^1.5, 5=alpha*ti^2 |
It calculates upper z-bounds and cumulative alpha-values for the repeated tests in a group sequential design. The correlation is assumed to be sqrt(t_i/t_j).
Use PocockBound and OBFBound for more exact bounds.
The result is a matrix.
time |
time of test |
up.bound |
upper z-bound |
cum.alpha |
cumulative alpha-value |
Kyun-Seop Bae [email protected]
Reboussin DM, DeMets DL, Kim K, Lan KKG. Computations for group sequential boundaries using the Lan-DeMets function method. Controlled Clinical Trials. 2000;21:190-207.
seqBound(ti=(1:5)/5) seqBound(ti=(1:5)/5, asf=2)seqBound(ti=(1:5)/5) seqBound(ti=(1:5)/5, asf=2)
Confidence interval with given upper bounds, times of tests, the last Z-value, and confidence level.
seqCI(bi, ti, Zval, conf.level=0.95)seqCI(bi, ti, Zval, conf.level=0.95)
bi |
upper bound z-values |
ti |
times for the tests. These should be in [0, 1]. |
Zval |
the last z-value from the observed data. This is not necessarily the planned final Z-value. |
conf.level |
confidence level |
It calculates the confidence interval with given upper bounds, times of tests, the last Z-value, and confidence level. It assumes a two-sided test. mvtnorm::pmvt (with noncentrality) is better than pmvnorm in calculating power, sample size, and confidence interval. However, Lan-DeMets used the multivariate normal rather than the multivariate noncentral t distribution. This function follows Lan-DeMets for consistency with previous results. For the theoretical background, see the reference.
confidence interval of the Z-value for the given confidence level.
Kyun-Seop Bae [email protected]
Reboussin DM, DeMets DL, Kim K, Lan KKG. Computations for group sequential boundaries using the Lan-DeMets function method. Controlled Clinical Trials. 2000;21:190-207.
seqCI(bi = c(2.53, 2.61, 2.57, 2.47, 2.43, 2.38), ti = c(.2292, .3333, .4375, .5833, .7083, .8333), Zval=2.82)seqCI(bi = c(2.53, 2.61, 2.57, 2.47, 2.43, 2.38), ti = c(.2292, .3333, .4375, .5833, .7083, .8333), Zval=2.82)
Apply a cryptographic digital signature to a (report) PDF, as an alternative to the
print-sign-scan workflow, and verify it. signPDFLM signs the PDF bytes with the
signer's private key (RSA or EC) using SHA-256 and writes a detached signature
sidecar ‘<pdf>.sig’ (and, by default, the signer's public key
‘<pdf>.pubkey.pem’). verifyPDFLM confirms, with the signer's public key
or certificate, that the PDF is byte-for-byte intact and was signed by that key
(tamper-evidence and non-repudiation). These functions require the openssl
package.
This is a detached signature (a separate ‘.sig’ file), not a PAdES signature embedded inside the PDF; embedding a visible in-viewer signature requires a dedicated PDF tool (e.g. Adobe Acrobat or ‘pyhanko’). The detached signature is cryptographically equivalent for integrity and non-repudiation.
signPDFLM(pdf, key, password = NULL, signer = "", role = "", sigFile = paste0(pdf, ".sig"), writePubkey = TRUE) verifyPDFLM(pdf, sigFile = paste0(pdf, ".sig"), pubkey = paste0(pdf, ".pubkey.pem"))signPDFLM(pdf, key, password = NULL, signer = "", role = "", sigFile = paste0(pdf, ".sig"), writePubkey = TRUE) verifyPDFLM(pdf, sigFile = paste0(pdf, ".sig"), pubkey = paste0(pdf, ".pubkey.pem"))
pdf |
path to the PDF file to sign or verify. |
key |
the signer's private key: a path to a PEM key file or an openssl
key object. Create one once with, e.g., |
password |
password for an encrypted private key, or |
signer |
name of the signer, recorded in the signature sidecar. Defaults to the login name. |
role |
optional role recorded in the sidecar, e.g. |
sigFile |
path of the signature sidecar to write ( |
writePubkey |
if |
pubkey |
the signer's public key or certificate used to verify: a PEM path (public key or X.509 certificate) or an openssl pubkey/cert object. |
signPDFLM invisibly returns a list with the sidecar path, the PDF sha256,
the signer, and the public-key fingerprint; it writes the ‘.sig’ sidecar
as a side effect. verifyPDFLM invisibly returns TRUE only if both the
signature is valid and the recorded hash matches, and prints a human-readable result.
Kyun-Seop Bae <[email protected]>
#key <- openssl::rsa_keygen() #openssl::write_pem(key, "signer_key.pem") #IQLM("sasLM-IQ-Report.pdf", performedBy = "Kyun-Seop Bae") #signPDFLM("sasLM-IQ-Report.pdf", "signer_key.pem", signer = "Kyun-Seop Bae") #verifyPDFLM("sasLM-IQ-Report.pdf")#key <- openssl::rsa_keygen() #openssl::write_pem(key, "signer_key.pem") #IQLM("sasLM-IQ-Report.pdf", performedBy = "Kyun-Seop Bae") #signPDFLM("sasLM-IQ-Report.pdf", "signer_key.pem", signer = "Kyun-Seop Bae") #verifyPDFLM("sasLM-IQ-Report.pdf")
Skewness with a conventional formula.
Skewness(y)Skewness(y)
y |
a vector of numerics |
It removes NA in the input vector.
Estimate of skewness
Kyun-Seop Bae [email protected]
Standard error of the skewness with a conventional formula.
SkewnessSE(y)SkewnessSE(y)
y |
a vector of numerics |
It removes NA in the input vector.
Standard error of the estimated skewness
Kyun-Seop Bae [email protected]
Performs an F test with a given slice term.
SLICE(Formula, Data, Term, By)SLICE(Formula, Data, Term, By)
Formula |
a conventional formula for a linear model |
Data |
a |
Term |
a factor name (not an interaction) to calculate the sum of squares and do an F test with least square means |
By |
a factor name to be used for slicing |
It performs an F test with a given slice term. It is similar to the SLICE option of SAS PROC GLM.
Returns the sum of squares and its F value and p-value. Row names are the levels of the slice term.
Df |
degrees of freedom |
Sum Sq |
sum of squares for the set of contrasts |
Mean Sq |
mean square |
F value |
F value for the F distribution |
Pr(>F) |
probability of a larger F value |
Kyun-Seop Bae [email protected]
SLICE(uptake ~ Type*Treatment, CO2, "Type", "Treatment") SLICE(uptake ~ Type*Treatment, CO2, "Treatment", "Type")SLICE(uptake ~ Type*Treatment, CO2, "Type", "Treatment") SLICE(uptake ~ Type*Treatment, CO2, "Treatment", "Type")
Sum of squares with ANOVA.
SS(x, rx, L, eps=1e-8)SS(x, rx, L, eps=1e-8)
x |
a result of |
rx |
a result of |
L |
linear hypothesis, a full matrix matching the information in |
eps |
Values less than this are considered as zero. |
It calculates the sum of squares and completes the ANOVA table.
ANOVA table |
a classical ANOVA table without the residual (Error) part. |
Kyun-Seop Bae [email protected]
Calculates a formula table for the expected mean square of Type III SS.
T3MS(Formula, Data, L0, eps=1e-8)T3MS(Formula, Data, L0, eps=1e-8)
Formula |
a conventional formula for a linear model |
Data |
a |
L0 |
a matrix of row linear contrasts; if missing, |
eps |
Values less than this are considered as zero. |
This is necessary for further hypothesis tests of nesting factors.
A coefficient matrix for Type III expected mean square
Kyun-Seop Bae [email protected]
T3MS(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata)T3MS(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata)
Hypothesis test of Type III SS using an error term other than MSE. This corresponds to SAS PROC GLM's RANDOM /TEST clause.
T3test(Formula, Data, H="", E="", eps=1e-8)T3test(Formula, Data, H="", E="", eps=1e-8)
Formula |
a conventional formula for a linear model |
Data |
a |
H |
Hypothesis term |
E |
Error term |
eps |
Values less than this are considered as zero. |
It tests a factor of Type III SS using some other term as an error term. Here the error term should not be MSE.
Returns one or more ANOVA table(s) of Type III SS.
Kyun-Seop Bae [email protected]
T3test(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, E=c("SEQ:SUBJ")) T3test(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, H="SEQ", E=c("SEQ:SUBJ"))T3test(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, E=c("SEQ:SUBJ")) T3test(log(CMAX) ~ SEQ/SUBJ + PRD + TRT, BEdata, H="SEQ", E=c("SEQ:SUBJ"))
This produces essentially the same result as t.test, except using summarized input (sufficient statistics).
tmtest(m1, s1, n1, m0, s0, n0, conf.level=0.95, nullHypo=0, var.equal=FALSE)tmtest(m1, s1, n1, m0, s0, n0, conf.level=0.95, nullHypo=0, var.equal=FALSE)
m1 |
mean of the first (test, active, experimental) group |
s1 |
sample standard deviation of the first group |
n1 |
sample size of the first group |
m0 |
mean of the second (reference, control, placebo) group |
s0 |
sample standard deviation of the second group |
n0 |
sample size of the second group |
conf.level |
confidence level |
nullHypo |
value for the difference of means under the null hypothesis |
var.equal |
assumption on the variance equality |
The default is the Welch t-test with the Satterthwaite approximation.
The output format is very similar to that of t.test.
Kyun-Seop Bae [email protected]
tmtest(5.4, 10.5, 3529, 5.1, 8.9, 5190) # NEJM 388;15 p1386 tmtest(5.4, 10.5, 3529, 5.1, 8.9, 5190, var.equal=TRUE)tmtest(5.4, 10.5, 3529, 5.1, 8.9, 5190) # NEJM 388;15 p1386 tmtest(5.4, 10.5, 3529, 5.1, 8.9, 5190, var.equal=TRUE)
Trimmed mean wrapping the mean function.
trimmedMean(y, Trim=0.05)trimmedMean(y, Trim=0.05)
y |
a vector of numerics |
Trim |
trimming proportion. Default is 0.05 |
It removes NA in the input vector.
The value of the trimmed mean
Kyun-Seop Bae [email protected]
Summarize a continuous dependent variable with or without independent variables.
tsum(Formula=NULL, Data=NULL, ColNames=NULL, MaxLevel=30, ...)tsum(Formula=NULL, Data=NULL, ColNames=NULL, MaxLevel=30, ...)
Formula |
a conventional formula |
Data |
a |
ColNames |
If there is no |
MaxLevel |
An independent variable with more levels than this will not be handled. |
... |
arguments to be passed to |
A convenient summarization function for a continuous variable. This is a wrapper function for tsum0, tsum1, tsum2, or tsum3.
A data.frame of descriptive summarization values.
Kyun-Seop Bae [email protected]
tsum(lh) t(tsum(CO2)) t(tsum(uptake ~ Treatment, CO2)) tsum(uptake ~ Type + Treatment, CO2) print(tsum(uptake ~ conc + Type + Treatment, CO2), digits=3)tsum(lh) t(tsum(CO2)) t(tsum(uptake ~ Treatment, CO2)) tsum(uptake ~ Type + Treatment, CO2) print(tsum(uptake ~ conc + Type + Treatment, CO2), digits=3)
Summarize a continuous dependent (y) variable without any independent (x) variable.
tsum0(d, y, e=c("Mean", "SD", "N"), repl=list(c("length"), c("n")))tsum0(d, y, e=c("Mean", "SD", "N"), repl=list(c("length"), c("n")))
d |
a |
y |
y variable name, a continuous variable |
e |
a vector of summary function names |
repl |
a list of strings to replace after summarization. The length of the list should be 2, and both elements should have the same length. |
A convenient summarization function for a continuous variable.
A vector of summarized values
Kyun-Seop Bae [email protected]
tsum0(CO2, "uptake") tsum0(CO2, "uptake", repl=list(c("mean", "length"), c("Mean", "n")))tsum0(CO2, "uptake") tsum0(CO2, "uptake", repl=list(c("mean", "length"), c("Mean", "n")))
Summarize a continuous dependent (y) variable with one independent (x) variable.
tsum1(d, y, u, e=c("Mean", "SD", "N"), ou="", repl=list(c("length"), ("n")))tsum1(d, y, u, e=c("Mean", "SD", "N"), ou="", repl=list(c("length"), ("n")))
d |
a |
y |
y variable name, a continuous variable |
u |
x variable name, upper side variable |
e |
a vector of summary function names |
ou |
order of levels of the upper side x variable |
repl |
a list of strings to replace after summarization. The length of the list should be 2, and both elements should have the same length. |
A convenient summarization function for a continuous variable with one x variable.
A data.frame of summarized values. Row names are from e names. Column names are from the levels of the x variable.
Kyun-Seop Bae [email protected]
tsum1(CO2, "uptake", "Treatment") tsum1(CO2, "uptake", "Treatment", e=c("mean", "median", "sd", "min", "max", "length"), ou=c("chilled", "nonchilled"), repl=list(c("median", "length"), c("med", "n")))tsum1(CO2, "uptake", "Treatment") tsum1(CO2, "uptake", "Treatment", e=c("mean", "median", "sd", "min", "max", "length"), ou=c("chilled", "nonchilled"), repl=list(c("median", "length"), c("med", "n")))
Summarize a continuous dependent (y) variable with two independent (x) variables.
tsum2(d, y, l, u, e=c("Mean", "SD", "N"), h=NULL, ol="", ou="", rm.dup=TRUE, repl=list(c("length"), c("n")))tsum2(d, y, l, u, e=c("Mean", "SD", "N"), h=NULL, ol="", ou="", rm.dup=TRUE, repl=list(c("length"), c("n")))
d |
a |
y |
y variable name, a continuous variable |
l |
x variable name to be shown on the left side |
u |
x variable name to be shown on the upper side |
e |
a vector of summary function names |
h |
a vector of summary function names for the horizontal subgroup. If |
ol |
order of levels of the left side x variable |
ou |
order of levels of the upper side x variable |
rm.dup |
if |
repl |
a list of strings to replace after summarization. The length of the list should be 2, and both elements should have the same length. |
A convenient summarization function for a continuous variable with two x variables; one on the left side, the other on the upper side.
A data.frame of summarized values. Column names are from the levels of u. Row names are basically from the levels of l.
Kyun-Seop Bae [email protected]
tsum2(CO2, "uptake", "Type", "Treatment") tsum2(CO2, "uptake", "Type", "conc") tsum2(CO2, "uptake", "Type", "Treatment", e=c("mean", "median", "sd", "min", "max", "length"), ou=c("chilled", "nonchilled"), repl=list(c("median", "length"), c("med", "n")))tsum2(CO2, "uptake", "Type", "Treatment") tsum2(CO2, "uptake", "Type", "conc") tsum2(CO2, "uptake", "Type", "Treatment", e=c("mean", "median", "sd", "min", "max", "length"), ou=c("chilled", "nonchilled"), repl=list(c("median", "length"), c("med", "n")))
Summarize a continuous dependent (y) variable with three independent (x) variables.
tsum3(d, y, l, u, e=c("Mean", "SD", "N"), h=NULL, ol1="", ol2="", ou="", rm.dup=TRUE, repl=list(c("length"), c("n")))tsum3(d, y, l, u, e=c("Mean", "SD", "N"), h=NULL, ol1="", ol2="", ou="", rm.dup=TRUE, repl=list(c("length"), c("n")))
d |
a |
y |
y variable name, a continuous variable |
l |
a vector of two x variable names to be shown on the left side. The length should be 2. |
u |
x variable name to be shown on the upper side |
e |
a vector of summary function names |
h |
a list of two vectors of summary function names for the first and second horizontal subgroups. If |
ol1 |
order of levels of the 1st left side x variable |
ol2 |
order of levels of the 2nd left side x variable |
ou |
order of levels of the upper side x variable |
rm.dup |
if |
repl |
a list of strings to replace after summarization. The length of the list should be 2, and both elements should have the same length. |
A convenient summarization function for a continuous variable with three x variables; two on the left side, the other on the upper side.
A data.frame of summarized values. Column names are from the levels of u. Row names are basically from the levels of l.
Kyun-Seop Bae [email protected]
tsum3(CO2, "uptake", c("Type", "Treatment"), "conc") tsum3(CO2, "uptake", c("Type", "Treatment"), "conc", e=c("mean", "median", "sd", "min", "max", "length"), h=list(c("mean", "sd", "length"), c("mean", "length")), ol2=c("chilled", "nonchilled"), repl=list(c("median", "length"), c("med", "n")))tsum3(CO2, "uptake", c("Type", "Treatment"), "conc") tsum3(CO2, "uptake", c("Type", "Treatment"), "conc", e=c("mean", "median", "sd", "min", "max", "length"), h=list(c("mean", "sd", "length"), c("mean", "length")), ol2=c("chilled", "nonchilled"), repl=list(c("median", "length"), c("med", "n")))
This is comparable to SAS PROC TTEST.
TTEST(x, y, conf.level=0.95)TTEST(x, y, conf.level=0.95)
x |
a vector of data from the first (test, active, experimental) group |
y |
a vector of data from the second (reference, control, placebo) group |
conf.level |
confidence level |
Be cautious when choosing the row to use in the output.
The output format is comparable to that of SAS PROC TTEST.
Kyun-Seop Bae [email protected]
TTEST(mtcars[mtcars$am==1, "mpg"], mtcars[mtcars$am==0, "mpg"])TTEST(mtcars[mtcars$am==1, "mpg"], mtcars[mtcars$am==0, "mpg"])
The estimate of the upper bound of the confidence limit using the t-distribution
UCL(y, conf.level=0.95)UCL(y, conf.level=0.95)
y |
a vector of numerics |
conf.level |
confidence level |
It removes NA in the input vector.
The estimate of the upper bound of the confidence limit using the t-distribution
Kyun-Seop Bae [email protected]
Returns descriptive statistics of a numeric vector.
UNIV(y, conf.level = 0.95)UNIV(y, conf.level = 0.95)
y |
a numeric vector |
conf.level |
confidence level for confidence limit |
A convenient and comprehensive function for descriptive statistics. NA is removed during the calculation. This is similar to SAS PROC UNIVARIATE.
nAll |
count of all elements in the input vector |
nNA |
count of |
nFinite |
count of finite numbers |
Mean |
mean excluding |
Variance |
variance excluding |
SD |
standard deviation excluding |
CV |
coefficient of variation in percent |
SEM |
standard error of the sample mean, the sample standard deviation divided by the square root of nFinite |
LowerCL |
lower confidence limit of mean |
UpperCL |
upper confidence limit of mean |
TrimmedMean |
trimmed mean with a trimming proportion of 1 - confidence level |
Min |
minimum value |
Q1 |
first quartile value with quantile type 2, the SAS default |
Median |
median value |
Q3 |
third quartile value with quantile type 2, the SAS default |
Max |
maximum value |
Range |
range of finite numbers, maximum - minimum |
IQR |
interquartile range with quantile type 2, the SAS default |
MAD |
median absolute deviation |
VarLL |
lower confidence limit of variance |
VarUL |
upper confidence limit of variance |
Skewness |
skewness |
SkewnessSE |
standard error of skewness |
Kurtosis |
kurtosis |
KurtosisSE |
standard error of kurtosis |
GeometricMean |
geometric mean, calculated only when all given values are positive. |
GeometricCV |
geometric coefficient of variation in percent, calculated only when all given values are positive. |
Kyun-Seop Bae [email protected]
UNIV(lh)UNIV(lh)
F-test for the ratio of two groups' variances. This is similar to var.test except using the summarized input.
vtest(v1, n1, v0, n0, ratio=1, conf.level=0.95)vtest(v1, n1, v0, n0, ratio=1, conf.level=0.95)
v1 |
sample variance of the first (test, active, experimental) group |
n1 |
sample size of the first group |
v0 |
sample variance of the second (reference, control, placebo) group |
n0 |
sample size of the second group |
ratio |
value for the ratio of variances under the null hypothesis |
conf.level |
confidence level |
For the confidence interval of one group, use the UNIV function.
The output format is very similar to that of var.test.
Kyun-Seop Bae [email protected]
vtest(10.5^2, 3529, 8.9^2, 5190) # NEJM 388;15 p1386 vtest(2.3^2, 13, 1.5^2, 11, conf.level=0.9) # Red book p240vtest(10.5^2, 3529, 8.9^2, 5190) # NEJM 388;15 p1386 vtest(2.3^2, 13, 1.5^2, 11, conf.level=0.9) # Red book p240
This is shown in SAS PROC REG as the Test of First and Second Moment Specification.
WhiteTest(rx)WhiteTest(rx)
rx |
a result of |
This is also called White's general test for heteroskedasticity.
Returns a direct test result by the more complex Theorem 2, not by the simpler Corollary 1.
Kyun-Seop Bae [email protected]
White H. A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity. Econometrica 1980;48(4):817-838.
WhiteTest(lm(mpg ~ disp, mtcars))WhiteTest(lm(mpg ~ disp, mtcars))
Write the standard R ‘MD5’ manifest into an installed package directory so that
checkMD5sums (and IQLM) can verify file integrity.
The manifest records the md5 checksum of every installed file. It is the same
‘MD5’ file that CRAN ships and that R's own installer writes; packages installed
from CRAN already contain it, but a local source install usually does not (which is
why the IQ integrity check then reports WARN). Run writeMD5LM() once,
immediately after installing the package, to record the trusted baseline; thereafter
checkMD5sums() detects any later modification of the installed files.
writeMD5LM(pkg = "sasLM", lib.loc = NULL)writeMD5LM(pkg = "sasLM", lib.loc = NULL)
pkg |
name of the installed package whose ‘MD5’ manifest should be written. |
lib.loc |
character vector of library paths to search for |
The manifest is written in the format used by checkMD5sums: one line per file,
<md5sum> *<relative-path>. The ‘MD5’ file itself is excluded. The
function uses only the exported md5sum. The installed package
directory must be writable.
Invisibly, the path to the ‘MD5’ file that was written. Called for its side effect.
Kyun-Seop Bae <[email protected]>
#writeMD5LM("sasLM") #sasLM::IQLM() # the file-integrity check now reports PASS#writeMD5LM("sasLM") #sasLM::IQLM() # the file-integrity check now reports PASS
This is similar to the two groups t-test, but using the standard normal (Z) distribution.
ztest(m1, s1, n1, m0, s0, n0, conf.level=0.95, nullHypo=0)ztest(m1, s1, n1, m0, s0, n0, conf.level=0.95, nullHypo=0)
m1 |
mean of the first (test, active, experimental) group |
s1 |
known standard deviation of the first group |
n1 |
sample size of the first group |
m0 |
mean of the second (reference, control, placebo) group |
s0 |
known standard deviation of the second group |
n0 |
sample size of the second group |
conf.level |
confidence level |
nullHypo |
value for the difference of means under the null hypothesis |
Use this only for known standard deviations (or variances) or very large sample sizes per group.
The output format is very similar to that of t.test.
Kyun-Seop Bae [email protected]
ztest(5.4, 10.5, 3529, 5.1, 8.9, 5190) # NEJM 388;15 p1386ztest(5.4, 10.5, 3529, 5.1, 8.9, 5190) # NEJM 388;15 p1386