Helmut ★★★ Vienna, Austria, 20150420 17:34 (2419 d 01:35 ago) Posting: # 14718 Views: 23,832 

Dear all, I closed this thread. Let’s continue over here. I would say that nowadays few people perform a replicated study without having referencescaling in mind. What would we need in bear?
PS: @ElMaestro. Remember that one of the referees of our reference datasetMS wanted to discuss the statistical model and we refused? Until somebody shows that what regulators want right now is crap (which will happen – at least partly) we should get the same results independent from the software used. I would be happy if a noncommercial one is amongst them. — Diftor heh smusma 🖖 Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
yjlee168 ★★★ Kaohsiung, Taiwan, 20150420 19:36 (2418 d 23:33 ago) @ Helmut Posting: # 14719 Views: 21,801 

Dear Helmut, » What would we need in bear? Wow! thank you so much. But that's too much. It must be a typo. May I suggest as "what we can do with R?" There are lots of R users/gurus in this Forum. » For the FDA that would be a mixedeffects model as given in the 2001 guidance and the progesterone guidance. I like this. » For the EMA we would need an “all effects fixed” model (Method A of the Q&A). One stupid question : if EMA needs an "all effects fixed" model, why do we use a mixed model to do analysis? I don't remember lme() right now, but lmer() must have at least a random variable. I tried lmer() last night and was pretty sure about that. Does EMA imply that we should use lm() to analyze replicate crossover dataset? » • FDA – RSABE: Again the progesterone guidance. » EMA – ABEL: The Q&A, Methods A & B. Crippled model to get CV_{wR}... » I would say, that’s the target. You should try to get funded from EMA or FDA for these projects. These are really great projects, indeed. — All the best,  Yungjin Lee bear v2.9.0: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
Helmut ★★★ Vienna, Austria, 20150421 01:02 (2418 d 18:07 ago) @ yjlee168 Posting: # 14725 Views: 21,675 

Hi Yungjin, » » For the EMA we would need an “all effects fixed” model (Method A of the Q&A). » if EMA needs an "all effects fixed" model, why do we use a mixed model to do analysis? Both the BEGL and the Q&A state that all effect should be fixed. EMA generally accepts a mixed model, but may ask for recalculation in borderline cases (e.g., extreme imbalance, CI close to the AR). Sometimes assessors ask for recalculation even if the CI was very narrow and the PE was close to 1. Happened to me once. The Q&A specifically mentions SAS’ Proc GLM , not Proc Mixed .» I don't remember lme() right now,… In lme() it should be possible to specify all effects fixed. Actually the random effect(s) are optional. If is doesn’t work, lm() should do.» Does EMA imply that we should use lm() to analyze replicate crossover dataset? No. Both the FDA and the EMA do not mandate any software or routine within. Quote from the Q&A: SAS (version 9.1, SAS Institute Inc, NC) was used in the previous computations. Results obtained by alternative, validated statistical programs are also acceptable except spreadsheets because outputs of spreadsheets are not suitable for secondary assessment. » You should try to get funded from EMA or FDA for these projects. These are really great projects, indeed. Gimme a break. I’m a lousy Rcoder. Agencies want to get results from validated software. They don’t pay for its development. — Diftor heh smusma 🖖 Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
yjlee168 ★★★ Kaohsiung, Taiwan, 20150421 23:41 (2417 d 19:29 ago) @ Helmut Posting: # 14730 Views: 21,916 

Dear Helmut, » ... The Q&A specifically mentions SAS’ Proc GLM , not Proc Mixed .glm() works in R. » In lme() it should be possible to specify all effects fixed. Actually the random effect(s) are optional. If is doesn’t work, lm() should do.lme() does not work without random effect. Am I doing anything wrong? > require(nlme) If I compared lme(), lm(), and glm(), I got > mod.lme<lme(log(Cmax) ~ seq + prd + drug, random=~drug 1subj, data=TotalData, method="REML") ###lme in bear Summary of mod.lme Linear mixedeffects model fit by REML Summary of mod.lm Call: Summary of glm Call: — All the best,  Yungjin Lee bear v2.9.0: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
Astea ★★ Russia, 20161104 00:13 (1855 d 17:56 ago) @ yjlee168 Posting: # 16771 Views: 18,400 

