d_labes ★★★ Berlin, Germany, 20090417 12:06 Posting: # 3554 Views: 61,992 

Replicate berlin buddy bears Dear all, especially the Ruses among you, (ok, the category could also has been Regulatives/guidelines) continuing my deeper look into R exemplified by the source code of bear I noticed the following: The key statement in bear for evaluating replicate BE studies (astonishing enough claimed as only 2sequence designs allowed?) is model < lme(y~seq+prd+drug, random=~1subj, data=blabla, method="REML") with y the PK parameter under evaluation (raw or logtransformed). This is the equivalent of SAS Proc MIXED code: Proc MIXED data=BlaBla ; The used model assumes equal within/between variabilities for Test and Reference and also no subjectbytreatmentinteraction, if I'm right. This is very stringent and consequent, if we talk about ABE. But in real world there are regulators out there! And some of them do not like consequent scientific assumptions / models. Thus it happens that the FDA model for evaluating ABE within replicate designs (discussed occasionally here in the forum, SEARCH!!!) must have different within/between variabilities for Test and Reference and also subjectbytreatmentinteraction included in the model. After reading again tons of Websites (seems nobody evaluates xover studies with R) I came up with the following attempt in R: model2 < lme(y~ drug + prd + seq, This code gives for the bear built in dataset the same covariance parameters as the SAS code from the FDA statistical guidance except the fact that SAS Proc MIXED uses variances and the lme() code uses standard deviations.  snip  SAS:
The Mixed Procedure R (have changed drug to tmt, prd to period, subj to subject) cat("G matrixn") any body out there who knows how to get more decimals here?
Residual 0.00335726  end snip  I was soooo happy: An old SAS dog, nothing learned else, has become a Ruser! But then it came to the freedom to some degree:  snip  SAS: now with the DDFM=SATTERTHWAITE (DDFM = denominator degree of freedom method) Type 3 Tests of Fixed Effects R: anova(model2,type="marginal")[2:4,]
 end snip  Ok, I have learned (by reading again tons of websites) that R has no package with Satterthwaite denominator degrees of freedom in mixed models. But all other available DDFM in Proc MIXED gave not a denDF=38 for the treatment effect. One gets always lower denDF in SAS and as a consequence of that wider 90% confidence intervals! Thus my choice of software would be lme() and R because with that it is easier to make my sponsors happy . But the question is: Are regulators evenly happy with that? Whats going on here? Do I miss anything in the R approach? Or is the denominator DF question also like type III SumOfSquares? (I have read some Websites indicating that). BTW: Parallel Bears (crossing at infinity according to the geometric axioms) have also same variance . But this is another story. — Regards, Detlew 
ElMaestro ★★★ Belgium?, 20090417 14:50 @ d_labes Posting: # 3556 Views: 56,814 

Hi dlabes, » Or is the denominator DF question also like type III SumOfSquares? (I have read some Websites indicating that). As far as I recall, the R default with no option to change it is to use a denDF which is similar to SAS' "containment" meffud. I you insist on having Satterthwaite's denDF or whatever, then R will not be your package of choice. One little problem is that FDA indirectly endorsed SAS some time ago, so many regulators (incl. EU) expect an output which can be replicated with SAS using FDA scripts. Then clearly R will fall short of the goal. I guess you can safely conclude this is another "type III SS" or LSMean phenomenon. Big bro dictates the standard and regulators have no clue. Do it the SAS way or do it wrong. OK, that was a provocation, sorry. One can go on and on; discussion of the proof that REML gives unbiased coefficients in unbalanced datasets comes next, I guess? Best regards EM. 
d_labes ★★★ Berlin, Germany, 20090417 15:07 @ ElMaestro Posting: # 3557 Views: 56,800 

