WHM
☆    

Korea,
2015-02-16 08:54
(3328 d 06:51 ago)

Posting: # 14425
Views: 5,631
 

 check code for Type III in R [🇷 for BE/BA]

Dear R users

I've used SAS proc glm or WNL for BE.
Recently, I've tried to get type III anova table (because of my manager who loves R). I've noticed that it is difficult to manage type III ANOVA Table.
I read prior posts and compared SAS proc glm, proc mixed and WNL, so I've known some problems and solutions. I've tried drop1() mentioned prior posts but it did now show all SAS showed.
firstly, I(actually my manager) want to take basic code using lm() for unbalanced data which is 2x2x2 crossover design.
I used two lm() models. one is treatment(TRT) behind period(PRD), the other is period(PRD) behind treatment(TRT).
so here is my code for type III anova table and sas results(WNL is also same result).
Can I use this code for type III anova table in 2x2x2 crossover design ?

----R code--------------------------------------------
lm.test1<-lm(log(Cmax)~TRT+PRD+SEQ+SID:SEQ, data=df)
lm.test2<-lm(log(Cmax)~PRD+TRT+SEQ+SID:SEQ, data=df)
anova.F<-anova(lm.test1)
# trt result to lm.test1 to make type III
anova.F[1,]<-anova(lm.test2)[2,]
#Tests of Hypotheses Using the Type III MS for subject(sequence) as an Error Term
anova.F[3,4]<-anova.F[3,3]/anova.F[4,3] 
anova.F[3,5]<-pf(anova.F[3,4],anova.F[3,1],anova.F[4,1],lower.tail =F)


----R result--------------------------------------------
Analysis of Variance Table

Response: log(Cmax)
          Df  Sum Sq Mean Sq F value    Pr(>F)   
TRT        1  0.1181 0.11805  1.3950    0.2436   
PRD        1  0.0765 0.07649  0.9039    0.3467   
SEQ        1  0.5889 0.58890  1.6448    0.2061   
SEQ:SID   46 16.4697 0.35804  4.2308 1.458e-06 ***
Residuals 46  3.8928 0.08463                     


----SAS result(proc glm)--------------------------------------------

Source DF Type III SS Mean Square F Value Pr > F
SEQ 1 0.58889907 0.58889907 6.96 0.0113
SID(SEQ) 46 16.46968620 0.35803666 4.23 <.0001
PRD 1 0.07648993 0.07648993 0.90 0.3467
TRT 1 0.11805116 0.11805116 1.39 0.2436

Tests of Hypotheses Using the Type III MS for subject(sequence) as an Error Term
Source DF Type III SS Mean Square F Value Pr > F
SEQ 1 0.58889907 0.58889907 1.64 0.2061


Best regards
WHM
ElMaestro
★★★

Denmark,
2015-02-16 13:17
(3328 d 02:27 ago)

@ WHM
Posting: # 14427
Views: 4,592
 

 check code for Type III in R

Hi WHM,

❝ I've tried drop1() mentioned prior posts but it did now show all SAS showed.


Your code -at a quick glance- looks ok, but this is said without having tried it on data.
You should not necessarily just think that whatever comes out of SAS is what you want to reproduce.

Type III means you will take out each of the various factors one by one wile keeping the others in the model. Your differences in residuals vs the full model residual correspond to the Type III SS for the individual factors. By that definition sequence will not have any appreciable SS (apart from the numerical precision) and precisely that is what you will find in R with drop1. Think about it this way: Sequence is a between-subject factor; you cannot force in an extra column for Sequence RT or TR in the design matrix since any such column can be constructed from the Subject columns and the intercepts with a little addition and subtraction.

SAS will (with the random bogus statement) goa step further and get you a type III SS for the Sequence. This is not plain type III, but perhaps called type III with a useful twist. I will leave the interpretation to you but note that R's output with drop1 is not 'wrong'.

Pass or fail!
ElMaestro
WHM
☆    

Korea,
2015-02-24 02:56
(3320 d 12:49 ago)

@ ElMaestro
Posting: # 14495
Views: 4,475
 

 check code for Type III in R

Hi ElMaestro,

❝ ❝ I've tried drop1() mentioned prior posts but it did now show all SAS showed.

(I wrote wrong word. I meant it did not show all SAS showed) anyway,,

First of all, Thank you for your quick reply and sorry for late my reply.

I run drop1() you recommended. the results are below.
----R code--------------------------------------------
>lm.test1<-lm(log(Cmax)~TRT+PRD+SEQ+SID:SEQ, data=df)
>drop1(lm.test1, ~.,  test=c('F'))
----R result-------------------------------------------
Single term deletions

Model:
log(Cmax) ~ TRT + PRD + SEQ + SID:SEQ
        Df Sum of Sq     RSS     AIC F value    Pr(>F)   
<none>                3.8928 -207.70                     
TRT      1    0.1181  4.0109 -206.83  1.3950    0.2436   
PRD      1    0.0765  3.9693 -207.83  0.9039    0.3467   
SEQ      0    0.0000  3.8928 -207.70                     
SEQ:SID 46   16.4697 20.3625 -140.86  4.2308 1.458e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------------------

the results are same with SAS result except for SEQ and residuals.
if i want to get the SS of SEQ and residuals, should i use other function? or can get from drop1()?

You mentioned below. but I'm sorry I did not exactly understand about sequence.
Could you explain that more in detail?

❝ By that definition sequence will not have any appreciable SS (apart from the numerical precision) and precisely that is what you will find in R with drop1. Think about it this way: Sequence is a between-subject factor; you cannot force in an extra column for Sequence RT or TR in the design matrix since any such column can be constructed from the Subject columns and the intercepts with a little addition and subtraction.



Best Regards
WHM
ElMaestro
★★★

Denmark,
2015-02-24 10:59
(3320 d 04:46 ago)

@ WHM
Posting: # 14498
Views: 4,479
 

 check code for Type III in R

Hi WHM,

❝ if i want to get the SS of SEQ and residuals, should i use other function? or can get from drop1()?


Specifically for
a: lm.test1<-lm(log(Cmax)~TRT+PRD+SEQ+SID:SEQ, data=df)

you will see a zero value for Sequence.
If you want to mimic SAS' output, then I think you could try
b: lm.test1<-lm(log(Cmax)~TRT+PRD+SEQ, data=df)
and then you grab the F-value and P-value you get for Sequence in "b" and insert it in the anova table that you got with "a".

I hope this revolves it.

Pass or fail!
ElMaestro
UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,638 registered users;
69 visitors (0 registered, 69 guests [including 5 identified bots]).
Forum time: 15:45 CET (Europe/Vienna)

Nothing shows a lack of mathematical education more
than an overly precise calculation.    Carl Friedrich Gauß

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