Dear all! Trying to understand lm and lme I found a simple and very clear instruction (Lme tutorial). Hope it would be useful for someone not familiar with such models. I try to evaluate Data set I from QA in R. The model is as follows: > mod.lme<lme(log(Cmax) ~ seq + prd + drug, random=~drugsubj/seq, data=QA) Then I try to interpret the data in order to estimate CI. I took 0.1465 as point estimation (exp(0,1465)=1,1578 and calculate CI as PE+ SE*t(0,1;219), where SE=0.04627. The result is 115,78: 107,26  124,97 (comparing with method B from QA might be 115,73: 107,17 124,97). Not ideal but seems to be similar. What is the reason for difference: rounding or different model, or am I do something wrong? Later I try to recalculate it using bear and get > Fixed effects: log(Cmax) ~ seq + prd + drug And finally if I try to use lme just from yjlee168's (post#14730), < mod.lme<lme(log(Cmax) ~ seq + prd + drug, random=~drug 1subj, data=TotalData, method="REML") I got another result:
Fixed effects: log(Cmax) ~ seq + prd + drug However it seemed to me strange, what model is in bear now? And I still wonder how does the model affect the Point Estimation? For QA data set 1 I've got manually T_{mean}=7,830303 and R_{mean}=7,676036, TR=0,1543 far from the estimated via lme PE. Trying to understand afermentioned Phoenix results I see noninteger df. How can they be achieved manually? — "Being in minority, even a minority of one, did not make you mad" 
mittyri ★★ Russia, 20161105 17:38 (1854 d 00:31 ago) @ Astea Posting: # 16777 Views: 18,218 

Dear Astea! so many questions I'll just put here my code as a solution. Yes, we cannot get SAS results, sorry... Please let me know if after investigation you still have some questions library(readxl) — Kind regards, Mittyri 
Astea ★★ Russia, 20161105 19:27 (1853 d 22:42 ago) @ mittyri Posting: # 16778 Views: 18,187 

mittyri!!! May I hug you?! This is what I actually wanted to get spending two nights in a desperate attempt to apply Alexandra Kuznetsova's lmerTest to BE! And the questions.. Still there are too many questions. Why didn't you provide the result of the code (that's the way of secret initiation?) First start with method B. My PE coincides fully with QA result (115.73), the CI calculated by your method slightly differs (107.26124.88). But! What is that function  confint? May be something wrong with it? Cause if I use PE+SE*t(0.1;219) where SE is 4.65e02 I get 1,07171,2797 totally the target!!! Please explain how to get Std.Error and Value with more digits? I tried to add "digits=8" but nothing happened... Why does the model have (1Subject) as random? I also tried (1Subject/Sequence) and the result (PE and CI) remains the same. Can we use it instead? The other question is what from the output we can use to estimate Rvariance in order to apply scABE? We are in front of the final step  from now we can make R tutorial to hold scABEL as did Helmut in WinNonLin, am I right? If it is so, Es ist fantastisch! — "Being in minority, even a minority of one, did not make you mad" 
mittyri ★★ Russia, 20161105 20:01 (1853 d 22:08 ago) @ Astea Posting: # 16779 Views: 18,165 

» mittyri!!! » May I hug you?! » Why didn't you provide the result of the code (that's the way of secret initiation?) I was so surprized by the results that was in hurry to share my code » Please explain how to get Std.Error and Value with more digits? Oh, you are right, please substitute confint string to FormulationT<summary(muddle.lmer)$coefficients["FormulationT",]
cat(round(exp(FormulationT[1]+FormulationT[2]*qt(0.05, FormulationT[3]))*100, 2), "", confint is not the right thing here. Here we go: lmer treats Sequence, Period and Formulation as factors, Random is /~1Subject/ (close to Method B), then [1] 0.1460882 PE is 115.73 SAS gives 115.73 what about LSMeans? PHX gives for Method B 7.67001367911898 and 7.81610190985527 and lmer gives Formulation Estimate what about CI? SAS gives 107.17124.97 107.17  124.97 » Why does the model have (1Subject) as random? I also tried (1Subject/Sequence) and the result (PE and CI) remains the same. Can we use it instead? Yes, we can. Subject in Sequence is SAS legacy. By the way as ElMaestro wrote, that's meaningless, because the reference grid (Gmatrix) doesn't change. It will be useful only if we are enumerating the subjects in each sequence from 1 (Seq ABAB: 1, 2, 3...; Seq BABA: 1, 2, 3...) » The other question is what from the output we can use to estimate Rvariance in order to apply scABE? No, we cannot. Sorry, I must confess I don't have time now to implement this, but I would do it as follows:  exclude all T data first and then  build a linear model with lm() using only reference data. If nobody does it during the next week I'll try to code it by myself. » We are in front of the final step  from now we can make R tutorial to hold scABEL as did Helmut in WinNonLin, am I right? If it is so, Es ist fantastisch! Ja, new feature to implement in BEAR, Yungjin will be happy — Kind regards, Mittyri 
Astea ★★ Russia, 20161106 11:50 (1853 d 06:19 ago) @ mittyri Posting: # 16781 Views: 18,180 

» If nobody does it during the next week I'll try to code it by myself. Challenge accepted! Here is my rough code: library(readxl) gives PE is 115.73 as in Phoenix (see Helmut's RSABE tutorial). The code needs to be tested via another data... — "Being in minority, even a minority of one, did not make you mad" 
mittyri ★★ Russia, 20161107 06:07 (1852 d 12:02 ago) @ Astea Posting: # 16782 Views: 18,289 

Hi Astea, » Here is my rough code: Nice work! Here's mine. I tried to prepare the dataframe for validation (Results) library(readxl) » as in Phoenix (see Helmut's RSABE tutorial). The code needs to be tested via another data... I suspect deviations from SAS, because some values are not equal in full precision... OK, that would be my next task — Kind regards, Mittyri 
yjlee168 ★★★ Kaohsiung, Taiwan, 20150420 21:34 (2418 d 21:36 ago) @ Helmut Posting: # 14720 Views: 21,693 