Hi ElMaestro, » As far as I recall, the R default with no option to change it is to use a denDF which is similar to SAS' "containment" meffud. Here are the SAS results with DDFM=CONTAIN Type 3 Tests of Fixed Effects Seems not to contain the R results , as already mention in my previous post. — Regards, Detlew 
ElMaestro ★★★ Belgium?, 20090418 03:45 @ d_labes Posting: # 3560 Views: 56,932 

Dear dlabes » I was soooo happy: An old SAS dog, nothing learned else, has become a Ruser! Amazing Grace, how sweet the sound, That saved a wretch like me. I once was lost but now am found, Was blind, but now I see. (John Newton, 17251807) You did the right thing by parking the power in the corner. It must have been a tough decision. Now go bash it with a baseball bat, bury it in your back yard and you will have my deepest respect. Sincerely, EM. 
d_labes ★★★ Berlin, Germany, 20090420 08:42 @ ElMaestro Posting: # 3569 Views: 56,663 

Dear ElMaestro » [...] Now go bash it with a baseball bat, bury it in your back yard and you will have my deepest respect. what a political incorrect suggestion ! It is the "Power to know"! — Regards, Detlew 
yjlee168 ★★ Kaohsiung, Taiwan, 20090419 21:38 @ d_labes Posting: # 3567 Views: 56,833 

Dear D. Labes, » » Replicate berlin buddy bears very cute. I like that! » The key statement in bear for evaluating replicate BE studies (astonishing enough claimed as only 2sequence designs allowed?) is I think we can easily program bears to handle ABE data from more than 2 seq and/or more than 2 tmt in the future. » [...] » tmtR 0.14901302 tmtR » tmtT 0.12216340 0.15 any body out there who knows how to get more decimals here?You can try "options (digits=n)" in R console to increase decimals, where n the number decimal you like to display. If not work, it can be the builtin option in the package. » [...] But all other available DDFM in Proc MIXED gave not a denDF=38 for the treatment effect. One gets always lower denDF in SAS and as a consequence of that wider 90% confidence intervals! But as I can remember, we had validated 90%CI obtained from bears with those from SAS with replicate study using lme(). I will post the results here later. Basically we got the equivalent 90%CIs with SAS. Indeed, we know DenDFs obtained from bears are different from those using SAS. However, 90% CIs are the same. — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
d_labes ★★★ Berlin, Germany, 20090420 08:30 @ yjlee168 Posting: # 3568 Views: 56,654 

Dear bears, dear all » You can try "options (digits=n)" in R console to increase decimals, where n the number decimal you like to splay. If not work, it can be the builtin option in the package. Of course I had tried this. Unfortunately all this doesn't help here. » But as I can remember, we had validated 90%CI obtained from bears with those from SAS with replicate study using lme(). I will post the results here later. Basically we got the equivalent 90%CIs with SAS. Indeed, we know DenDFs obtained from bears are different from those using SAS. However, 90% CIs are the same. Don't worry about this, because SAS and bear gave the same result (also with DDFM=SATTERTHWAITE or KR) as long as lme(y~tmt+period+sequence, random=1subject, data=blabla) and the Proc MIXED equivalent of that is used. The difference comes only into effect if one tries to implement the FDA model (which is mandatory I think, at least for studies aimed to the FDA regulated market). My questions remains: Do I do something wrong here? If yes, what? If not, what are the chances that regulators accept the bear / R handling of replicate studies? — Regards, Detlew 
yjlee168 ★★ Kaohsiung, Taiwan, 20090421 14:03 @ d_labes Posting: # 3575 Views: 56,687 

dear D. Labes, » My questions remains: » Do I do something wrong here? If yes, what? I don't know why you used TYPE=UN in your SAS code. In the FDA 2001 Guidance (Appendix E, p.37), it says: "...In the Random statement, TYPE=FA0(2) could possibly be replaced by TYPE=CSH. This guidance recommends that TYPE=UN not be used, as it could result in an invalid (i.e., not nonnegative definite) estimated covariance matrix..." » If not, what are the chances that regulators accept the bear / R handling of replicate studies? As long as bear complies with the FDA Guidance (2001) and is statistically correct, it should be acceptable, I guess. Edit: Link corrected for FDA's new site. [Helmut] — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
d_labes ★★★ Berlin, Germany, 20090421 16:11 @ yjlee168 Posting: # 3578 Views: 56,646 

