Medium rare. [🇷 for BE/BA]
Hi ElMaestro,
I prefer medium rare over well done (aka quick and dirty or clean and never)!
Below my more R-ish code. Changes:
Gives:
❝ Job done?
I prefer medium rare over well done (aka quick and dirty or clean and never)!
Below my more R-ish code. Changes:
- Reads CSV files in any format (from the web or local system, arbitrary column separator [default comma] and decimal separator [default period]).
- Vectorizes the data set.
- Checks the design and assigns a variable '
des
' ("partial
" or "full
") for further use. In principle all full and partial replicates are supported.
- Get rid of missing geometric means.
- In
Obj.F12()
returns a double instead of a 1:1 matrix. Useful if you decide to switch the optimizer (some require a vector).
- Direct access of D instead of
subset()
. Subset is convenient in development but leads to endless troubles when building a package.
optim(control = list(reltol = 1e-12, trace = 0,...)
- Information about optimization.
- G-matrix.
- I wonder why this code works at all.
InObj.F12()
the line
B <- -0.5*log(det(t(X) %*% solve(CovM) %*% X))
is the reason for the warnings becauset(X)
is not a square matrix. ThereforeB
is NaN.
- Adapt for full replicates (
varWT
already estimated).
- Adapt for the conventional TR|RT – mixed model also required for the FDA? Currently stops execution for such a data set.
- Most important: Get the PE (from package
lmerTest
together with Satterthwaite’s DFs) and calculate the CI.
- Implement the FDA’s RSABE (should be easy).
###############################
# Section 1: Household things #
###############################
GetData <- function(path, sep, dec) {
# read data from a CSV-file
if (missing(sep)) sep <- ","
if (!sep %in% c(",", ";", "t"))
stop("Column separator must be any of ',', ';', 't'.")
if (sep == "t") sep <- "\t"
if (missing(dec)) dec <- "."
if (!dec %in% c(".", ","))
stop("Decimal separator must be '.', or ','.")
# strip quotes and skip eventual commentary lines
D <- read.csv(path, sep = sep, dec = dec,
quote = "", comment.char = "#")
names(D) <- c("Subj", "Per", "Seq", "Trt", "Y")
cols <- c("Subj", "Per", "Seq", "Trt")
D[cols] <- lapply(D[cols], factor)
return(invisible(D))
}
CreateX <- function(D) {
# first the two treatments
TrtT <- as.numeric(as.character(D$Trt) == "T")
TrtR <- as.numeric(as.character(D$Trt) == "R")
X <- cbind(TrtT, TrtR)
R <- qr(X)$rank
# now the periods
for (p in unique(D$Per)) {
v <- as.numeric(D$Per == p)
XX <- data.frame(X, v)
names(XX)[ncol(XX)] <- paste0("P", p)
rnk <- qr(XX)$rank
if (rnk > R) {
X <- XX
R <- rnk
}
}
for (q in unique(D$Seq)) {
v <- as.numeric(D$Seq == q)
XX <- data.frame(X, v)
names(XX)[ncol(XX)] <- paste0("Q", q)
rnk <- qr(XX)$rank
if (rnk > R) {
X <- XX
R <- rnk
}
}
return(as.matrix(X))
}
Create.CovM <- function(Params) {
# block diagonal covariance matrix
varT <- Params[1]
varBR <- Params[2]
varWR <- Params[3]
varRT <- Params[4]
Nobs <- length(D$Y)
V <- matrix(0, ncol = Nobs, nrow = Nobs)
for (iRow in 1:Nobs) {
for (iCol in 1:Nobs) {
if (iCol == iRow) { # diagonal elements
if (D$Trt[iRow] == "T") {
V[iRow, iCol] <- V[iRow, iCol] + varT
} else {
V[iRow, iCol] <- V[iRow, iCol] + varWR + varBR
}
} else { # off diagonal elements
if (D$Subj[iRow] == D$Subj[iCol]) {
if (D$Trt[iCol] == D$Trt[iRow]) {
V[iRow, iCol] <- V[iRow, iCol] + varBR
} else {
V[iRow, iCol] <- V[iRow, iCol] + varRT
}
}
}
}
}
return(V)
}
############################
# Section 2: Matrix things #
############################
Obj.F12 <- function(Pars) {
# this is line 3 of page 10 of:
# http://people.csail.mit.edu/xiuming/docs/tutorials/reml.pdf
CovM <- Create.CovM(Pars)
A <- -0.5*log(det(CovM))
# TODO: NaNs produced in the next line cause t(X) is not a square matrix
B <- -0.5*log(det(t(X) %*% solve(CovM) %*% X))
est.b <- solve(t(X) %*% solve(CovM) %*% X) %*% t(X) %*% solve(CovM) %*% y
tmp <- y - X %*% est.b
C <- -0.5 *(t(tmp) %*% solve(CovM) %*% tmp)
C <- as.numeric(C) # make a scalar for other solvers
return(A+B+C)
}
########################
# Section 3: Execution #
########################
Some.Initial.Guesses <- function() {
# design (all possible ones, error handling, full or partial)
subjs <- unique(as.numeric(D$Subj))
trts <- sort(unique(as.character(D$Trt)), decreasing = TRUE)
if (!length(trts) == 2) stop("Only two treatments supported.")
if (sum(!trts %in% c("T", "R")) !=0)
stop("treatments must be coded as 'T' and 'R'.")
seqs <- sort(unique(as.character(D$Seq)), decreasing = TRUE) # T first
type <- paste(seqs, collapse = "|") # identifier
nseq <- length(unique(as.character(D$Seq)))
nper <- length(unique(as.numeric(D$Per)))
nsubj <- length(subjs)
if (nchar(type) == 19) { # 4-period 4-sequence full replicate designs
if (nper != 4) stop("4 periods required in this full replicate design.")
if (nseq != 4) stop("4 sequences required in this full replicate design.")
des <- "full"
}
if (nchar(type) == 9) { # 4-period full replicate designs
if (nper != 4) stop("4 periods required in this full replicate design.")
if (nseq != 2) stop("2 sequences required in this full replicate design.")
des <- "full"
}
if (nchar(type) == 7) { # 3-period replicates
if (type %in% c("TRT|RTR", "TRR|RTT")) {
if (nper != 3) stop("3 periods required in this full replicate design.")
if (nseq != 2) stop("2 sequences required in this full replicate design.")
des <- "full"
}
if (type == "TRR|RTR") {
if (nper != 3) stop("3 periods required in the extra-reference design.")
if (nseq != 2) stop("2 sequences required in the extra-reference design.")
des <- "partial"
}
}
if (nchar(type) == 11) { # Balaam's design or 3-sequence partial
if (!type == "TR|RT|TT|RR") { # Balamm's
if (nper != 3) stop("3 periods required in this partial replicate design.")
if (nseq != 3) stop("3 sequences required in this partial replicate design.")
des <- "partial"
} else { # TRT|RTR|RRT
if (nper != 2) stop("2 periods required in Balaam's design.")
if (nseq != 4) stop("4 sequences required in Balaam's design.")
des <- "full"
}
}
if (type == "TR|RT") stop("TR|RT design not supported yet.")
gm.T <- gm.R <- numeric(nsubj) # vector of geometric means
for (subj in seq_along(subjs)) {
gm.T[subj] <- mean(log(D$Y[which(D$Subj == subj & D$Trt == "T")]))
gm.R[subj] <- mean(log(D$Y[which(D$Subj == subj & D$Trt == "R")]))
}
# get rid of missings
gm.T <- gm.T[!is.na(gm.T)]
gm.R <- gm.R[!is.na(gm.R)]
if (des == "partial") { # guess varT (between + within)
m <- lm(log(Y) ~ Seq+Per, data = D[D$Trt == "T", ])
varT <- anova(m)["Residuals", "Mean Sq"]
} else { # guess VarWT
m <- lm(log(Y) ~ Seq+Per+Subj, data = D[D$Trt == "T", ])
varWT <- anova(m)["Residuals", "Mean Sq"]
}
# guess VarWR
m <- lm(log(Y) ~ Seq+Per+Subj, data = D[D$Trt == "R", ])
varWR <- anova(m)["Residuals", "Mean Sq"]
# guess varBR
varBR <- var(gm.R)-varWR
# TODO: full replicates
varRT <- 0.5*(varT+varBR)
rslt <- c(varT, varBR, varWR, varRT)
return(rslt)
}
MyREML <- function() {
Ini <- Some.Initial.Guesses()
F <- optim(par = Ini, fn = Obj.F12,
control = list(reltol = 1e-12, trace = 0, fnscale = -1))
Nms <- c("var_T", "var_bR", "var_wR", "covar_bTbR")
X <- data.frame(Nms, F$par, Ini)
names(X) <- c("Component", "Estimate", "Initial")
cat("Maximum :", F$value,
"\nEvaluations :", F$counts[["function"]],
"\nConvergence code:", F$convergence, "\n")
return(X)
}
# Data has to have five columns in this order:
# Subject (numeric), Period (numeric), Sequence (character),
# Treatment (character), Response (numeric)
# Column separator 'sep' comma (default), semicolon, or t (tab)
# Decimal separator 'dec' period (default) or comma
# path <- "https://bebac.at/downloads/ds01.csv" # EMA full replicate
path <- "https://bebac.at/downloads/ds02.csv" # EMA partial replicate
D <- NULL # nullify eventual previous one
D <- GetData(path)
X <- CreateX(D)
y <- log(D$Y)
# the workhorse
res <- MyREML() # patience!
G.matrix <- data.frame(c(NA, res$Estimate[4]),
c(res$Estimate[4], res$Estimate[2]))
names(G.matrix) <- 1:2
print(res, row.names = FALSE); cat("G matrix:\n"); print(G.matrix)
Gives:
Maximum : 75.98722
Evaluations : 197
Convergence code: 0
Warning message:
In log(det(t(X) %*% solve(CovM) %*% X)) : NaNs produced
Component Estimate Initial
var_T 0.07049482 0.06841615
var_bR 0.03621594 0.02776728
var_wR 0.01324652 0.01240137
covar_bTbR 0.04563735 0.04809171
G matrix:
1 2
1 NA 0.04563735
2 0.04563735 0.03621594
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz
The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz
The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Complete thread:
- Semireplicated + REML in R ElMaestro 2020-07-10 13:19
- Semireplicated + REML in R Helmut 2020-07-10 19:03
- Semireplicated + REML in R ElMaestro 2020-07-10 19:24
- Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57
- Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43
- Braveheart! Helmut 2020-07-11 15:38
- Who can help? ElMaestro 2020-07-12 12:28
- Who can help? ElMaestro 2020-07-12 13:26
- Update II ElMaestro 2020-07-12 21:36
- Update III ElMaestro 2020-07-12 21:46
- Final update today ElMaestro 2020-07-12 22:27
- Medium rare.Helmut 2020-07-13 13:52
- took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18
- Will take much more hours still… Helmut 2020-07-13 15:34
- Negative determinant ElMaestro 2020-07-14 03:22
- Are we loosers? Helmut 2020-07-14 13:58
- "we"? Loosers? ElMaestro 2020-07-14 15:07
- Misunderstanding? Helmut 2020-07-14 15:32
- "we"? Loosers? ElMaestro 2020-07-14 15:07
- Are we loosers? Helmut 2020-07-14 13:58
- took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18
- Medium rare.Helmut 2020-07-13 13:52
- Braveheart! ElMaestro 2020-07-13 10:13
- Braveheart! Helmut 2020-07-13 14:16
- Braveheart! PharmCat 2020-07-15 14:19
- Braveheart! ElMaestro 2020-08-02 17:39
- Who can help? ElMaestro 2020-07-12 12:28
- Braveheart! Helmut 2020-07-11 15:38
- Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43
- Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57
- Semireplicated + REML in R ElMaestro 2020-07-10 19:24
- We were all blind (except Detlew) Helmut 2020-07-15 14:27
- It is the opposite way around for me ElMaestro 2020-07-15 16:25
- Desultory thoughts Helmut 2020-07-15 17:33
- FDA RSABE is ISC d_labes 2020-07-15 18:13
- FDA RSABE is ISC Helmut 2020-07-16 11:11
- Desultory thoughts ElMaestro 2020-07-15 23:06
- Desultory thoughts Helmut 2020-07-16 10:59
- FDA RSABE is ISC d_labes 2020-07-15 18:13
- Desultory thoughts Helmut 2020-07-15 17:33
- Phoenix - which template? mittyri 2020-07-19 00:42
- FDA RSABE Project template_ v1.4.phxproj Helmut 2020-07-19 01:45
- It is the opposite way around for me ElMaestro 2020-07-15 16:25
- "By popular demand": likelihood ElMaestro 2020-07-24 10:07
- And by the way.... ElMaestro 2020-07-24 12:52
- And by the way.... PharmCat 2020-08-03 14:24
- Not understood ElMaestro 2020-08-03 22:55
- Not understood PharmCat 2020-08-05 01:41
- Not understood ElMaestro 2020-08-05 08:13
- Not understood PharmCat 2020-08-05 16:37
- Open issues ElMaestro 2020-08-06 21:11
- Open issues PharmCat 2020-08-07 00:02
- Open issues ElMaestro 2020-08-07 07:49
- Open issues PharmCat 2020-08-07 11:29
- Still can't make it work ElMaestro 2020-08-07 13:08
- Still can't make it work PharmCat 2020-08-07 15:42
- Still can't make it work ElMaestro 2020-08-07 16:44
- Still can't make it work PharmCat 2020-08-07 18:14
- Still can't make it work ElMaestro 2020-08-07 18:23
- And now it works ElMaestro 2020-08-07 21:31
- Still can't make it work PharmCat 2020-08-08 01:08
- Speed improvement ElMaestro 2020-08-08 12:25
- Speed improvement PharmCat 2020-08-08 17:27
- Speed improvement ElMaestro 2020-08-08 18:10
- Speed improvement PharmCat 2020-08-09 18:22
- Some tests... PharmCat 2020-08-10 11:48
- Speed improvement ElMaestro 2020-08-08 12:25
- Still can't make it work ElMaestro 2020-08-07 18:23
- Still can't make it work PharmCat 2020-08-07 18:14
- Still can't make it work ElMaestro 2020-08-07 16:44
- Still can't make it work PharmCat 2020-08-07 15:42
- Still can't make it work ElMaestro 2020-08-07 13:08
- Open issues PharmCat 2020-08-07 11:29
- Open issues ElMaestro 2020-08-07 07:49
- Open issues PharmCat 2020-08-07 00:02
- Open issues ElMaestro 2020-08-06 21:11
- Not understood PharmCat 2020-08-05 16:37
- Not understood ElMaestro 2020-08-05 08:13
- Not understood PharmCat 2020-08-05 01:41
- Not understood ElMaestro 2020-08-03 22:55
- And by the way.... PharmCat 2020-08-03 14:24
- And by the way.... ElMaestro 2020-07-24 12:52
- Semireplicated + REML in R Helmut 2020-07-10 19:03