Bioequivalence and Bioavailability Forum

ElMaestro
★★★

Denmark,
2017-10-07 01:01
(2772 d 22:17 ago)

Posting: # 17871
Views: 19,339

Goodness of fits: one model, different datasets [General Statistics]

Hi all,

I have some different datasets that are not large at all, and I am doing some model fits on them and from that I am extracting a correlation measure like r squared etc. For the sake of simplicity let us just assume there in a single model like a linear regression without weights.

I would like to compare apples and pears, or at least the fits of three different datasets A, B and C. These datasets do not have the same number of data points, so I cannot fairly compare e.g. r squared of A vs B vs C directly.

Akaike and Schwarz are presumably not the way to go, I think, as I am not varying the model but the dataset, so to say. Kolmogorov-Smirnoff would potentially be useful if I had a boatload of points, which I don't anyway. I am very poor at explaining what I think I am looking for :crying:

but I would call it a "fit likelihood" or "correlation statistic that is sample size corrected" :lookaround:

. Google and Wikipedia aren't my friends in this regard (although on all other matters, including politics, religion, science and baking recipes G. and W. are always providing the right answers).

Does anyone here know of a handy statistic that allows a fair comparison of goodness of fits across datasets with unequal sizes, given a single model??

Muchas gracias.

—
Pass or fail!
ElMaestro

nobody
nothing

2017-10-07 18:03
(2772 d 05:16 ago)

@ ElMaestro
Posting: # 17872
Views: 17,942

Goodness of fits: one model, different datasets

Post reply

Hi Danmark!

❝ ... I am very poor at explaining what I think I am looking for :crying: but I would call it a »"fit likelihood" or "correlation statistic that is sample size corrected"...

Ehhm, whut? Maybe you are to fixated for statical testing? How about some oldfashioned "visual inspection of fit" paired with "back-calculated value of calibration samples from calibration function, including relative deviation"?

Just saying ;-)

No idea what you want to compare in the end, so...

—
Kindest regards, nobody

ElMaestro
★★★

Denmark,
2017-10-07 20:06
(2772 d 03:13 ago)

@ nobody
Posting: # 17873
Views: 17,989

Experimental setup, details

Post reply

Hi nobody,

❝ Ehhm, whut? Maybe you are to fixated for statical testing? How about some oldfashioned "visual inspection of fit" paired with "back-calculated value of calibration samples from calibration function, including relative deviation"?

❝ Just saying ;-)

Your idea is absolutely good. But there are complications...

❝ No idea what you want to compare in the end, so...

Let us say I have a cattle prod, a long-haired Austrian, and a microphone that is rigged up to my computer for recording.

The cattle prod is a wonderful tool, it makes persuasion such an easy task. My cattle prod has a lot of knobs and buttons and sliders that allow me to tweak combinations of:
- Tesla coil zap modifier strength
- Evil discharge combobulator intensity
- Ion stream barbaric voltage gain
- Apocalyptic wolfram anode ray modulation
...and so forth.

I am zapping the Austrian is his sorry rear and recording the intensity (dB) of his immediate...shall we say.... verbal reaction. Let me add, each experiment takes no more than a few seconds but there is a bit of recovery time before the next experiment can be run and that isn't related to latency of the prod. Surprisingly, under the present experimental conditions my sensitive microphone also records a lot of cursing and swearing after each zap. I am not sure why or what the source is, but it is a matter I am looking into as I am not sure such behavioral excursions are good for my poor sensitive electronics.

Anyways, I know from the past that with the default prod settings there is a decent -though not overwhelmingly good- correlation of the reaction recorded in dB and the prod's output measured in GJ per zap. I have done that experiment with 16 evenly dispersed levels between 0 and 17500 GJ/zap and I got an r squared of 0.9583.

The question is if I can improve the correlation by playing a little around with the seetings mentioned above. So I acquired via Amazon an ACME Automated Experimental Prod Analyzer Kit (only $ 598 w. 2 years of warranty) so that I can automate the experiments comprising of changing a lot of setting, actuating the prod and recording the resulting reaction.