Dear Yungjin, » I don't know why you used TYPE=UN in your SAS code. [...] This was because I could not identify the syntax to implement such things (FA0(2) or CSH) within lme() and I had the guess that the lme() output "Random effects: corresponds eventually to SAS UN parameterization. BTW could not find any hint up to now what this really means: General positivedefinite, LogCholesky parametrization. Questions for that are always answered: Buy the book of Pinheiro/Bates "MixedEffects Models in S and SPLUS", Springer (2000), there it is described. But recent I am a little short in money. We have a financial crisis . » [...] As long as bear complies with the FDA Guidance (2001) and is statistically correct, it should be acceptable, I guess. Sure! As long as ... But the FDA statistical guidance explicitly states the model with different inter/intrasubject variabilities (Proc MIXED code in appendix G, I think) in evaluating ABE in replicate studies. And this is not the model implemented at present in bear for evaluation of replicate studies. I myself are convinced that the assumptions of no subjectbytreatment interaction and equal within subject variabilities are reasonable if we talk about ABE. Although with these assumptions an advantage of replicate designs vanishes, namely that we are able to study different variabilities. But unfortunately statistically correct or reasonable is in regulators view not always correct or reasonable. Not so astonishing because statisticians are them selfs of such different opinions if I think about Type III SumOfSqares or now the question of denominator degrees of freedom. — Regards, Detlew 
Helmut ★★★ Vienna, Austria, 20090422 00:46 @ d_labes Posting: # 3580 Views: 57,095 

Dear D. Labes, » [...] BTW could not find any hint up to now what this really means: General positivedefinite, LogCholesky parametrization. Questions for that are always answered: Buy the book of Pinheiro/Bates "MixedEffects Models in S and SPLUS", Springer (2000), there it is described. I have it in my shelf, but — Cheers, Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
d_labes ★★★ Berlin, Germany, 20090422 11:20 @ Helmut Posting: # 3581 Views: 56,603 

Dear Helmut, thanks very much for the links. » [ I.] too stupid to understand >90% of it's contents » [II.] too short of time to dig into it […] Me too. But my percentage for [I.] is higher . Much to much "matratzen" in it for an Initiated. — Regards, Detlew 
yjlee168 ★★ Kaohsiung, Taiwan, 20090423 08:08 @ d_labes Posting: # 3587 Views: 56,702 

Dear D. Labes, Thank you for your message. I find that currently bear's model dealing with replicate ABE data is equivalent to SAS Proc MIXED as the follows: PROC MIXED; (Ref.: Jones, B and Kenward MG. Design and Analysis of CrossOver Trials, 2nd ed., Chapman & Hall/CRC, 2003). If using TYPE setting, either FA0(2) or CSH as recommended by FDA, we will not get the exactly same 90%CIs from running bear using lme(). I guess we need to find out how to comply with here. — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
d_labes ★★★ Berlin, Germany, 20090423 09:48 @ yjlee168 Posting: # 3588 Views: 56,835 