Dear Helmut, » I closed this thread. Since you have closed that thread, I like to thank you for this post first (very nice codes). But do you think that we should also present df, upper/lower CI (90% or 95%), and p values stuffs for lsmean? Or just keep it simple as you suggest? Thanks in advanced. — All the best,  Yungjin Lee bear v2.9.0: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear Download link (updated) > here 
Helmut ★★★ Vienna, Austria, 20150421 01:15 (2418 d 17:54 ago) @ yjlee168 Posting: # 14726 Views: 21,591 

Hi Yungjin, » I like to thank you for this post first Welcome. » (very nice codes). I would never (!) use such a code given the nice vectorization possible in R. It was only intended as an example to get the means without a model. » But do you think that we should also present df, upper/lower CI (90% or 95%), and p values stuffs for lsmean? IMHO, the df makes sense since it tells us whether we pulled the right stuff from the model. The CI is nice to know but not necessary. If the CIs overlap there is no significant difference between treatments. So what? No agency asks for it. I never understood the purpose of the t and pvalues. What is tested here? A difference to zero or what? » Or just keep it simple as you suggest? Please don’t think about using this awful code! — Diftor heh smusma 🖖 Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
Astea ★★ Russia, 20161102 23:43 (1856 d 18:27 ago) @ Helmut Posting: # 16767 Views: 18,292 

Dear all! While reading tonns of posts I stil can't understand how do SAS has got such a unique position? The method C for replicate studies (that is random effect with interactions PROC MIXED) in R (not only in bear but also in other packages) doesn't exist, does it? As I understand only method A is available in bear and only for partial replicate designs 2x2x4, 2x2x6.. But what about fully replicated 2x4x4? VStus (in that post) mentioned that » bear's lm.mod() was not confused by having more than 2 periods and 2 sequences Also a question appears what should we use as a reference variance in order to apply scaling of confidence intervals (fully replicated design allows us even to estimate intrasubject variability of reference drug for two pairs of R that is we can desintegrate 2x4x4 for two 2x2x4 and estimate R variance independently for both. How will choosing one of them affect the scaling and possibility to fail?) And the thing that is totally out of my mind: why should point estimation be affected by the way of analysis? As it was once shown by Helmut for balanced studies PE is just a means over periods and sequences. Why can't we do analogous in the case of inbalanced designs? — "Being in minority, even a minority of one, did not make you mad" 
VStus ★ Poland, 20161106 11:34 (1853 d 06:35 ago) @ Astea Posting: # 16780 Views: 18,069 

