Examp2.5.3.1 from Duchateau, L. and Janssen, P. and Rowlands, G. J. (1998).<em>Linear Mixed Models. An Introduction with applications in Veterinary Research</em>. International Livestock Research Institute.

Examp2.5.3.1 is used for inspecting probability distribution and to define a plausible process through linear models and generalized linear models.

References

Duchateau, L. and Janssen, P. and Rowlands, G. J. (1998).Linear Mixed Models. An Introduction with applications in Veterinary Research. International Livestock Research Institute.

Examples

#-------------------------------------------------------------
## Example 2.5.3.1 p-70
#-------------------------------------------------------------
 # PROC GLM DATA=ex125;
 # CLASS drug dose region;
 # MODEL pcv=region drug region*drug dose drug*dose;
 # RANDOM region drug*region;
 # RUN;


 # PROC MIXED DATA=ex125;
 # CLASS drug dose region;
 # MODEL pcv=drug dose drug*dose / ddfm=satterth;
 # RANDOM region drug*region;
 # ESTIMATE 'drug dif' drug -1 1 drug*dose -0.5 -0.5 0.5 0.5;
 # ESTIMATE 'Samorin mean' INTERCEPT 1 drug 0 1 dose 0.5 0.5
 #                             drug*dose 0 0 0.5 0.5;
 # ESTIMATE 'Samorin HvsL' dose 1 -1 drug*dose 0 0 1 -1;
 # ESTIMATE 'Samorin high' INTERCEPT 1 drug 0 1 dose 1 0
 #                             drug*dose 0 0 1 0;
 # RUN;

library(lme4)
str(ex125)
#> 'data.frame':	24 obs. of  4 variables:
#>  $ Region: int  1 1 1 1 2 2 2 2 3 3 ...
#>  $ Drug  : Factor w/ 2 levels "BERENIL","samorin": 1 1 2 2 1 1 2 2 1 1 ...
#>  $ dose  : Factor w/ 2 levels "h","l": 1 2 1 2 1 2 1 2 1 2 ...
#>  $ Pcv   : num  22.6 21.8 19.1 16.4 29 28.8 25.3 18.2 24 23.7 ...
ex125$Region1 <- factor(ex125$Region)
 fm2.11 <-
  aov(
      formula     = Pcv ~ Region1 + Drug + Error(Drug:Region1) + dose + dose:Drug
    , data        = ex125
    , projections = FALSE
    , qr          = TRUE
    , contrasts   = NULL
  #  , ...
    )
#> Warning: Error() model is singular
 summary(fm2.11)
#> 
#> Error: Drug:Region1
#>           Df Sum Sq Mean Sq F value   Pr(>F)    
#> Region1    5 117.37   23.47   8.178 0.018857 *  
#> Drug       1 185.37  185.37  64.580 0.000483 ***
#> Residuals  5  14.35    2.87                     
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Error: Within
#>           Df Sum Sq Mean Sq F value   Pr(>F)    
#> dose       1  45.65   45.65  21.775 0.000886 ***
#> Drug:dose  1  15.20   15.20   7.251 0.022594 *  
#> Residuals 10  20.96    2.10                     
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 fm2.12 <-
  lmerTest::lmer(
         formula    = Pcv ~ dose*Drug + (1|Region/Drug)
       , data       = ex125
       , REML       = TRUE
       , control    = lmerControl()
       , start      = NULL
       , verbose    = 0L
    #  , subset
    #  , weights
    #  , na.action
    #  , offset
       , contrasts  = list(dose = "contr.SAS", Drug = "contr.SAS")
       , devFunOnly = FALSE
    #  , ...
       )
 summary(fm2.12)
#> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
#> lmerModLmerTest]
#> Formula: Pcv ~ dose * Drug + (1 | Region/Drug)
#>    Data: ex125
#> 
#> REML criterion at convergence: 92.4
#> 
#> Scaled residuals: 
#>      Min       1Q   Median       3Q      Max 
#> -1.71077 -0.51628  0.08414  0.70276  1.16129 
#> 
#> Random effects:
#>  Groups      Name        Variance Std.Dev.
#>  Drug:Region (Intercept) 0.387    0.6221  
#>  Region      (Intercept) 5.151    2.2695  
#>  Residual                2.096    1.4479  
#> Number of obs: 24, groups:  Drug:Region, 12; Region, 6
#> 
#> Fixed effects:
#>                   Estimate Std. Error      df t value Pr(>|t|)    
#> (Intercept)        17.1333     1.1280  8.2080  15.189 2.69e-07 ***
#> doseh               4.3500     0.8359 10.0000   5.204 0.000399 ***
#> DrugBERENIL         7.1500     0.9098 11.8185   7.859 4.96e-06 ***
#> doseh:DrugBERENIL  -3.1833     1.1822 10.0000  -2.693 0.022594 *  
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Correlation of Fixed Effects:
#>             (Intr) doseh  DBEREN
#> doseh       -0.371              
#> DrugBERENIL -0.403  0.459       
#> ds:DBERENIL  0.262 -0.707 -0.650
 anova(object = fm2.12, ddf = "Satterthwaite")
#> Type III Analysis of Variance Table with Satterthwaite's method
#>           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#> dose       45.65   45.65     1    10 21.7755 0.0008856 ***
#> Drug      135.39  135.39     1     5 64.5796 0.0004826 ***
#> dose:Drug  15.20   15.20     1    10  7.2507 0.0225945 *  
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(multcomp)
#> Loading required package: mvtnorm
#> Loading required package: survival
#> Loading required package: TH.data
#> Loading required package: MASS
#> 
#> Attaching package: ‘TH.data’
#> The following object is masked from ‘package:MASS’:
#> 
#>     geyser
Contrasts1 <-
          matrix(c(
                  1, 0.5, 0, 0
                , 0, 0, -1, -0.5
                , 1, 1, 0, 0
                , 0, 1, 0, 0
                )
               , ncol = 4
               , byrow = TRUE
               , dimnames = list(
                  c("C1", "C2", "C3", "C4")
                , rownames(summary(fm2.12)$coef)
               )
              )

Contrasts1
#>    (Intercept) doseh DrugBERENIL doseh:DrugBERENIL
#> C1           1   0.5           0               0.0
#> C2           0   0.0          -1              -0.5
#> C3           1   1.0           0               0.0
#> C4           0   1.0           0               0.0
summary(glht(fm2.12, linfct=Contrasts1))
#> 
#> 	 Simultaneous Tests for General Linear Hypotheses
#> 
#> Fit: lmerTest::lmer(formula = Pcv ~ dose * Drug + (1 | Region/Drug), 
#>     data = ex125, REML = TRUE, control = lmerControl(), start = NULL, 
#>     verbose = 0L, contrasts = list(dose = "contr.SAS", Drug = "contr.SAS"), 
#>     devFunOnly = FALSE)
#> 
#> Linear Hypotheses:
#>         Estimate Std. Error z value Pr(>|z|)    
#> C1 == 0  19.3083     1.0477  18.429  < 1e-07 ***
#> C2 == 0  -5.5583     0.6917  -8.036  < 1e-07 ***
#> C3 == 0  21.4833     1.1280  19.046  < 1e-07 ***
#> C4 == 0   4.3500     0.8359   5.204 7.21e-07 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> (Adjusted p values reported -- single-step method)
#>

Examp2.5.3.1 from Duchateau, L. and Janssen, P. and Rowlands, G. J. (1998).Linear Mixed Models. An Introduction with applications in Veterinary Research. International Livestock Research Institute.

References

See also

Examples

Contents

Author