Dear Yungjin, where does the code come from? I cannot find it in Jones / Kenward. And sorry it does not work for me (Syntax errors in the random statement! Maybe there is something missing, Intercept or sub(seq)?). And if you have working code, please guide me to the different intrasubject variabilities (resulting always from the SAS REPEATED statement) in the results from bear's lme() statement. Up to now I have not found any equivalent in bear's output. And what about the DDFM=blabla? I have learned so far (in reading tons of Web sites) that there is no such thing as different denominator degrees of freedom methods in lme()? I had been informed further, that the successor of lme()  lmer()  does not give any DenDF and pvalues for the effects in the model (fixed or random) by design, because the inventor and developer D. Bates does not belief that the test statistics follow any of the usually applied distribution (F distr.). See here f.i. A consequent but strange countenance. Since SAS Proc MIXED uses these DDFM also in calculation of CI's for the treatment effect, there must be any difference. BTW: Kenward/Jones in their Proc MIXED code always use DDFM=KR, i.e. the Kenward/Roger method. Not so astonishing I think . BTW2: Let me state explicitly that these discussions are not for proving you / bear wrong. I am seriously on the way to change "The power to know" to R, but only if I can fulfill the needs in statistics for bioequivalence studies (regulators bullet proof), my everyday bread and butter work. And currently I suspect that this is not feasible for replicate studies. Or is anybody out there to prove me that black is white? Anybody out there who has evaluated a replicate study with R (not necessarily bear) which was approved by a regulatory agency? — Regards, Detlew 
yjlee168 ★★ Kaohsiung, Taiwan, 20090423 10:46 (edited by yjlee168 on 20090423 12:27) @ d_labes Posting: # 3590 Views: 56,729 

Dear D. Labes, Sorry about the code. Here is the correct one. proc mixed data=blabla; I've changed the parameter names to fit our discussion. » where does the code come from? » I cannot find it in Jones / Kenward. Ch. 7. » And what about the DDFM=blabla? I have learned so far (in reading tons of Web sites) that there is no such thing as different denominator degrees of freedom methods in lme()? Nope, we cannot find equivalent one, either. » BTW: Kenward/Jones in their Proc MIXED code always use DDFM=KR, i.e. the Kenward/Roger method. Not so astonishing I think . Yes, you're right. » BTW2: Let me state explicitly that these discussions are not for proving you / bear wrong. I am seriously on the way to change "The power to know" to R, but only if I can fulfill the needs in statistics for bioequivalence studies (regulators bullet proof), my everyday bread and butter work. Wow. You can choose both, can't you? If you can't, better stay with Power to know for now. » And currently I suspect that this is not feasible for replicate studies. Or is anybody out there to prove me that black is white? no luck so far. — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
d_labes ★★★ Berlin, Germany, 20090423 15:07 @ yjlee168 Posting: # 3594 Views: 56,609 

Dear Yungjin, » » where does the code come from? » » I cannot find it in Jones / Kenward. » » Ch. 7. Ok, this is the code for a classical 2x2x2 crossover, as I had stated above in my opening of this thread. » Wow. You can choose both, can't you? If you can't, better stay with Power to know for now. Of course I could use both. But why? Why should an old lazy SAS dog take the extra effort in learning R (or any other Statistical system) if he had to stay with SAS for part of his work ? And if he had all things he needed in his everyday work already coded in that infamous "Power to know". — Regards, Detlew 
Helmut ★★★ Vienna, Austria, 20090423 15:25 @ d_labes Posting: # 3595 Views: 56,662 

Dear D. Labes and all bears! » Why should an old lazy SAS dog take the extra effort in learning R (or any other Statistical system) if he had to stay with SAS for part of his work ? Select at least one:
— Cheers, Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
d_labes ★★★ Berlin, Germany, 20090423 16:06 @ Helmut Posting: # 3596 Views: 56,531 

Dear Helmut! » Select at least one: My favorite one is Lorenzian. Unfortunately I have not heard about earlier. Now I argue it is to late for me to stay young . Additionally unfortunately it was my habit not to have a breakfast. My second one is: Masochism — Regards, Detlew 
Helmut ★★★ Vienna, Austria, 20090423 16:30 @ d_labes Posting: # 3597 Views: 56,718 