Dear Astea, » As I understand only method A is available in bear and only for partial replicate designs 2x2x4, 2x2x6.. But what about fully replicated 2x4x4? VStus (in that post) mentioned that I meant 3 and 4arm nonreplicate studies in my post, exploratory pilot comparing 2 or 3 pilot formulations versus the Reference. lm() or PROC GLM is enough in this case. Regards, VStus 
StatR ☆ Turkey, 20170203 13:53 (1764 d 04:16 ago) @ VStus Posting: # 17007 Views: 17,009 

Hi dear all, First of all, I am sorry if my question is too stupid. I could not find the answer no matter how much I investigate. I am new in the area of bioequivalence, and recently I am interested in individual bioequivalence. Also, I am not so familiar with the lme4. I am trying to extract the variance components (sigmaBT^2, sigmaBR^2, sigmaWT^2, sigmaWR^2 and rho) of the model (Restricted Maximum Likelihood Estimation). My data is as follows: head(mydata) I am trying to do it by lmer function as follows: lmer(Yijk~Formula+Period+Sequence+(1Sequence)) I know, the estimated value of sigmaWR^2 = 1.269e01 and the mean difference FtFr = 0.178645. But how can I get the estimated values of sigmaBT^2, sigmaBR^2, sigmaWT^2 and rho ? I appreciate any help. 
VStus ★ Poland, 20170203 15:47 (1764 d 02:23 ago) (edited by VStus on 20170203 16:34) @ StatR Posting: # 17009 Views: 16,791 

Dear StatR, Sorry for confusion, deleted my incorrect reponse. summary(lmerdata) seems to be already provided.What is Yijk ? Yijk seems to be a PK paramenter of interest... Residual you've extracted seems to be pooled both formulations. Regards, VStus PS: Helmut, thank you for options(digits=12) ! 
StatR ☆ Turkey, 20170203 17:12 (1764 d 00:57 ago) @ VStus Posting: # 17011 Views: 16,703 

» Yijk seems to be a PK paramenter of interest... Residual you've extracted seems to be pooled both formulations. Dear VStur, Thanks for your reply. You are right, Yijk is PK characteristic of interest. The design is 2 sequence 2 treatment and 4 period design(TRTR/RTRT). More clearly, I want to extract between subject variances (sigmaBT^2, sigmaBR^2), within subject variances (sigmaWT^2, sigmaWR^2) under the test and reference formulation, also subjectbyformulation interaction (sigmaD^2 = sigmaBT^2 + sigmaBR^2  2*rho*sigmaBT*sigmaBR). Is it possible to obtain all of these estimates from the summary(lmer_model), or do I need to do an extra step to obtain those values? Thanks in advance Edit: Full quote removed. Please delete everything from the text of the original poster which is not necessary in understanding your answer; see also this post #5! [Helmut] 
d_labes ★★★ Berlin, Germany, 20170207 11:16 (1760 d 06:53 ago) @ StatR Posting: # 17029 Views: 16,597 

Dear StatR, » ... More clearly, I want to extract between subject variances (sigmaBT^2, sigmaBR^2), within subject variances (sigmaWT^2, sigmaWR^2) under the test and reference formulation, also subjectbyformulation interaction (sigmaD^2 = sigmaBT^2 + sigmaBR^2  2*rho*sigmaBT*sigmaBR). What is needed for that goal is to formulate a model which have these variancecovariance parameters. Your model is one with assuming sigmaWT = sigmaWR and sigmaBT=sigmaBR, rho=1. Thus you have only two variancecovariance parameters in your output above: Random effects: Unfortunately there is currently no way to formulate a model similar to the FDA code for replicate designs within package lme4 / lmer() . See f.i. this post for more details on the why.In the older package nlme / lme() you may formulate the model similar to the FDA code. See that post.But nevertheless you can't get the 'correct' CI from that R solution since the degrees of freedom for the contrast TR are different to that of SAS Proc MIXED. Hope this helps. — Regards, Detlew 
StatR ☆ Turkey, 20170207 11:36 (1760 d 06:33 ago) @ d_labes Posting: # 17030 Views: 16,518 

» What is needed for that goal is to formulate a model which have these variancecovariance parameters. » Your model is one with assuming sigmaWT = sigmaWR and sigmaBT=sigmaBR, rho=1. » Thus you have only two variancecovariance parameters Dear d_labes, Thank you very much for the explanation. At least, now, I know what should I do. Best. 
StatR ☆ Turkey, 20170208 08:41 (1759 d 09:28 ago) @ d_labes Posting: # 17033 Views: 16,482 

