############################Split-plot experiment##########################
library(lmerTest) #this is the package suggested
## Loading required package: Matrix
## Loading required package: lme4
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(nlme) #for random effects, degrees of freedom has some problem
##
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
##
## lmList
require(ggplot2)
## Loading required package: ggplot2
require(reshape2) #to change long format to wide format or change from wide format to long format
## Loading required package: reshape2
library(heplots) #test equivalence of covariance matrix
## Loading required package: car
library(multcomp) #for multiple comparison
## 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
library(xtable)
library(lsmeans)
## Loading required package: estimability
##
## Attaching package: 'lsmeans'
## The following object is masked from 'package:lmerTest':
##
## lsmeans
#hrate is in long format
hrate<-read.table("C:/jenn/teaching/stat579/data/heartrate.txt",header=TRUE)
nrow(hrate) #120
## [1] 120
hrate
## drug subject time rate
## 1 a 11 time1 81
## 2 a 11 time2 81
## 3 a 11 time3 82
## 4 a 11 time4 82
## 5 a 12 time1 82
## 6 a 12 time2 83
## 7 a 12 time3 80
## 8 a 12 time4 81
## 9 a 13 time1 81
## 10 a 13 time2 77
## 11 a 13 time3 80
## 12 a 13 time4 80
## 13 a 14 time1 84
## 14 a 14 time2 86
## 15 a 14 time3 85
## 16 a 14 time4 85
## 17 a 15 time1 88
## 18 a 15 time2 90
## 19 a 15 time3 88
## 20 a 15 time4 86
## 21 a 16 time1 83
## 22 a 16 time2 82
## 23 a 16 time3 86
## 24 a 16 time4 85
## 25 a 17 time1 85
## 26 a 17 time2 83
## 27 a 17 time3 87
## 28 a 17 time4 86
## 29 a 18 time1 81
## 30 a 18 time2 85
## 31 a 18 time3 86
## 32 a 18 time4 85
## 33 a 19 time1 87
## 34 a 19 time2 89
## 35 a 19 time3 87
## 36 a 19 time4 82
## 37 a 20 time1 77
## 38 a 20 time2 75
## 39 a 20 time3 73
## 40 a 20 time4 77
## 41 b 21 time1 76
## 42 b 21 time2 83
## 43 b 21 time3 85
## 44 b 21 time4 79
## 45 b 22 time1 75
## 46 b 22 time2 81
## 47 b 22 time3 85
## 48 b 22 time4 73
## 49 b 23 time1 75
## 50 b 23 time2 82
## 51 b 23 time3 80
## 52 b 23 time4 77
## 53 b 24 time1 68
## 54 b 24 time2 73
## 55 b 24 time3 72
## 56 b 24 time4 69
## 57 b 25 time1 78
## 58 b 25 time2 87
## 59 b 25 time3 86
## 60 b 25 time4 77
## 61 b 26 time1 81
## 62 b 26 time2 85
## 63 b 26 time3 81
## 64 b 26 time4 74
## 65 b 27 time1 67
## 66 b 27 time2 73
## 67 b 27 time3 75
## 68 b 27 time4 66
## 69 b 28 time1 68
## 70 b 28 time2 73
## 71 b 28 time3 73
## 72 b 28 time4 66
## 73 b 29 time1 68
## 74 b 29 time2 75
## 75 b 29 time3 79
## 76 b 29 time4 69
## 77 b 30 time1 73
## 78 b 30 time2 78
## 79 b 30 time3 80
## 80 b 30 time4 70
## 81 p 1 time1 80
## 82 p 1 time2 77
## 83 p 1 time3 73
## 84 p 1 time4 69
## 85 p 2 time1 64
## 86 p 2 time2 66
## 87 p 2 time3 68
## 88 p 2 time4 71
## 89 p 3 time1 75
## 90 p 3 time2 73
## 91 p 3 time3 73
## 92 p 3 time4 69
## 93 p 4 time1 72
## 94 p 4 time2 70
## 95 p 4 time3 74
## 96 p 4 time4 73
## 97 p 5 time1 74
## 98 p 5 time2 74
## 99 p 5 time3 71
## 100 p 5 time4 67
## 101 p 6 time1 71
## 102 p 6 time2 71
## 103 p 6 time3 72
## 104 p 6 time4 70
## 105 p 7 time1 76
## 106 p 7 time2 78
## 107 p 7 time3 74
## 108 p 7 time4 71
## 109 p 8 time1 73
## 110 p 8 time2 68
## 111 p 8 time3 64
## 112 p 8 time4 64
## 113 p 9 time1 76
## 114 p 9 time2 73
## 115 p 9 time3 74
## 116 p 9 time4 76
## 117 p 10 time1 77
## 118 p 10 time2 78
## 119 p 10 time3 77
## 120 p 10 time4 73
hrate$time<-factor(hrate$time)
hrate$drug<-factor(hrate$drug)
hrate$subject<-factor(hrate$subject)
boxplot(rate~drug, data=hrate, main="boxplot of rate vs drug")
#use ggplot for rate vs drug for each time period
p <- ggplot(hrate, aes(factor(time), rate)) + geom_boxplot(aes(fill = drug))
p
## d2 is a data set in wide format, with 30 variables and 6 columns, drug,
#subject, time1-time4, you can change d2 to the long format, such as the data
#hrate with 120 obsn and 4 columns drug, subject, time, and rate
d2<-read.table("C:/jenn/teaching/stat579/data/splitd2.txt",header=TRUE)
d2.long <- melt(d2,
# id.vars: ID variables
# all variables to keep but not split apart on
id.vars=c("subject","drug"),
# measure.vars: The source columns
# (if unspecified then all other variables are measure.vars)
measure.vars = c("time1", "time2", "time3", "time4"),
# variable.name: Name of the destination column identifying each
# original column that the measurement came from
variable.name = "time",
# value.name: column name for values in table
value.name = "rate"
)
nrow(d2)
## [1] 30
d2
## drug subject time1 time2 time3 time4
## 1 p 11 80 77 73 69
## 2 p 12 64 66 68 71
## 3 p 13 75 73 73 69
## 4 p 14 72 70 74 73
## 5 p 15 74 74 71 67
## 6 p 16 71 71 72 70
## 7 p 17 76 78 74 71
## 8 p 18 73 68 64 64
## 9 p 19 76 73 74 76
## 10 p 20 77 78 77 73
## 11 a 21 81 81 82 82
## 12 a 22 82 83 80 81
## 13 a 23 81 77 80 80
## 14 a 24 84 86 85 85
## 15 a 25 88 90 88 86
## 16 a 26 83 82 86 85
## 17 a 27 85 83 87 86
## 18 a 28 81 85 86 85
## 19 a 29 87 89 87 82
## 20 a 30 77 75 73 77
## 21 b 1 76 83 85 79
## 22 b 2 75 81 85 73
## 23 b 3 75 82 80 77
## 24 b 4 68 73 72 69
## 25 b 5 78 87 86 77
## 26 b 6 81 85 81 74
## 27 b 7 67 73 75 66
## 28 b 8 68 73 73 66
## 29 b 9 68 75 79 69
## 30 b 10 73 78 80 70
nrow(d2.long)
## [1] 120
d2.long[1:20,]
## subject drug time rate
## 1 11 p time1 80
## 2 12 p time1 64
## 3 13 p time1 75
## 4 14 p time1 72
## 5 15 p time1 74
## 6 16 p time1 71
## 7 17 p time1 76
## 8 18 p time1 73
## 9 19 p time1 76
## 10 20 p time1 77
## 11 21 a time1 81
## 12 22 a time1 82
## 13 23 a time1 81
## 14 24 a time1 84
## 15 25 a time1 88
## 16 26 a time1 83
## 17 27 a time1 85
## 18 28 a time1 81
## 19 29 a time1 87
## 20 30 a time1 77
d2p<-d2[1:10,3:6]
d2a<-d2[11:20,3:6]
d2b<-d2[21:30,3:6]
cov(d2p)
## time1 time2 time3 time4
## time1 18.6222222 15.066667 8.222222 0.7333333
## time2 15.0666667 17.066667 11.111111 3.0666667
## time3 8.2222222 11.111111 13.333333 8.8888889
## time4 0.7333333 3.066667 8.888889 11.3444444
cov(d2a)
## time1 time2 time3 time4
## time1 10.544444 13.56667 12.93333 6.766667
## time2 13.566667 22.54444 18.62222 10.122222
## time3 12.933333 18.62222 21.82222 12.822222
## time4 6.766667 10.12222 12.82222 8.988889
cov(d2b)
## time1 time2 time3 time4
## time1 24.10000 25.11111 19.17778 18.88889
## time2 25.11111 28.22222 23.11111 22.33333
## time3 19.17778 23.11111 24.93333 19.00000
## time4 18.88889 22.33333 19.00000 22.00000
##test if the three covariance matrices are equivalent across the groups
tres <- boxM(d2[, 3:6], d2[, "drug"])
tres
##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: d2[, 3:6]
## Chi-Sq (approx.) = 24.408, df = 20, p-value = 0.225
summary(tres)
## Summary for Box's M-test of Equality of Covariance Matrices
##
## Chi-Sq: 24.40793
## df: 20
## p-value: 0.225
##
## log of Covariance determinants:
## a b p pooled
## 5.821413 6.824833 7.334394 7.807921
##
## Eigenvalues:
## a b p pooled
## 1 56.2796579 89.0611189 41.3959082 61.387858
## 2 5.1989142 5.7388204 15.3936986 8.675867
## 3 1.7697531 4.0060892 2.6799066 2.785441
## 4 0.6516748 0.4495271 0.8971533 1.658242
##
## Statistics based on eigenvalues:
## a b p pooled
## product 337.4486511 920.4226490 1532.0986435 2460.0118833
## sum 63.9000000 99.2555556 60.3666667 74.5074074
## precision 0.4329614 0.3759879 0.6341546 0.9144028
## max 56.2796579 89.0611189 41.3959082 61.3878580
##checking interactions #seems interaction is significant
interaction.plot(hrate$drug,hrate$time,hrate$rate)
##fit full model with two treatments
##note that there are 30 subjects with 1-30, identifiable betweent the groups
##adding nested structure or not doesn't make difference
myfit<-lmer(rate~drug*time+(1|subject),data=hrate)
summary(myfit)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: rate ~ drug * time + (1 | subject)
## Data: hrate
##
## REML criterion at convergence: 571.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39109 -0.52260 -0.05669 0.50612 2.41993
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 13.864 3.723
## Residual 4.763 2.182
## Number of obs: 120, groups: subject, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.290e+01 1.365e+00 4.057e+01 60.741 < 2e-16 ***
## drugb -1.000e+01 1.930e+00 4.057e+01 -5.181 6.41e-06 ***
## drugp -9.100e+00 1.930e+00 4.057e+01 -4.715 2.85e-05 ***
## timetime2 2.000e-01 9.760e-01 8.100e+01 0.205 0.8382
## timetime3 5.000e-01 9.760e-01 8.100e+01 0.512 0.6099
## timetime4 1.057e-14 9.760e-01 8.100e+01 0.000 1.0000
## drugb:timetime2 5.900e+00 1.380e+00 8.100e+01 4.274 5.20e-05 ***
## drugp:timetime2 -1.200e+00 1.380e+00 8.100e+01 -0.869 0.3872
## drugb:timetime3 6.200e+00 1.380e+00 8.100e+01 4.492 2.32e-05 ***
## drugp:timetime3 -2.300e+00 1.380e+00 8.100e+01 -1.666 0.0995 .
## drugb:timetime4 -9.000e-01 1.380e+00 8.100e+01 -0.652 0.5162
## drugp:timetime4 -3.500e+00 1.380e+00 8.100e+01 -2.536 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drugb drugp timtm2 timtm3 timtm4 drgb:2 drgp:2 drgb:3
## drugb -0.707
## drugp -0.707 0.500
## timetime2 -0.358 0.253 0.253
## timetime3 -0.358 0.253 0.253 0.500
## timetime4 -0.358 0.253 0.253 0.500 0.500
## drugb:tmtm2 0.253 -0.358 -0.179 -0.707 -0.354 -0.354
## drugp:tmtm2 0.253 -0.179 -0.358 -0.707 -0.354 -0.354 0.500
## drugb:tmtm3 0.253 -0.358 -0.179 -0.354 -0.707 -0.354 0.500 0.250
## drugp:tmtm3 0.253 -0.179 -0.358 -0.354 -0.707 -0.354 0.250 0.500 0.500
## drugb:tmtm4 0.253 -0.358 -0.179 -0.354 -0.354 -0.707 0.500 0.250 0.500
## drugp:tmtm4 0.253 -0.179 -0.358 -0.354 -0.354 -0.707 0.250 0.500 0.250
## drgp:3 drgb:4
## drugb
## drugp
## timetime2
## timetime3
## timetime4
## drugb:tmtm2
## drugp:tmtm2
## drugb:tmtm3
## drugp:tmtm3
## drugb:tmtm4 0.250
## drugp:tmtm4 0.500 0.500
print(xtable(anova(myfit))) #this is for the latex version to enter the table
## % latex table generated in R 3.2.5 by xtable 1.8-2 package
## % Tue Mar 27 00:04:06 2018
## \begin{table}[ht]
## \centering
## \begin{tabular}{lrrrrrr}
## \hline
## & Sum Sq & Mean Sq & NumDF & DenDF & F.value & Pr($>$F) \\
## \hline
## drug & 192.89 & 96.44 & 2.00 & 27.00 & 20.25 & 0.0000 \\
## time & 222.29 & 74.10 & 3.00 & 81.00 & 15.56 & 0.0000 \\
## drug:time & 320.13 & 53.36 & 6.00 & 81.00 & 11.20 & 0.0000 \\
## \hline
## \end{tabular}
## \end{table}
anova(myfit)
## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## drug 192.89 96.443 2 27 20.247 4.249e-06 ***
## time 222.29 74.097 3 81 15.556 4.444e-08 ***
## drug:time 320.13 53.356 6 81 11.201 4.539e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
myfit2<-lmer(rate~drug*time+(1|subject/drug),data=hrate)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
summary(myfit2)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: rate ~ drug * time + (1 | subject/drug)
## Data: hrate
##
## REML criterion at convergence: 571.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39109 -0.52260 -0.05669 0.50612 2.41993
##
## Random effects:
## Groups Name Variance Std.Dev.
## drug:subject (Intercept) 6.773 2.602
## subject (Intercept) 7.091 2.663
## Residual 4.763 2.182
## Number of obs: 120, groups: drug:subject, 30; subject, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.290e+01 1.365e+00 4.057e+01 60.741 < 2e-16 ***
## drugb -1.000e+01 1.930e+00 4.057e+01 -5.181 6.41e-06 ***
## drugp -9.100e+00 1.930e+00 4.057e+01 -4.715 2.85e-05 ***
## timetime2 2.000e-01 9.760e-01 8.100e+01 0.205 0.8382
## timetime3 5.000e-01 9.760e-01 8.100e+01 0.512 0.6099
## timetime4 2.968e-14 9.760e-01 8.100e+01 0.000 1.0000
## drugb:timetime2 5.900e+00 1.380e+00 8.100e+01 4.274 5.20e-05 ***
## drugp:timetime2 -1.200e+00 1.380e+00 8.100e+01 -0.869 0.3872
## drugb:timetime3 6.200e+00 1.380e+00 8.100e+01 4.492 2.32e-05 ***
## drugp:timetime3 -2.300e+00 1.380e+00 8.100e+01 -1.666 0.0995 .
## drugb:timetime4 -9.000e-01 1.380e+00 8.100e+01 -0.652 0.5162
## drugp:timetime4 -3.500e+00 1.380e+00 8.100e+01 -2.536 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drugb drugp timtm2 timtm3 timtm4 drgb:2 drgp:2 drgb:3
## drugb -0.707
## drugp -0.707 0.500
## timetime2 -0.358 0.253 0.253
## timetime3 -0.358 0.253 0.253 0.500
## timetime4 -0.358 0.253 0.253 0.500 0.500
## drugb:tmtm2 0.253 -0.358 -0.179 -0.707 -0.354 -0.354
## drugp:tmtm2 0.253 -0.179 -0.358 -0.707 -0.354 -0.354 0.500
## drugb:tmtm3 0.253 -0.358 -0.179 -0.354 -0.707 -0.354 0.500 0.250
## drugp:tmtm3 0.253 -0.179 -0.358 -0.354 -0.707 -0.354 0.250 0.500 0.500
## drugb:tmtm4 0.253 -0.358 -0.179 -0.354 -0.354 -0.707 0.500 0.250 0.500
## drugp:tmtm4 0.253 -0.179 -0.358 -0.354 -0.354 -0.707 0.250 0.500 0.250
## drgp:3 drgb:4
## drugb
## drugp
## timetime2
## timetime3
## timetime4
## drugb:tmtm2
## drugp:tmtm2
## drugb:tmtm3
## drugp:tmtm3
## drugb:tmtm4 0.250
## drugp:tmtm4 0.500 0.500
## convergence code: 0
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
anova(myfit2)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## drug 192.89 96.443 2 27 20.247 4.249e-06 ***
## time 222.29 74.097 3 81 15.556 4.444e-08 ***
## drug:time 320.13 53.356 6 81 11.201 4.539e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#same result if using d2.long data
myfit3<-lmer(rate~drug*time+(1|subject),data=d2.long)
summary(myfit3)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: rate ~ drug * time + (1 | subject)
## Data: d2.long
##
## REML criterion at convergence: 571.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39109 -0.52260 -0.05669 0.50612 2.41993
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 13.864 3.723
## Residual 4.763 2.182
## Number of obs: 120, groups: subject, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.290e+01 1.365e+00 4.057e+01 60.741 < 2e-16 ***
## drugb -1.000e+01 1.930e+00 4.057e+01 -5.181 6.41e-06 ***
## drugp -9.100e+00 1.930e+00 4.057e+01 -4.715 2.85e-05 ***
## timetime2 2.000e-01 9.760e-01 8.100e+01 0.205 0.8382
## timetime3 5.000e-01 9.760e-01 8.100e+01 0.512 0.6099
## timetime4 1.057e-14 9.760e-01 8.100e+01 0.000 1.0000
## drugb:timetime2 5.900e+00 1.380e+00 8.100e+01 4.274 5.20e-05 ***
## drugp:timetime2 -1.200e+00 1.380e+00 8.100e+01 -0.869 0.3872
## drugb:timetime3 6.200e+00 1.380e+00 8.100e+01 4.492 2.32e-05 ***
## drugp:timetime3 -2.300e+00 1.380e+00 8.100e+01 -1.666 0.0995 .
## drugb:timetime4 -9.000e-01 1.380e+00 8.100e+01 -0.652 0.5162
## drugp:timetime4 -3.500e+00 1.380e+00 8.100e+01 -2.536 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drugb drugp timtm2 timtm3 timtm4 drgb:2 drgp:2 drgb:3
## drugb -0.707
## drugp -0.707 0.500
## timetime2 -0.358 0.253 0.253
## timetime3 -0.358 0.253 0.253 0.500
## timetime4 -0.358 0.253 0.253 0.500 0.500
## drugb:tmtm2 0.253 -0.358 -0.179 -0.707 -0.354 -0.354
## drugp:tmtm2 0.253 -0.179 -0.358 -0.707 -0.354 -0.354 0.500
## drugb:tmtm3 0.253 -0.358 -0.179 -0.354 -0.707 -0.354 0.500 0.250
## drugp:tmtm3 0.253 -0.179 -0.358 -0.354 -0.707 -0.354 0.250 0.500 0.500
## drugb:tmtm4 0.253 -0.358 -0.179 -0.354 -0.354 -0.707 0.500 0.250 0.500
## drugp:tmtm4 0.253 -0.179 -0.358 -0.354 -0.354 -0.707 0.250 0.500 0.250
## drgp:3 drgb:4
## drugb
## drugp
## timetime2
## timetime3
## timetime4
## drugb:tmtm2
## drugp:tmtm2
## drugb:tmtm3
## drugp:tmtm3
## drugb:tmtm4 0.250
## drugp:tmtm4 0.500 0.500
anova(myfit3)
## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## drug 192.89 96.443 2 27 20.247 4.249e-06 ***
## time 222.29 74.097 3 81 15.556 4.444e-08 ***
## drug:time 320.13 53.356 6 81 11.201 4.539e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#however, if within each group p, a and b, subjects are listed as 1-10, 1-10, 1-10 respectively
#you will need to add the nested structure, otherwise, R will consider subject 1 from group p, a and b as the
#same subject
#d3 is constructed as "subjects are listed as 1-10, 1-10, 1-10 respectively within the three groups"
d3<-read.table("C:/jenn/teaching/stat579/data/splitd3.txt",header=TRUE)
d3.long <- melt(d3,
# id.vars: ID variables
# all variables to keep but not split apart on
id.vars=c("subject","drug"),
# measure.vars: The source columns
# (if unspecified then all other variables are measure.vars)
measure.vars = c("time1", "time2", "time3", "time4"),
# variable.name: Name of the destination column identifying each
# original column that the measurement came from
variable.name = "time",
# value.name: column name for values in table
value.name = "rate"
)
nrow(d3)
## [1] 30
d3
## drug subject time1 time2 time3 time4
## 1 p 1 80 77 73 69
## 2 p 2 64 66 68 71
## 3 p 3 75 73 73 69
## 4 p 4 72 70 74 73
## 5 p 5 74 74 71 67
## 6 p 6 71 71 72 70
## 7 p 7 76 78 74 71
## 8 p 8 73 68 64 64
## 9 p 9 76 73 74 76
## 10 p 10 77 78 77 73
## 11 a 1 81 81 82 82
## 12 a 2 82 83 80 81
## 13 a 3 81 77 80 80
## 14 a 4 84 86 85 85
## 15 a 5 88 90 88 86
## 16 a 6 83 82 86 85
## 17 a 7 85 83 87 86
## 18 a 8 81 85 86 85
## 19 a 9 87 89 87 82
## 20 a 10 77 75 73 77
## 21 b 1 76 83 85 79
## 22 b 2 75 81 85 73
## 23 b 3 75 82 80 77
## 24 b 4 68 73 72 69
## 25 b 5 78 87 86 77
## 26 b 6 81 85 81 74
## 27 b 7 67 73 75 66
## 28 b 8 68 73 73 66
## 29 b 9 68 75 79 69
## 30 b 10 73 78 80 70
d3.long[1:60,]
## subject drug time rate
## 1 1 p time1 80
## 2 2 p time1 64
## 3 3 p time1 75
## 4 4 p time1 72
## 5 5 p time1 74
## 6 6 p time1 71
## 7 7 p time1 76
## 8 8 p time1 73
## 9 9 p time1 76
## 10 10 p time1 77
## 11 1 a time1 81
## 12 2 a time1 82
## 13 3 a time1 81
## 14 4 a time1 84
## 15 5 a time1 88
## 16 6 a time1 83
## 17 7 a time1 85
## 18 8 a time1 81
## 19 9 a time1 87
## 20 10 a time1 77
## 21 1 b time1 76
## 22 2 b time1 75
## 23 3 b time1 75
## 24 4 b time1 68
## 25 5 b time1 78
## 26 6 b time1 81
## 27 7 b time1 67
## 28 8 b time1 68
## 29 9 b time1 68
## 30 10 b time1 73
## 31 1 p time2 77
## 32 2 p time2 66
## 33 3 p time2 73
## 34 4 p time2 70
## 35 5 p time2 74
## 36 6 p time2 71
## 37 7 p time2 78
## 38 8 p time2 68
## 39 9 p time2 73
## 40 10 p time2 78
## 41 1 a time2 81
## 42 2 a time2 83
## 43 3 a time2 77
## 44 4 a time2 86
## 45 5 a time2 90
## 46 6 a time2 82
## 47 7 a time2 83
## 48 8 a time2 85
## 49 9 a time2 89
## 50 10 a time2 75
## 51 1 b time2 83
## 52 2 b time2 81
## 53 3 b time2 82
## 54 4 b time2 73
## 55 5 b time2 87
## 56 6 b time2 85
## 57 7 b time2 73
## 58 8 b time2 73
## 59 9 b time2 75
## 60 10 b time2 78
#in this case, without nested structure, R analyze as a completely randomized block design
myfit4<-lmer(rate~drug*time+(1|subject),data=d3.long)
summary(myfit4)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: rate ~ drug * time + (1 | subject)
## Data: d3.long
##
## REML criterion at convergence: 644.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3605 -0.6893 0.1336 0.7202 1.7950
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 2.391 1.546
## Residual 16.236 4.029
## Number of obs: 120, groups: subject, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.290e+01 1.365e+00 9.143e+01 60.741 < 2e-16 ***
## drugb -1.000e+01 1.802e+00 9.900e+01 -5.549 2.40e-07 ***
## drugp -9.100e+00 1.802e+00 9.900e+01 -5.050 2.02e-06 ***
## timetime2 2.000e-01 1.802e+00 9.900e+01 0.111 0.9119
## timetime3 5.000e-01 1.802e+00 9.900e+01 0.277 0.7820
## timetime4 -5.256e-14 1.802e+00 9.900e+01 0.000 1.0000
## drugb:timetime2 5.900e+00 2.548e+00 9.900e+01 2.315 0.0227 *
## drugp:timetime2 -1.200e+00 2.548e+00 9.900e+01 -0.471 0.6388
## drugb:timetime3 6.200e+00 2.548e+00 9.900e+01 2.433 0.0168 *
## drugp:timetime3 -2.300e+00 2.548e+00 9.900e+01 -0.903 0.3690
## drugb:timetime4 -9.000e-01 2.548e+00 9.900e+01 -0.353 0.7247
## drugp:timetime4 -3.500e+00 2.548e+00 9.900e+01 -1.373 0.1727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drugb drugp timtm2 timtm3 timtm4 drgb:2 drgp:2 drgb:3
## drugb -0.660
## drugp -0.660 0.500
## timetime2 -0.660 0.500 0.500
## timetime3 -0.660 0.500 0.500 0.500
## timetime4 -0.660 0.500 0.500 0.500 0.500
## drugb:tmtm2 0.467 -0.707 -0.354 -0.707 -0.354 -0.354
## drugp:tmtm2 0.467 -0.354 -0.707 -0.707 -0.354 -0.354 0.500
## drugb:tmtm3 0.467 -0.707 -0.354 -0.354 -0.707 -0.354 0.500 0.250
## drugp:tmtm3 0.467 -0.354 -0.707 -0.354 -0.707 -0.354 0.250 0.500 0.500
## drugb:tmtm4 0.467 -0.707 -0.354 -0.354 -0.354 -0.707 0.500 0.250 0.500
## drugp:tmtm4 0.467 -0.354 -0.707 -0.354 -0.354 -0.707 0.250 0.500 0.250
## drgp:3 drgb:4
## drugb
## drugp
## timetime2
## timetime3
## timetime4
## drugb:tmtm2
## drugp:tmtm2
## drugb:tmtm3
## drugp:tmtm3
## drugb:tmtm4 0.250
## drugp:tmtm4 0.500 0.500
anova(myfit4)
## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## drug 2438.47 1219.23 2 99 75.095 < 2.2e-16 ***
## time 222.29 74.10 3 99 4.564 0.004886 **
## drug:time 320.13 53.36 6 99 3.286 0.005444 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#adding nested structure
myfit5<-lmer(rate~drug*time+(1|subject/drug),data=d3.long) #subject nested in drug
summary(myfit5)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: rate ~ drug * time + (1 | subject/drug)
## Data: d3.long
##
## REML criterion at convergence: 571.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39109 -0.52260 -0.05669 0.50612 2.41993
##
## Random effects:
## Groups Name Variance Std.Dev.
## drug:subject (Intercept) 13.864 3.723
## subject (Intercept) 0.000 0.000
## Residual 4.763 2.182
## Number of obs: 120, groups: drug:subject, 30; subject, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.290e+01 1.365e+00 4.057e+01 60.741 < 2e-16 ***
## drugb -1.000e+01 1.930e+00 4.057e+01 -5.181 6.41e-06 ***
## drugp -9.100e+00 1.930e+00 4.057e+01 -4.715 2.85e-05 ***
## timetime2 2.000e-01 9.760e-01 8.100e+01 0.205 0.8382
## timetime3 5.000e-01 9.760e-01 8.100e+01 0.512 0.6099
## timetime4 -7.512e-14 9.760e-01 8.100e+01 0.000 1.0000
## drugb:timetime2 5.900e+00 1.380e+00 8.100e+01 4.274 5.20e-05 ***
## drugp:timetime2 -1.200e+00 1.380e+00 8.100e+01 -0.869 0.3872
## drugb:timetime3 6.200e+00 1.380e+00 8.100e+01 4.492 2.32e-05 ***
## drugp:timetime3 -2.300e+00 1.380e+00 8.100e+01 -1.666 0.0995 .
## drugb:timetime4 -9.000e-01 1.380e+00 8.100e+01 -0.652 0.5162
## drugp:timetime4 -3.500e+00 1.380e+00 8.100e+01 -2.536 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drugb drugp timtm2 timtm3 timtm4 drgb:2 drgp:2 drgb:3
## drugb -0.707
## drugp -0.707 0.500
## timetime2 -0.358 0.253 0.253
## timetime3 -0.358 0.253 0.253 0.500
## timetime4 -0.358 0.253 0.253 0.500 0.500
## drugb:tmtm2 0.253 -0.358 -0.179 -0.707 -0.354 -0.354
## drugp:tmtm2 0.253 -0.179 -0.358 -0.707 -0.354 -0.354 0.500
## drugb:tmtm3 0.253 -0.358 -0.179 -0.354 -0.707 -0.354 0.500 0.250
## drugp:tmtm3 0.253 -0.179 -0.358 -0.354 -0.707 -0.354 0.250 0.500 0.500
## drugb:tmtm4 0.253 -0.358 -0.179 -0.354 -0.354 -0.707 0.500 0.250 0.500
## drugp:tmtm4 0.253 -0.179 -0.358 -0.354 -0.354 -0.707 0.250 0.500 0.250
## drgp:3 drgb:4
## drugb
## drugp
## timetime2
## timetime3
## timetime4
## drugb:tmtm2
## drugp:tmtm2
## drugb:tmtm3
## drugp:tmtm3
## drugb:tmtm4 0.250
## drugp:tmtm4 0.500 0.500
anova(myfit5)
## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## drug 192.89 96.443 2 27 20.247 4.249e-06 ***
## time 222.29 74.097 3 81 15.556 4.444e-08 ***
## drug:time 320.13 53.356 6 81 11.201 4.539e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##wrong degrees of freedom from nlme package again, note that the degrees of
#freedom for testing drug effect is 18 instead of 27
myfit6<-lme(rate~drug*time,random=~1|subject/drug,data=d3.long)
summary(myfit6)
## Linear mixed-effects model fit by REML
## Data: d3.long
## AIC BIC logLik
## 601.2025 641.4345 -285.6013
##
## Random effects:
## Formula: ~1 | subject
## (Intercept)
## StdDev: 0.0004930058
##
## Formula: ~1 | drug %in% subject
## (Intercept) Residual
## StdDev: 3.723383 2.182492
##
## Fixed effects: rate ~ drug * time
## Value Std.Error DF t-value p-value
## (Intercept) 82.9 1.3648023 81 60.74140 0.0000
## drugb -10.0 1.9301219 18 -5.18102 0.0001
## drugp -9.1 1.9301219 18 -4.71473 0.0002
## timetime2 0.2 0.9760401 81 0.20491 0.8382
## timetime3 0.5 0.9760401 81 0.51227 0.6099
## timetime4 0.0 0.9760401 81 0.00000 1.0000
## drugb:timetime2 5.9 1.3803292 81 4.27434 0.0001
## drugp:timetime2 -1.2 1.3803292 81 -0.86936 0.3872
## drugb:timetime3 6.2 1.3803292 81 4.49168 0.0000
## drugp:timetime3 -2.3 1.3803292 81 -1.66627 0.0995
## drugb:timetime4 -0.9 1.3803292 81 -0.65202 0.5162
## drugp:timetime4 -3.5 1.3803292 81 -2.53563 0.0131
## Correlation:
## (Intr) drugb drugp timtm2 timtm3 timtm4 drgb:2 drgp:2
## drugb -0.707
## drugp -0.707 0.500
## timetime2 -0.358 0.253 0.253
## timetime3 -0.358 0.253 0.253 0.500
## timetime4 -0.358 0.253 0.253 0.500 0.500
## drugb:timetime2 0.253 -0.358 -0.179 -0.707 -0.354 -0.354
## drugp:timetime2 0.253 -0.179 -0.358 -0.707 -0.354 -0.354 0.500
## drugb:timetime3 0.253 -0.358 -0.179 -0.354 -0.707 -0.354 0.500 0.250
## drugp:timetime3 0.253 -0.179 -0.358 -0.354 -0.707 -0.354 0.250 0.500
## drugb:timetime4 0.253 -0.358 -0.179 -0.354 -0.354 -0.707 0.500 0.250
## drugp:timetime4 0.253 -0.179 -0.358 -0.354 -0.354 -0.707 0.250 0.500
## drgb:3 drgp:3 drgb:4
## drugb
## drugp
## timetime2
## timetime3
## timetime4
## drugb:timetime2
## drugp:timetime2
## drugb:timetime3
## drugp:timetime3 0.500
## drugb:timetime4 0.500 0.250
## drugp:timetime4 0.250 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.39108657 -0.52260081 -0.05669221 0.50611907 2.41992750
##
## Number of Observations: 120
## Number of Groups:
## subject drug %in% subject
## 10 30
anova(myfit6)
## numDF denDF F-value p-value
## (Intercept) 1 81 11833.060 <.0001
## drug 2 18 20.247 <.0001
## time 3 81 15.556 <.0001
## drug:time 6 81 11.201 <.0001
#We could check of if the random effect is important or not, by the step
#function of lmerTest: (usually, people would prefer to include the random effect
#in the model whether they are signnificant or not)
mystep <- step(myfit) #subject is significant, need to keep
mystep$rand.table
## Chi.sq Chi.DF elim.num p.value
## subject 78.7765 1 kept 0
mystep$anova.table
## Sum Sq Mean Sq NumDF DenDF F.value elim.num Pr(>F)
## drug 192.8851 96.44257 2 27.00002 20.24713 kept 4.249203e-06
## time 222.2917 74.09722 3 81.00001 15.55595 kept 4.444219e-08
## drug:time 320.1333 53.35556 6 81.00001 11.20145 kept 4.539026e-09
mystep$diffs.lsmeans.table #mutlitple comparisons
## Estimate Standard Error DF t-value
## drug a - b 7.2000 1.7352 27.0 4.15
## drug a - p 10.8500 1.7352 27.0 6.25
## drug b - p 3.6500 1.7352 27.0 2.10
## time time1 - time2 -1.7667 0.5635 81.0 -3.14
## time time1 - time3 -1.8000 0.5635 81.0 -3.19
## time time1 - time4 1.4667 0.5635 81.0 2.60
## time time2 - time3 -0.0333 0.5635 81.0 -0.06
## time time2 - time4 3.2333 0.5635 81.0 5.74
## time time3 - time4 3.2667 0.5635 81.0 5.80
## drug:time a time1 - b time1 10.0000 1.9301 40.6 5.18
## drug:time a time1 - p time1 9.1000 1.9301 40.6 4.71
## drug:time a time1 - a time2 -0.2000 0.9760 81.0 -0.20
## drug:time a time1 - b time2 3.9000 1.9301 40.6 2.02
## drug:time a time1 - p time2 10.1000 1.9301 40.6 5.23
## drug:time a time1 - a time3 -0.5000 0.9760 81.0 -0.51
## drug:time a time1 - b time3 3.3000 1.9301 40.6 1.71
## drug:time a time1 - p time3 10.9000 1.9301 40.6 5.65
## drug:time a time1 - a time4 0.0000 0.9760 81.0 0.00
## drug:time a time1 - b time4 10.9000 1.9301 40.6 5.65
## drug:time a time1 - p time4 12.6000 1.9301 40.6 6.53
## drug:time b time1 - p time1 -0.9000 1.9301 40.6 -0.47
## drug:time b time1 - a time2 -10.2000 1.9301 40.6 -5.28
## drug:time b time1 - b time2 -6.1000 0.9760 81.0 -6.25
## drug:time b time1 - p time2 0.1000 1.9301 40.6 0.05
## drug:time b time1 - a time3 -10.5000 1.9301 40.6 -5.44
## drug:time b time1 - b time3 -6.7000 0.9760 81.0 -6.86
## drug:time b time1 - p time3 0.9000 1.9301 40.6 0.47
## drug:time b time1 - a time4 -10.0000 1.9301 40.6 -5.18
## drug:time b time1 - b time4 0.9000 0.9760 81.0 0.92
## drug:time b time1 - p time4 2.6000 1.9301 40.6 1.35
## drug:time p time1 - a time2 -9.3000 1.9301 40.6 -4.82
## drug:time p time1 - b time2 -5.2000 1.9301 40.6 -2.69
## drug:time p time1 - p time2 1.0000 0.9760 81.0 1.02
## drug:time p time1 - a time3 -9.6000 1.9301 40.6 -4.97
## drug:time p time1 - b time3 -5.8000 1.9301 40.6 -3.00
## drug:time p time1 - p time3 1.8000 0.9760 81.0 1.84
## drug:time p time1 - a time4 -9.1000 1.9301 40.6 -4.71
## drug:time p time1 - b time4 1.8000 1.9301 40.6 0.93
## drug:time p time1 - p time4 3.5000 0.9760 81.0 3.59
## drug:time a time2 - b time2 4.1000 1.9301 40.6 2.12
## drug:time a time2 - p time2 10.3000 1.9301 40.6 5.34
## drug:time a time2 - a time3 -0.3000 0.9760 81.0 -0.31
## drug:time a time2 - b time3 3.5000 1.9301 40.6 1.81
## drug:time a time2 - p time3 11.1000 1.9301 40.6 5.75
## drug:time a time2 - a time4 0.2000 0.9760 81.0 0.20
## drug:time a time2 - b time4 11.1000 1.9301 40.6 5.75
## drug:time a time2 - p time4 12.8000 1.9301 40.6 6.63
## drug:time b time2 - p time2 6.2000 1.9301 40.6 3.21
## drug:time b time2 - a time3 -4.4000 1.9301 40.6 -2.28
## drug:time b time2 - b time3 -0.6000 0.9760 81.0 -0.61
## drug:time b time2 - p time3 7.0000 1.9301 40.6 3.63
## drug:time b time2 - a time4 -3.9000 1.9301 40.6 -2.02
## drug:time b time2 - b time4 7.0000 0.9760 81.0 7.17
## drug:time b time2 - p time4 8.7000 1.9301 40.6 4.51
## drug:time p time2 - a time3 -10.6000 1.9301 40.6 -5.49
## drug:time p time2 - b time3 -6.8000 1.9301 40.6 -3.52
## drug:time p time2 - p time3 0.8000 0.9760 81.0 0.82
## drug:time p time2 - a time4 -10.1000 1.9301 40.6 -5.23
## drug:time p time2 - b time4 0.8000 1.9301 40.6 0.41
## drug:time p time2 - p time4 2.5000 0.9760 81.0 2.56
## drug:time a time3 - b time3 3.8000 1.9301 40.6 1.97
## drug:time a time3 - p time3 11.4000 1.9301 40.6 5.91
## drug:time a time3 - a time4 0.5000 0.9760 81.0 0.51
## drug:time a time3 - b time4 11.4000 1.9301 40.6 5.91
## drug:time a time3 - p time4 13.1000 1.9301 40.6 6.79
## drug:time b time3 - p time3 7.6000 1.9301 40.6 3.94
## drug:time b time3 - a time4 -3.3000 1.9301 40.6 -1.71
## drug:time b time3 - b time4 7.6000 0.9760 81.0 7.79
## drug:time b time3 - p time4 9.3000 1.9301 40.6 4.82
## drug:time p time3 - a time4 -10.9000 1.9301 40.6 -5.65
## drug:time p time3 - b time4 0.0000 1.9301 40.6 0.00
## drug:time p time3 - p time4 1.7000 0.9760 81.0 1.74
## drug:time a time4 - b time4 10.9000 1.9301 40.6 5.65
## drug:time a time4 - p time4 12.6000 1.9301 40.6 6.53
## drug:time b time4 - p time4 1.7000 1.9301 40.6 0.88
## Lower CI Upper CI p-value
## drug a - b 3.6397 10.7603 0.0003
## drug a - p 7.2897 14.4103 0.0000
## drug b - p 0.0897 7.2103 0.0449
## time time1 - time2 -2.8879 -0.6454 0.0024
## time time1 - time3 -2.9212 -0.6788 0.0020
## time time1 - time4 0.3454 2.5879 0.0110
## time time2 - time3 -1.1546 1.0879 0.9530
## time time2 - time4 2.1121 4.3546 0.0000
## time time3 - time4 2.1454 4.3879 0.0000
## drug:time a time1 - b time1 6.1008 13.8992 0.0000
## drug:time a time1 - p time1 5.2008 12.9992 0.0000
## drug:time a time1 - a time2 -2.1420 1.7420 0.8382
## drug:time a time1 - b time2 0.0008 7.7992 0.0500
## drug:time a time1 - p time2 6.2008 13.9992 0.0000
## drug:time a time1 - a time3 -2.4420 1.4420 0.6099
## drug:time a time1 - b time3 -0.5992 7.1992 0.0950
## drug:time a time1 - p time3 7.0008 14.7992 0.0000
## drug:time a time1 - a time4 -1.9420 1.9420 1.0000
## drug:time a time1 - b time4 7.0008 14.7992 0.0000
## drug:time a time1 - p time4 8.7008 16.4992 0.0000
## drug:time b time1 - p time1 -4.7992 2.9992 0.6435
## drug:time b time1 - a time2 -14.0992 -6.3008 0.0000
## drug:time b time1 - b time2 -8.0420 -4.1580 0.0000
## drug:time b time1 - p time2 -3.7992 3.9992 0.9589
## drug:time b time1 - a time3 -14.3992 -6.6008 0.0000
## drug:time b time1 - b time3 -8.6420 -4.7580 0.0000
## drug:time b time1 - p time3 -2.9992 4.7992 0.6435
## drug:time b time1 - a time4 -13.8992 -6.1008 0.0000
## drug:time b time1 - b time4 -1.0420 2.8420 0.3592
## drug:time b time1 - p time4 -1.2992 6.4992 0.1854
## drug:time p time1 - a time2 -13.1992 -5.4008 0.0000
## drug:time p time1 - b time2 -9.0992 -1.3008 0.0102
## drug:time p time1 - p time2 -0.9420 2.9420 0.3086
## drug:time p time1 - a time3 -13.4992 -5.7008 0.0000
## drug:time p time1 - b time3 -9.6992 -1.9008 0.0045
## drug:time p time1 - p time3 -0.1420 3.7420 0.0688
## drug:time p time1 - a time4 -12.9992 -5.2008 0.0000
## drug:time p time1 - b time4 -2.0992 5.6992 0.3566
## drug:time p time1 - p time4 1.5580 5.4420 0.0006
## drug:time a time2 - b time2 0.2008 7.9992 0.0398
## drug:time a time2 - p time2 6.4008 14.1992 0.0000
## drug:time a time2 - a time3 -2.2420 1.6420 0.7594
## drug:time a time2 - b time3 -0.3992 7.3992 0.0772
## drug:time a time2 - p time3 7.2008 14.9992 0.0000
## drug:time a time2 - a time4 -1.7420 2.1420 0.8382
## drug:time a time2 - b time4 7.2008 14.9992 0.0000
## drug:time a time2 - p time4 8.9008 16.6992 0.0000
## drug:time b time2 - p time2 2.3008 10.0992 0.0026
## drug:time b time2 - a time3 -8.2992 -0.5008 0.0280
## drug:time b time2 - b time3 -2.5420 1.3420 0.5405
## drug:time b time2 - p time3 3.1008 10.8992 0.0008
## drug:time b time2 - a time4 -7.7992 -0.0008 0.0500
## drug:time b time2 - b time4 5.0580 8.9420 0.0000
## drug:time b time2 - p time4 4.8008 12.5992 0.0001
## drug:time p time2 - a time3 -14.4992 -6.7008 0.0000
## drug:time p time2 - b time3 -10.6992 -2.9008 0.0011
## drug:time p time2 - p time3 -1.1420 2.7420 0.4148
## drug:time p time2 - a time4 -13.9992 -6.2008 0.0000
## drug:time p time2 - b time4 -3.0992 4.6992 0.6807
## drug:time p time2 - p time4 0.5580 4.4420 0.0123
## drug:time a time3 - b time3 -0.0992 7.6992 0.0558
## drug:time a time3 - p time3 7.5008 15.2992 0.0000
## drug:time a time3 - a time4 -1.4420 2.4420 0.6099
## drug:time a time3 - b time4 7.5008 15.2992 0.0000
## drug:time a time3 - p time4 9.2008 16.9992 0.0000
## drug:time b time3 - p time3 3.7008 11.4992 0.0003
## drug:time b time3 - a time4 -7.1992 0.5992 0.0950
## drug:time b time3 - b time4 5.6580 9.5420 0.0000
## drug:time b time3 - p time4 5.4008 13.1992 0.0000
## drug:time p time3 - a time4 -14.7992 -7.0008 0.0000
## drug:time p time3 - b time4 -3.8992 3.8992 1.0000
## drug:time p time3 - p time4 -0.2420 3.6420 0.0854
## drug:time a time4 - b time4 7.0008 14.7992 0.0000
## drug:time a time4 - p time4 8.7008 16.4992 0.0000
## drug:time b time4 - p time4 -2.1992 5.5992 0.3836
##since interaction is significant, let's do multiple comparisons
lsmeans(myfit,pairwise~drug:time)
## $lsmeans
## drug time lsmean SE df lower.CL upper.CL
## a time1 82.9 1.364802 40.57 80.14285 85.65715
## b time1 72.9 1.364802 40.57 70.14285 75.65715
## p time1 73.8 1.364802 40.57 71.04285 76.55715
## a time2 83.1 1.364802 40.57 80.34285 85.85715
## b time2 79.0 1.364802 40.57 76.24285 81.75715
## p time2 72.8 1.364802 40.57 70.04285 75.55715
## a time3 83.4 1.364802 40.57 80.64285 86.15715
## b time3 79.6 1.364802 40.57 76.84285 82.35715
## p time3 72.0 1.364802 40.57 69.24285 74.75715
## a time4 82.9 1.364802 40.57 80.14285 85.65715
## b time4 72.0 1.364802 40.57 69.24285 74.75715
## p time4 70.3 1.364802 40.57 67.54285 73.05715
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## a,time1 - b,time1 1.000000e+01 1.9301219 40.57 5.181 0.0004
## a,time1 - p,time1 9.100000e+00 1.9301219 40.57 4.715 0.0015
## a,time1 - a,time2 -2.000000e-01 0.9760401 81.00 -0.205 1.0000
## a,time1 - b,time2 3.900000e+00 1.9301219 40.57 2.021 0.6775
## a,time1 - p,time2 1.010000e+01 1.9301219 40.57 5.233 0.0003
## a,time1 - a,time3 -5.000000e-01 0.9760401 81.00 -0.512 1.0000
## a,time1 - b,time3 3.300000e+00 1.9301219 40.57 1.710 0.8534
## a,time1 - p,time3 1.090000e+01 1.9301219 40.57 5.647 0.0001
## a,time1 - a,time4 -1.057053e-14 0.9760401 81.00 0.000 1.0000
## a,time1 - b,time4 1.090000e+01 1.9301219 40.57 5.647 0.0001
## a,time1 - p,time4 1.260000e+01 1.9301219 40.57 6.528 <.0001
## b,time1 - p,time1 -9.000000e-01 1.9301219 40.57 -0.466 1.0000
## b,time1 - a,time2 -1.020000e+01 1.9301219 40.57 -5.285 0.0003
## b,time1 - b,time2 -6.100000e+00 0.9760401 81.00 -6.250 <.0001
## b,time1 - p,time2 1.000000e-01 1.9301219 40.57 0.052 1.0000
## b,time1 - a,time3 -1.050000e+01 1.9301219 40.57 -5.440 0.0002
## b,time1 - b,time3 -6.700000e+00 0.9760401 81.00 -6.864 <.0001
## b,time1 - p,time3 9.000000e-01 1.9301219 40.57 0.466 1.0000
## b,time1 - a,time4 -1.000000e+01 1.9301219 40.57 -5.181 0.0004
## b,time1 - b,time4 9.000000e-01 0.9760401 81.00 0.922 0.9987
## b,time1 - p,time4 2.600000e+00 1.9301219 40.57 1.347 0.9670
## p,time1 - a,time2 -9.300000e+00 1.9301219 40.57 -4.818 0.0011
## p,time1 - b,time2 -5.200000e+00 1.9301219 40.57 -2.694 0.2654
## p,time1 - p,time2 1.000000e+00 0.9760401 81.00 1.025 0.9967
## p,time1 - a,time3 -9.600000e+00 1.9301219 40.57 -4.974 0.0007
## p,time1 - b,time3 -5.800000e+00 1.9301219 40.57 -3.005 0.1443
## p,time1 - p,time3 1.800000e+00 0.9760401 81.00 1.844 0.7886
## p,time1 - a,time4 -9.100000e+00 1.9301219 40.57 -4.715 0.0015
## p,time1 - b,time4 1.800000e+00 1.9301219 40.57 0.933 0.9983
## p,time1 - p,time4 3.500000e+00 0.9760401 81.00 3.586 0.0266
## a,time2 - b,time2 4.100000e+00 1.9301219 40.57 2.124 0.6095
## a,time2 - p,time2 1.030000e+01 1.9301219 40.57 5.336 0.0002
## a,time2 - a,time3 -3.000000e-01 0.9760401 81.00 -0.307 1.0000
## a,time2 - b,time3 3.500000e+00 1.9301219 40.57 1.813 0.8014
## a,time2 - p,time3 1.110000e+01 1.9301219 40.57 5.751 0.0001
## a,time2 - a,time4 2.000000e-01 0.9760401 81.00 0.205 1.0000
## a,time2 - b,time4 1.110000e+01 1.9301219 40.57 5.751 0.0001
## a,time2 - p,time4 1.280000e+01 1.9301219 40.57 6.632 <.0001
## b,time2 - p,time2 6.200000e+00 1.9301219 40.57 3.212 0.0914
## b,time2 - a,time3 -4.400000e+00 1.9301219 40.57 -2.280 0.5061
## b,time2 - b,time3 -6.000000e-01 0.9760401 81.00 -0.615 1.0000
## b,time2 - p,time3 7.000000e+00 1.9301219 40.57 3.627 0.0332
## b,time2 - a,time4 -3.900000e+00 1.9301219 40.57 -2.021 0.6775
## b,time2 - b,time4 7.000000e+00 0.9760401 81.00 7.172 <.0001
## b,time2 - p,time4 8.700000e+00 1.9301219 40.57 4.507 0.0028
## p,time2 - a,time3 -1.060000e+01 1.9301219 40.57 -5.492 0.0001
## p,time2 - b,time3 -6.800000e+00 1.9301219 40.57 -3.523 0.0432
## p,time2 - p,time3 8.000000e-01 0.9760401 81.00 0.820 0.9996
## p,time2 - a,time4 -1.010000e+01 1.9301219 40.57 -5.233 0.0003
## p,time2 - b,time4 8.000000e-01 1.9301219 40.57 0.414 1.0000
## p,time2 - p,time4 2.500000e+00 0.9760401 81.00 2.561 0.3188
## a,time3 - b,time3 3.800000e+00 1.9301219 40.57 1.969 0.7104
## a,time3 - p,time3 1.140000e+01 1.9301219 40.57 5.906 <.0001
## a,time3 - a,time4 5.000000e-01 0.9760401 81.00 0.512 1.0000
## a,time3 - b,time4 1.140000e+01 1.9301219 40.57 5.906 <.0001
## a,time3 - p,time4 1.310000e+01 1.9301219 40.57 6.787 <.0001
## b,time3 - p,time3 7.600000e+00 1.9301219 40.57 3.938 0.0145
## b,time3 - a,time4 -3.300000e+00 1.9301219 40.57 -1.710 0.8534
## b,time3 - b,time4 7.600000e+00 0.9760401 81.00 7.787 <.0001
## b,time3 - p,time4 9.300000e+00 1.9301219 40.57 4.818 0.0011
## p,time3 - a,time4 -1.090000e+01 1.9301219 40.57 -5.647 0.0001
## p,time3 - b,time4 6.547974e-14 1.9301219 40.57 0.000 1.0000
## p,time3 - p,time4 1.700000e+00 0.9760401 81.00 1.742 0.8432
## a,time4 - b,time4 1.090000e+01 1.9301219 40.57 5.647 0.0001
## a,time4 - p,time4 1.260000e+01 1.9301219 40.57 6.528 <.0001
## b,time4 - p,time4 1.700000e+00 1.9301219 40.57 0.881 0.9990
##
## P value adjustment: tukey method for comparing a family of 12 estimates
##diagnostics
fit.lmer<-predict(myfit)
resid.lmer<-resid(myfit)
par(mfrow=c(2,2))
plot(resid.lmer~fit.lmer,main="fitted values vs. residuals",xlab="fitted",ylab="residuals")
abline(0,0,lty=2)
qqnorm(resid.lmer)
qqline(resid.lmer)
hist(resid.lmer,main="histogram of residuals",xlab="residuals")
#random effects normality
par(mfrow=c(1,2))
alphahat<-ranef(myfit)
hist(alphahat$subject[,1])
abline(0,0,lty=2)
qqnorm(alphahat$subject[,1])
qqline(alphahat$subject[,1])