Dear D. Labes! » My favorite one is Lorenzian. Yeah, I like it also very much. » Unfortunately I have not heard about earlier. I was quoting according to my vanishing mem_{ory} from an Interview by Franz Kreuzer aired by Austrian (ORF) and German TV (ZDF) in May 1980. It was repeated on 01 March this year on 3SAT – unfortunatelly I missed it… The entire text was published in 1981 (Konrad Lorenz, Franz Kreuzer: Leben ist Lernen, Piper, ISBN 3492102239). I found a PDF on the net (don't ask me where, it's illegal) and could not locate the quote I gave… So I've heard the quote somewhere else… » Now I argue it is to late for me to stay young . Come on! Another quote from the mentioned book: <div lang="de" xml:lang="de"> K: Wie ertragen Sie das Altern? Wie haben Sie gemerkt, daß Sie zu altern begannen, worunter leiden Sie, welche Gegenwehr leisten Sie? Gibt es irgendeinen Gewinn, der die Bürde des Alterns leichter macht? — Cheers, Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
yjlee168 ★★ Kaohsiung, Taiwan, 20090423 19:56 @ d_labes Posting: # 3598 Views: 56,576 

Dear D. Labes, » Ok, this is the code for a classical 2x2x2 crossover, as I had stated above in my opening of this thread. Yes, but just let you know that we used the SAS code to validate bear when doing replicate data analysis. So... it is not appropriate to do so? » Of course I could use both. But why? As "the Power to know ", everything should be easy for him to get to know it. We all can see and appreciate your sincere and longterm contributions to this Forum. So why not keep The Power to Know? And R? why not? You just use two or more statistical languages to speak BE/BA. That's why this Forum is so amazing and so interesting because of the members like you and many others. » Why should an old lazy SAS dog take the extra effort in learning No switch, please. Be a SAS old dog AND also a R new dog. Thank you so much. — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
Helmut ★★★ Vienna, Austria, 20090423 20:03 @ yjlee168 Posting: # 3599 Views: 56,591 

Dear Yungjin, what do you think about a special category bear of the forum? — Cheers, Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
yjlee168 ★★ Kaohsiung, Taiwan, 20090423 20:19 @ Helmut Posting: # 3600 Views: 56,622 

Dear Helmut, Thank you for your kindly offering. How about the category R for BE/BA? I don't think bear should be the category name because of there are so many R tools applied in BE/BA or appearing in this Forum. Bear is just a small piece of them. Then we can expect or encourage more users to use R in the community in the future. What do you think? Many thanks. — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
Helmut ★★★ Vienna, Austria, 20090423 20:37 @ yjlee168 Posting: # 3601 Views: 56,641 

Dear Yungjin, yeah that's a good idea. Up to new horizons! I moved all respective threads from the category Software. If you want to set a link at your homepage to the new category, use https://forum.bebac.at/?category=19. — Cheers, Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
yjlee168 ★★ Kaohsiung, Taiwan, 20090424 09:19 @ Helmut Posting: # 3604 Views: 56,514 

Dear Helmut. » If you want to set a link at your homepage to the new category, use https://forum.bebac.at/?category=19. Done! A very nice and informative category. Thank you so much. — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
d_labes ★★★ Berlin, Germany, 20090424 08:27 @ yjlee168 Posting: # 3603 Views: 56,382 

Dear Yungjin, » » Ok, this is the code for a classical 2x2x2 crossover, as I had stated above in my opening of this thread. » Yes, but just let you know that we used the SAS code to validate bear when doing replicate data analysis. So… it is not appropriate to do so? If we talk about FDA view, definitely NO. If we talk about reference scaling in ABE also definitely NO. If we talk about ABE in the strictest sense, Ok. — Regards, Detlew 
yjlee168 ★★ Kaohsiung, Taiwan, 20090424 09:26 @ d_labes Posting: # 3605 Views: 56,454 

Dear D. Labes, » If we talk about FDA view, definitely NO. » If we talk about reference scaling in ABE also definitely NO. » If we talk about ABE in the strictest sense, Ok. Thank you for your message. We will continue to search and figure out how to do REML and DDFM (as SAS Proc MIXED & Random options) with R to comply with FDA Guidance. — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
ElMaestro ★★★ Belgium?, 20100603 05:21 @ d_labes Posting: # 5419 Views: 55,873 

