parallel designs (Welch-Satterthwaite in R) [🇷 for BE/BA]

posted by Helmut Homepage – Vienna, Austria, 2007-04-18 19:05 (6635 d 18:02 ago) – Posting: # 675
Views: 12,667

Hi everybody!

I prepared some code in [image] (v2.4.1) to be used with an example data set (actually data from a 2×2 cross-over, but only data of period 1 are used).
Save the data-file to the bin folder below the main [image]-application (you may use any other location, but then you will have to tell [image] where to look for it (in [image]’s GUI: File > Change directory...).
  disply <- function() {
    cat("Test/Reference with 95% Confidence Limits (90% CI)", fill = TRUE)
      round(tbldiff, 4) }
# read in data
  resp <- read.csv("24par.txt", header = T)
# calculate natural logs
  resp$logAUC <- log(resp$AUC)
  relevel(resp$treatment, ref = "Reference" )
# simple plot, period 1 data only
  def.par <- par(no.readonly = TRUE) # save default
  layout(matrix(c(1,2), byrow = FALSE, ncol=2 ))
  plot(AUC ~ treatment, data = resp, subset = (period == 1), log = "y" )
  plot(logAUC ~ treatment, data = resp, subset = (period == 1 ))
  par(def.par)                       # reset to default
# two sample t-test, equal variances; NOT recommended (anticonservative)!
  result <- t.test(logAUC ~ treatment,
      data = resp,
      subset = (period == 1),
      var.equal = TRUE,
      conf.level = 0.90)
# original output in log-domain
  result
# extract results from list and presentation in untransformed domain
  tbldiff <- matrix(
      c(as.numeric(exp(diff(result$estimate))),
      sort(as.numeric(exp(-result$conf.int)))),
      byrow = TRUE, nrow = 1)
  dimnames(tbldiff) <- list("Ratio", c("Point Estimate", "Lower CL", "Upper CL" ))
  disply()

# two sample t-test (Welch-Satterthwaite), unequal variances
  result <- t.test(logAUC ~ treatment,
    data = resp,
    subset = (period == 1),
    var.equal = FALSE,  # note: This is the default in R!
    conf.level = 0.90)
  result
  tbldiff <- matrix(
    c(as.numeric(exp(diff(result$estimate))),
    sort(as.numeric(exp(-result$conf.int)))),
    byrow = TRUE, nrow = 1)
  dimnames(tbldiff) <- list("Ratio", c("Point Estimate", "Lower CL", "Upper CL" ))
  disply()


Results should be:
        Two Sample t-test

data:  logAUC by treatment
t = 1.1123, df = 22, p-value = 0.278
alternative hypothesis: true difference in means is not equal to 0
90 percent confidence interval:
 -0.09702467  0.45390930
sample estimates:
mean in group Reference      mean in group Test
               3.562273                3.383831

Test/Reference with 95% Confidence Limits (90% CI)
      Point Estimate Lower CL Upper CL
Ratio         0.8366   0.6351   1.1019


        Welch Two Sample t-test

data:  logAUC by treatment
t = 1.1123, df = 21.431, p-value = 0.2783
alternative hypothesis: true difference in means is not equal to 0
90 percent confidence interval:
 -0.0973465  0.4542311
sample estimates:
mean in group Reference      mean in group Test
               3.562273                3.383831

Test/Reference with 95% Confidence Limits (90% CI)
      Point Estimate Lower CL Upper CL
Ratio         0.8366   0.6349   1.1022


Have fun, and don’t trust in commercial software!
[image]

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Thread locked

Complete thread:

UA Flag
Activity
 Admin contact
23,424 posts in 4,927 threads, 1,671 registered users;
28 visitors (0 registered, 28 guests [including 20 identified bots]).
Forum time: 13:07 CEST (Europe/Vienna)

Truth and clarity are complementary.    Niels Bohr

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5