So I have generated 9629 datasets, some of which have 16 data points and some have 15, 14 or 13 because of missing values (in practice this means a relative lack of contact between the tip of the prod when actuated and the Austrian's butt).
I need to find the settings for which I am most likely to be able give a zap that triggers a 138 dB (+/- 0.1 dB) response. Hence I need the correlation to be as good as possible, including a correction for sample size. With 9629 datasets I cannot visually inspect every single set and determine the best fit. I must find a solution in R.

This is my grave predicament. "Help me noone kenobi, you are my only hope".

—
Pass or fail!
ElMaestro

ElMaestro
★★★

Denmark,
2017-10-07 21:07
(2772 d 02:11 ago)

@ ElMaestro
Posting: # 17874
Views: 18,103

Visualization

Post reply

Dear all,

for anyone interested in the setup please follow this link.

Around 0:42 into the educational journey you will see an early analog version of a cattle prod (actually back then the concenpt was only semi-developed so some may call it only a cat prod), and a few seconds later you will witness the natural consequence of the prod's encounter with it's intended biological target. I believe the reaction here amounts to only roughly 115 dB, arguably because of the relative lack of advanced settings and the choice of population.

—
Pass or fail!
ElMaestro

Helmut
★★★

Vienna, Austria,
2017-10-08 19:17
(2771 d 04:02 ago)

@ ElMaestro
Posting: # 17875
Views: 17,879

multiple regression?

Post reply

Hi ElMaestro,

nice setup. Given

❝ […] a lot of knobs and buttons and sliders that allow me to tweak combinations of:

❝ - Tesla coil zap modifier strength

❝ - Evil discharge combobulator intensity

❝ - Ion stream barbaric voltage gain

❝ - Apocalyptic wolfram anode ray modulation

❝ ...and so forth.

I’m I right that you are aiming at multiple regression (i.e., >1 regressor and and one regressand)? See there for an arsenal of methods. For the same number of regressors the adjusted R² takes the number of data points into account.

x1 <- rnorm(10) y1 <- x1*2+rnorm(10, 0, 0.1) muddle1 <- lm(y1 ~ x1) R2.adj1 <- summary(muddle1)$adj.r.squared x2 <- x1[-1] # drop 1 y2 <- y1[-1] # drop 1 muddle2 <- lm(y2 ~ x2) R2.adj2 <- summary(muddle2)$adj.r.squared cat(R2.adj1, R2.adj2, "\n")

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

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

ElMaestro ★★★ Denmark, 2017-10-08 19:30 (2771 d 03:48 ago) @ Helmut Posting: # 17876 Views: 17,893	just y=ax+b Post reply
	Hi Hötzi, I simply have a lot of regression lines each of the type y=ax+b. I just need to find the best one and the only complication is that not all lines have the same number of points. No multiple regression and no fancy stuff like weights or funky residual distributions. — Pass or fail! ElMaestro

Helmut
★★★

Vienna, Austria,
2017-10-08 19:35
(2771 d 03:43 ago)

@ ElMaestro
Posting: # 17877
Views: 17,865

just y=ax+b

Post reply

Hi ElMeastro,

❝ I simply have a lot of regression lines each of the type y=ax+b. I just need to find the best one and the only complication is that not all lines have the same number of points.

See my edit above. Do you have the raw data or only the parameters a₁, a₂, … a_i, b₁, b₂, … b_i, and the sizes? In the former case easy (max. R²_adj), in the latter no idea.

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

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

ElMaestro ★★★ Denmark, 2017-10-08 19:50 (2771 d 03:29 ago) @ Helmut Posting: # 17878 Views: 17,901	just y=ax+b Post reply
	Hi again, and thanks for looking into it. It is as simple as this: `x1=c(1,2,3,4,5) y1=c(10.5, 11.4, 12.6, 13.3, 14.6) x2=c(1,2,3,4,5) y2=c(10.3, 11.4, NA, 13.5, 14.6) L1=lm(y1~x1) L2=lm(y2~x2) summary(L1) summary(L2)` I would like to know how to judge if the correlation is quantitatively appearing to be best for the first line or for the second line. — Pass or fail! ElMaestro

nobody
nothing

2017-10-08 22:26
(2771 d 00:52 ago)

@ ElMaestro
Posting: # 17879
Views: 18,018

just y=ax+b

Post reply

❝ x1=c(1,2,3,4,5)

❝ y1=c(10.5, 11.4, 12.6, 13.3, 14.6)

❝

❝ x2=c(1,2,3,4,5)

❝ y2=c(10.3, 11.4, NA, 13.5, 14.6)

❝

13.3 vs. 13.5? By intention?

"RMSE

The RMSE is the square root of the variance of the residuals. It indicates the absolute fit of the model to the data–how close the observed data points are to the model’s predicted values. Whereas R-squared is a relative measure of fit, RMSE is an absolute measure of fit. As the square root of a variance, RMSE can be interpreted as the standard deviation of the unexplained variance, and has the useful property of being in the same units as the response variable. Lower values of RMSE indicate better fit. RMSE is a good measure of how accurately the model predicts the response, and is the most important criterion for fit if the main purpose of the model is prediction."

...educated guess, no idea if sensitive for number of calibrators :-D

http://www.theanalysisfactor.com/assessing-the-fit-of-regression-models/

..didn't read it, maybe you want to have a look:

http://orfe.princeton.edu/~jqfan/papers/01/ant.pdf

—
Kindest regards, nobody

yjlee168 ★★★ Kaohsiung, Taiwan, 2017-10-08 23:28 (2770 d 23:50 ago) (edited on 2017-10-09 08:35) @ ElMaestro Posting: # 17880 Views: 18,013	ANCOVA with R? Post reply
	Dear Elmaestro, I am not quite sure if ANCOVA is what you need. Here is a 10-min video (title: ANCOVA in R) from Youtube. — All the best, -- Yung-jin Lee bear v2.9.2:- created by Hsin-ya Lee & Yung-jin Lee Kaohsiung, Taiwan https://www.pkpd168.com/bear Download link (updated) -> here

DavidManteigas
★

Portugal,
2017-10-09 12:34
(2770 d 10:45 ago)

@ ElMaestro
Posting: # 17881
Views: 17,852

just y=ax+b

Post reply

Dear El Maestro,

That is changing the logic of statistics :-D

Find the right data for your model instead of finding the right model for your data.

For those models, I think R squared is the way to go, although different sample sizes might have an influence depending on the range of samples that you have. How "small" are the samples? Between 30 and 50 observations per sample? Would methods like cross-validation work on the datasets?
What about calculating confidence intervals for the adjusted R squared of each model and use a Forest plot of those CI, one line for each dataset? This might give you an idea of the impact of the different sample sizes in this particular problem and to decide whether just looking at R square would be enough. Depending on the number of datasets, visual inspection and influence statistics might help on deciding which is the best data for the model. :-D

Regards,
David

nobody
nothing

2017-10-09 12:45
(2770 d 10:33 ago)

@ DavidManteigas
Posting: # 17882
Views: 17,920

just y=ax+b

Post reply

Best guess: in case of unweighted linear regression the highest calibrators always have the highest impact. For your example: Take 138 dBA as the mean of the regression function, as predictions will be best arround the mean.

If there is heteroscedastic data and weighted regression, the calibrators arround the weighted mean will have least impact on the calibration function.

I still have not fully understood the scope of the survey. Is it on best practice of study planning for calibration experiments?

—
Kindest regards, nobody

Helmut
★★★

Vienna, Austria,
2017-10-10 20:15
(2769 d 03:03 ago)

@ ElMaestro
Posting: # 17885
Views: 17,873

just y=ax+b

Post reply

Hi ElMaestro,

❝ I would like to know how to judge if the correlation is quantitatively appearing to be best for the first line or for the second line.

Sumfink along this?

x1 <- c(1, 2, 3, 4, 5) y1 <- c(10.5, 11.4, 12.6, 13.3, 14.6) n1 <- length(x1) x2 <- c(1, 2, 3, 4, 5) y2 <- c(10.5, 11.4, NA, 13.3, 14.6) n2 <- length(x2) x3 <- c(1, 2, 4, 5) y3 <- c(10.5, 11.4, 13.3, 14.6) n3 <- length(x3) raw <- matrix(c(rep(1, n1), rep(2, n2), rep(3, n3), x1, x2, x3, y1, y2, y3, rep(n1, each=n1), rep(n2, each=n2), rep(n3, each=n3)), nrow=n1+n2+n3, ncol=4, byrow=FALSE, dimnames=list(NULL, c("set", "x", "y", "n"))) raw <- as.data.frame(raw); attach(raw) set <- unique(set) res <- data.frame(set=set, n=c(n1, n2, n3), NAs=NA, int=NA, slope=NA, RSE=NA, R2.adj=NA) for (j in seq_along(set)) { mod <- lm(y ~ x, data=raw, set==j) sum.mod <- summary(mod) res[j, "NAs"] <- sum(is.na(raw$y[raw$set==j])) res[j, "int"] <- coef(mod)[[1]] res[j, "slope"] <- coef(mod)[[2]] res[j, "RSE"] <- sum.mod$sigma res[j, "R2.adj"] <- sum.mod$adj.r.squared } print(raw, row.names=FALSE) print(res[with(res, order(-R2.adj, +n)), ], row.names=FALSE) detach(raw)

Fits linear models to the sets and gives results sorted by R²_adj (descending) and n (ascending). y-values of set 2 and 3 are identical. Shows that NA gives the same result as a completely missing row due to the default na.action=na.omit in lm().

set n NAs int slope RSE R2.adj 3 4 0 9.42 1.01 0.1565248 0.9928293 2 5 1 9.42 1.01 0.1565248 0.9928293 1 5 0 9.45 1.01 0.1494434 0.9912998

Edit: Toys for boys:

a <- 9.45 b <- 1.01 sets <- 16 n <- 5 x <- 1:n sd <- 0.25 pct <- 20 # percentage of NAs y <- NULL for (j in 1:sets) { y <- c(y, a + b*x + rnorm(n, 0, sd)) } raw <- matrix(c(rep(1:sets, each=n), rep(x, sets), y, rep(n, n*sets)), nrow=n*sets, ncol=4, byrow=FALSE, dimnames=list(NULL, c("set", "x", "y", "n"))) loc <- floor(runif(as.integer(pct*sets*n/100), min=1, max=n*sets)) raw[loc, "y"] <- NA raw <- as.data.frame(raw); attach(raw) set <- 1:sets res <- data.frame(set=set, n=rep(n, sets), NAs=NA, used=NA, int=NA, slope=NA, RSE=NA, R2.adj=NA) op <- par(no.readonly=TRUE) split.screen(c(4, 4)) # adjust for the number of sets! for (j in seq_along(set)) { mod <- lm(y ~ x, data=raw, set==j) sum.mod <- summary(mod) res[j, "NAs"] <- sum(is.na(raw$y[raw$set==j])) res[j, "used"] <- res[j, "n"] - res[j, "NAs"] res[j, "int"] <- coef(mod)[[1]] res[j, "slope"] <- coef(mod)[[2]] res[j, "RSE"] <- sum.mod$sigma res[j, "R2.adj"] <- sum.mod$adj.r.squared screen(j) par(mar=c(1.5, 1.5, 1, 0.1), mgp=c(2, 0.5, 0), line=2, tcl=-0.25) if (res[j, "used"] >= 2) { plot(raw$x[raw$set==j], raw$y[raw$set==j], las=1, xlab="", ylab="", ylim=c((a+b*x[1])*0.95, (a+b*x[n])/0.95), cex.axis=0.65, cex.main=0.65, font.main=1, main=paste("Rsq.adj =", signif(res[j, "R2.adj"], 6))) ifelse (res[j, "used"] == 2, col <- "red", col <- "blue") abline(coef=coef(mod), col=col) legend("topleft", cex=0.65, box.lty=0, x.intersp=0, legend=paste0("set ", j, ": ", res[j, "used"], " used")) } } res <- res[with(res, order(-R2.adj, +n)), ] screen(head(res$set, 1)) par(mar=c(1.5, 1.5, 1, 0.1)) legend("right", cex=0.75, box.lty=0, x.intersp=0, legend="\"best\"") screen(tail(res$set, 1)) par(mar=c(1.5, 1.5, 1, 0.1)) legend("right", cex=0.75, box.lty=0, x.intersp=0, legend="\"worst\"") close.screen(all=TRUE) par(op) detach(raw) print(res, row.names=FALSE)

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

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

Goodness of fits: one model, different datasets [General Sta­tis­tics]

Goodness of fits: one model, different datasets

Experimental setup, details

Visualization

multiple regression?

just y=ax+b

just y=ax+b

just y=ax+b

just y=ax+b

ANCOVA with R?

just y=ax+b

just y=ax+b

just y=ax+b

Goodness of fits: one model, different datasets [General Statistics]