Ahoy d_labes, there was still some low hanging fruit to pick in this old thread... » >getVarCov(model2) » G matrix » Random effects variance covariance matrix » tmtR tmtT » tmtR 0.0222050 0.0027286 » tmtT 0.0027286 0.0149240 » Standard Deviations: 0.14901 0.12216 (and blah) » Random effects: » Formula: ~tmt  1  subject » Structure: General positivedefinite, LogCholesky parametrization » StdDev Corr » tmtR 0.14901302 tmtR » tmtT 0.12216340 0.15 any body out there who knows how to get» more decimals here? Corr(A,B)=cov(A,B)/sqrt(Var(A)*Var(B)) Here = 0.0027286/(sqrt(0.0222050*0.0149240)), I guess the object returned by getVarCov() contains all the further decimals anyone can dream of. Best regards EM. 
yjlee168 ★★ Kaohsiung, Taiwan, 20100603 08:37 (edited by yjlee168 on 20100603 12:23) @ ElMaestro Posting: # 5421 Views: 55,926 

Dear ElMaestro, Great. Your post reminds me to modify R codes in bear, as suggested in the very first post of this thread (adding option weights=varIdent(form = ~ 1  drug)) by D. Labes. Thank you so much. » Corr(A,B)=cov(A,B)/sqrt(Var(A)*Var(B)) » Here = 0.0027286/(sqrt(0.0222050*0.0149240)), I guess the object returned by getVarCov() contains all the further decimals anyone can dream of. It sounds possible only if can we modify the codes for lme() originated from nlme package. The d_labes' output was obtained from lme() directly. note: the question for d_labes was deleted. Sorry about this. Forget how to use bbcode of strike out. It should be used such as — All the best,  Yungjin Lee bear v2.8.7: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
ElMaestro ★★★ Belgium?, 20100604 01:59 @ yjlee168 Posting: # 5433 Views: 55,832 

Dear Bears, » It sounds possible only if can we modify the codes for lme() originated from nlme package. The d_labes' output was obtained from lme() directly. Yes, the output has a default number of decimals that is a little too low and which nobody seems to be able to change with easy commands. I meant to show how to get more decimals. Once done, I guess you can probably manipulate the lme object directly (e.g. manually replace "0.15" by "0.1498xxxxx") before printing. This does not require a modification of the lme() internals. Best regards EM. — This week's list of things I absolutely detest: Corona virus, the which function in R, WIAWIA interfaces for scanning under Windows 10, the Bee Gees, the smell of my fridge.Best regards, ElMaestro 
ElMaestro ★★★ Belgium?, 20110127 21:50 @ d_labes Posting: # 6502 Views: 54,799 