Hi again d_labes, I tried to find the correct model by the post you referred. » In the older package nlme / lme() you may formulate the model similar to the FDA code. See that post.Unfortunately, I am confused after reading tons of post, and I could not figure out what is "trt", "tmt" ... Simply, my data is as follows Subject Formula Period Sequence Yijk Where Yijk is the PK characteristic. The first sequence represents "TRTR" and the second sequence: "RTRT". In the formula 1 = T and 0 = R. In the referred post you say that model2 < lme(y~ drug + prd + seq, In your model y = Yijk, drug=Formula, prd=Period, seq=Sequence, but what is tmt? Also, is the my data format correct for the lme? Or, Should I change Formula as 1 = T, 2 = R and Sequence as 1 = TRTR, 2 = RTRT? Could you please make me clear. Best. 
d_labes ★★★ Berlin, Germany, 20170208 10:13 (1759 d 07:57 ago) @ StatR Posting: # 17034 Views: 16,582 

Dear StatR, sorry for confusing you. I was a little bit sloppy in the post and mixed the factor names in Bear with mine. To clarify my naming convention and coding: tmt is treatment coded as "T" or "R" period is period no sequence is sequence coded as "TRTR" or "RTRT" (or whatever you have) subject is subject no. y is logtransformed PK metric (Cmax or AUC or whatever you like to evaluate) Don't forget to make them as.factor()! Then call muddle < lme(y ~ tmt + period + sequence, If you prefer your numeric coding of tmt = Formula and Sequence: Don't forget to make them as.factor()! Then call muddle < lme(y ~ Formula + Period + Sequence, Using the dataset II from the EMA Q&A you will get the following output with my convention: Linear mixedeffects model fit by REML From that it is read (note that lme() works with SD instead of variances): s2wT = 0.09296929^2 Hope this helps. Homework: Identify the variancecovariance terms with your coding . — Regards, Detlew 
StatR ☆ Turkey, 20170208 10:19 (1759 d 07:50 ago) @ d_labes Posting: # 17035 Views: 16,475 

Dear d_labes, The only thing that I say: Wow! I do not know how I thanks for this reply. Thank you! » sorry for confusing you. I was a little bit sloppy in the post and mixed the factor names in Bear with mine. » Homework: Identify the variancecovariance terms with your coding . 
Helmut ★★★ Vienna, Austria, 20170208 10:33 (1759 d 07:36 ago) @ StatR Posting: # 17036 Views: 16,559 

Hi StatR, » Subject Formula Period Sequence Yijk » 1 1 1 1 4.933014 » 1 0 2 1 5.006829 » 1 1 3 1 4.854162 » 1 0 4 1 5.081755 » ... » » Where Yijk is the PK characteristic. The first sequence represents "TRTR" and the second sequence: "RTRT". In the formula 1 = T and 0 = R. Only a general remark about data structures in R. I suggest to code categorial (classification) variables as characters instead of numeric (Treatment, Sequence). It not only helps others to understand what you have meant without requiring an explanation (1=TRTR, 2=RTRT and 1=T, 0=R) but also models (e.g., lm() , lme() , lmer() ) will automatically treat them as factors. Avoid Formula and use Treatment or Formulation instead. Personally I prefer Treatment over Formulation for generality (imagine a foodeffect study where we have the same formulation).Convert numeric variables (except the PKresponse) to factors as well:
» Unfortunately, I am confused after reading tons of post, and I could not figure out what is "trt", "tmt" ... Always the treatment is meant. To add confusion in package bear Yungjin uses drug . — Diftor heh smusma 🖖 Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
StatR ☆ Turkey, 20170208 10:49 (1759 d 07:20 ago) @ Helmut Posting: # 17037 Views: 16,510 

Hi Helmut, Many thans for your suggestions. I think you are right, otherwise it's confusing. I got the message Best » Only a general remark about data structures in R. I suggest to code categorial (classification) variables as characters instead of numeric (Treatment, Sequence). It not only helps others to understand what you have meant without requiring an explanation (1=TRTR, 2=RTRT and 1=T, 0=R) but also models (e.g., lm() , lme() , lmer() ) will automatically treat them as factors. Avoid Formula and use Treatment or Formulation instead. Personally I prefer Treatment over Formulation for generality... 