Dear dlabes and others, » After reading again tons of Websites (seems nobody evaluates xover studies with R) I came up with the following attempt in R: » model2 < lme(y~ drug + prd + seq, » random= ~ tmt1subj, » #different within variabilities » weights=varIdent(form = ~ 1  drug), » data=BlaBla, method="REML") I am starting to believe lme does in some cases not really reach the goal, so I switched to lmer (package lme4). Syntactically something like:
ElModel=lmer(logCmax ~ per+seq+trt + (trt1subj)) which achieves the same. It will occasionally depending on the dataset actually fit the data better (that is, to a higher likelihood and thus with better sigma estimates) than lme can. And also: Who's got a 3period dataset (ref. replicated, but not testreplicated) with a known s_{WR} that I can look at? Mucho gracias — This week's list of things I absolutely detest: Corona virus, the which function in R, WIAWIA interfaces for scanning under Windows 10, the Bee Gees, the smell of my fridge.Best regards, ElMaestro 
d_labes ★★★ Berlin, Germany, 20110128 09:02 @ ElMaestro Posting: # 6503 Views: 54,971 

Dear ElMaestro! » I am starting to believe lme does in some cases not really reach the goal, so I switched to lmer (package lme4). Syntactically something like: » ElModel=lmer(logCmax ~ per+seq+trt + (trt1subj)) » which achieves the same. Really Using the newer lmer() instead of lme() is another possibility, considered by me also before coming up with my R implementation attempt of the FDA Proc MIXED code. But it lacks the separate estimability of the intrasubject error variance for Test and Reference, I think. As far as I know (from the R help lists) there are no plans to introduce somefink like weights in lmer().Therefore it is not so useful in the framework of scABE, I think. I could be wrong, but ... — Regards, Detlew 
ElMaestro ★★★ Belgium?, 20110128 09:42 @ d_labes Posting: # 6505 Views: 54,902 

Ahoy d_labes, » But it lacks the separate estimability of the intrasubject error variance for Test and Reference, I think. As far as I know (from the R help lists) there are no plans to introduce somefink like weights in lmer(). Here are two partial outputs from lme and lmer: With lme: Logrestrictedlikelihood: 98.12818 (...) Random effects: Formula: ~trt  1  subj Structure: General positivedefinite, LogCholesky parametrization StdDev Corr trt1 0.149013060 trt1 trt2 0.122163414 0.15 Residual 0.007643612 With lmer: AIC BIC logLik deviance REMLdev 170.2 150.0 95.12 229.3 190.2 (...) Random effects: Groups Name Variance Std.Dev. Corr subj trt1 2.2218e02 0.1490555 trt2 1.4913e02 0.1221188 0.150 Residual 3.3117e05 0.0057548 Note the slightly better likelihood found by lmer. Changes quite a bit (%wise, I mean) on the residual sigma. — This week's list of things I absolutely detest: Corona virus, the which function in R, WIAWIA interfaces for scanning under Windows 10, the Bee Gees, the smell of my fridge.Best regards, ElMaestro 
d_labes ★★★ Berlin, Germany, 20110128 11:37 (edited by d_labes on 20110128 11:56) @ ElMaestro Posting: # 6506 Views: 54,773 

Old pirate, where are the intraindividual (error or residual) variances or sd for Test and Reference? You have only one residual. The trt1 and trt2 variances (and corr) in your output are the between subject covariance parameter. — Regards, Detlew 
ElMaestro ★★★ Belgium?, 20110128 12:14 @ d_labes Posting: # 6507 Views: 54,952 

Dear d_labes, » where are the intraindividual (error or residual) variances or sd for Test and Reference? » You have only one residual. » The trt1 and trt2 variances (and corr) in your output are the between subject covariance parameter. Interesting. Let us discuss for a moment the covariance matrix V (not Z and not G, and at this point also not discussing which syntax SAS applies for the given problem (leaves disussions about CSH, FA(0) out here)). V holds the key. Further, let us assume for simplicity we only look at the first subject which is the first 4x4 of V where the data are ordered as TRTR (i.e. per1per2per3per4 for the first subject; in the original dataset above they were ordered differently but the order in this context does not change the point): Data in per1 and per3 correlate within the subject as (s)he gets T (within variance VT). Data in per2 and per4 correlate within the subject as (s)he gets R (within variance VR). To this we add a common error variance on the diagonal (V0). Having no particular mathematical understanding, the first 4x4 of V I think will be: V0+VT 0 VT 0 Both of the models discussed above correspond to a fit with this V, as far as I can see. — This week's list of things I absolutely detest: Corona virus, the which function in R, WIAWIA interfaces for scanning under Windows 10, the Bee Gees, the smell of my fridge.Best regards, ElMaestro 
d_labes ★★★ Berlin, Germany, 20110128 14:12 (edited by d_labes on 20110128 15:32) @ ElMaestro Posting: # 6508 Views: 54,852 

Dear ElMaestro, I think your V needs to be corrected to reflect the FDA model for replicate crossover designs. See Patterson et al "REPLICATE DESIGNS AND AVERAGE, INDIVIDUAL, AND POPULATION BIOEQUIVALENCE" GSK BDS Technical Report 2002 – 01 online resource It has to be composed from these two parts: (sequence TRTR, hope I got all these nasty subscripts right ) betweensubject:
and withinsubject diagonal matrix:
I'm not so familiar with the socalled V matrix but I guess its the sum of both?^{*see PS} Thus you have 5 variancecovariance parameters to estimate:
v_{BT}, v_{BT} and cov_{TR}, latter here given as Corr(elation) and one residual error term which is presumable the mean of v_{WT} and v_{WR} or other spoken common withinsubject error (assuming v_{WT} = v_{WR}=v_{e}). Have a look at the very beginning of this thread to see that FDA's Proc MIXED and my lme() attempt are giving the same covariance parameter values (except that Proc MIXED gives variances/covariances and lme() stddev's for the between part and an error term and factors for T/R for the withinsubject part, i.e. you have to calculate f.i. s_{WR} as error(residual) term multiplied by the factor for R. ^{*}PS: Regarding the matratzes we had it already here. [edit] link to technical report corrected — Regards, Detlew 
ElMaestro ★★★ Belgium?, 20110128 15:20 @ d_labes Posting: # 6509 Views: 54,766 

Dear d_labes, thanks for spending time to educate me. You might be right as you usually are, but it makes no sense to me. Let us for the time being forget the covariance of T and R. I am not disputing there should be a term for it, though. With the intrasubject s_{WR} and s_{WT} only on the diagonal of V my alarm bells are ringing. Imagine a refreplicated design; following your thinking we could fit an s_{WT} on the diagonal of V even though Test isn't replicated. This sounds very wrong. Apparently others have got fits and s_{WT} values for 2trt, 3seq, 3per designs (see e.g. here). This is fishy business. V can be considered as ZGZ^{T} plus R for the model afficionados , but we need not consider it here. The intrasubject issue can be modeled in either one of them, as long as the resulting V is the right one. The link you pasted doesn't work here. I'd be happy to read more about it though. Edit: Link corrected (David changed the domain name and folder from www.boomer.org/pkin/ to www.pharmpk.com/) [Helmut] — This week's list of things I absolutely detest: Corona virus, the which function in R, WIAWIA interfaces for scanning under Windows 10, the Bee Gees, the smell of my fridge.Best regards, ElMaestro 
d_labes ★★★ Berlin, Germany, 20110128 16:10 @ ElMaestro Posting: # 6510 Views: 54,751 

Dear Öberster Größter Meister, » With the intrasubject s_{WR} and s_{WT} only on the diagonal of V my alarm bells are ringing. Imagine a refreplicated design; following your thinking we could fit an s_{WT} on the diagonal of V even though Test isn't replicated. This sounds very wrong. Apparently others have got fits and s_{WT} values for 2trt, 3seq, 3per designs (see e.g. here). This is fishy business. Well noticed . We had discussed this peculiarity of the partial replicate design also here in the forum if you remember. See f.i. this post and followings. But be warned: Rather lengthy and nitpickin' . As you can see there we get s^{2}_{WT}, but rather different values between different software. And, as I have noticed, you get also different values if you change optimization parameters in SAS Proc MIXED. Seems the REML optimisation stops by chance. Thus that value is not real, the model is overspecified, but the REML or ML estimation can't of course decide how the model has to be formulated without an explicit value for it. It's the same as if we attempt to fit a 2x2 design with distinct intrasubject error terms for Test and Reference. We could not be successful because only (s^{2}_{WT}+s^{2}_{WR})*0.5 is estimable. Meanwhile I have the opinion that we can't handle this evaluation via such a detailed model of the different sources of variation. I at any rate was not successful to implement a reasonable simpler model in Proc MIXED. Maybe this is one reason why the FDA code given in the progesterone guidance resorts to evaluation via intrasubject contrast? BTW: Link corrected. — Regards, Detlew 