Lab 2 Longitudinal Models
ICPSR Summer Program at York University
Georges Monette
May 2018
(Updated: May 28 2018 15:04)
Contents:
- LME model
- Hausman test:
- Adjusting for time
- Diagnostics: Level 1
- Diagnostics for heteroskedasticity
- Diagnostics for autocorrelation
- Diagnostics: Level 2
- Dropping observations
- Modeling autocorrelation
- Modeling heteroskedasticity
- Interpreting different kinds of residual plots
- Visualizing the impact of model selection
- Displaying data and fitted values together
#
# Schizophrenia patients were assessed at 6 annual checkups.
#
# At each checkup, the type of drug currently prescribed by the treating
# physician was recorded. There were 3 types of drugs: Typical, Atypical and
# Clozapine.
#
# Researchers are interested in using this observational data to study the
# effectiveness of the newest drug, Clozapine, particularly comparing it with
# 'Typical' drugs, the established standard.
#
# There are a number of 'outcome' measures (i.e. measures of the
# severity of the illness) but only one, 'neg', is of interest
# for now.
#
library(spida2) # devtools::install_github('gmonette/spida2')
library(p3d) # devtools::install_github('gmonette/p3d')| Loading required package: rgl|
| Attaching package: 'rgl'| The following object is masked from 'package:spida2':
|
| %>%|
| Attaching package: 'p3d'| The following objects are masked from 'package:spida2':
|
| cell, center, ConjComp, dell, disp, ell, ell.conj, ellbox,
| ellplus, ellpt, ellptc, ellpts, ellptsc, elltan, elltanc,
| na.include, uvlibrary(car)| Loading required package: carData|
| Attaching package: 'car'| The following object is masked from 'package:p3d':
|
| Identify3dlibrary(lattice)
library(latticeExtra)| Loading required package: RColorBrewerinteractive <- FALSE
# ?Drugs
some( Drugs )| age1 yrsill1 pos neg gen total year drug status Sex
| 22 26.2 9 7 12 25 44 1 Typical Married F
| 326 30.5 10 20 19 39 78 6 Clozapine Single F
| 330 31.1 12 10 20 33 63 6 Atypical Separated/Divorced F
| 164 31.8 13 7 11 18 36 3 Clozapine Single M
| 188 31.3 6 8 14 42 64 4 Typical Single M
| 298 31.3 6 9 7 24 40 6 Typical Single M
| 294 33.4 14 13 22 37 72 6 Atypical Single M
| 128 31.3 15 17 12 27 56 3 Typical Single M
| 238 31.3 15 18 12 28 58 5 Clozapine Single M
| 75 37.7 15 26 21 43 90 2 Clozapine Single M
| Subj
| 22 F10
| 326 F54
| 330 F55
| 164 M29
| 188 M13
| 298 M13
| 294 M43
| 128 M37
| 238 M37
| 75 M47# Which drug seems best to reduce 'neg' symptoms
xyplot( neg ~ year | Subj, Drugs)
# smaller strips
xyplot( neg ~ year | Subj, Drugs,
par.strip.text = list(cex = .7))
dd <- Drugs
dd$id <- reorder( dd$Subj, dd$neg) # orders 'id's in order of average 'neg' symptoms
# choose nice color
pals('green')
pals('red')
display.brewer.all()
cols <- brewer.pal(3,'Dark2')
cols| [1] "#1B9E77" "#D95F02" "#7570B3"pal(cols)
| red green blue
| #1B9E77 27 158 119
| #D95F02 217 95 2
| #7570B3 117 112 179gd(col = cols) # uses ggplot2-like style
gd(pch = c(16,16), cex = 1)
xyplot( neg ~ year | id, dd , groups = drug,
auto.key = list(columns=3),
par.strip.text = list(cex = .7))
#
# QUESTION:
#
# Can we perform OLS fits on each cluster?
# Note that the data are balanced with respect to time
# but not with respect to drugs.
#
# Also note that Clozapine is more frequently given later in the study
#
# compare drugs:
# Pooling the data
fit.lm <- lm(neg ~ drug, dd)
summary(fit.lm)|
| Call:
| lm(formula = neg ~ drug, data = dd)
|
| Residuals:
| Min 1Q Median 3Q Max
| -9.299 -5.233 -1.233 3.767 18.767
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 15.27778 0.71557 21.351 <2e-16 ***
| drugClozapine 1.02107 0.96737 1.056 0.292
| drugTypical -0.04507 0.86250 -0.052 0.958
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 6.072 on 315 degrees of freedom
| Multiple R-squared: 0.005996, Adjusted R-squared: -0.0003149
| F-statistic: 0.9501 on 2 and 315 DF, p-value: 0.3878Ld <- Ldiff(fit.lm, 'drug')
wald(fit.lm, Ld)| numDF denDF F.value p.value
| 1 2 315 0.9501109 0.3878
|
| Estimate Std.Error DF t-value p-value Lower 0.95
| Clozapine - <ref> 1.021073 0.967366 315 1.055519 0.29200 -0.882243
| Typical - <ref> -0.045073 0.862500 315 -0.052259 0.95836 -1.742063
| Typical - Clozapine -1.066146 0.809706 315 -1.316707 0.18889 -2.659262
|
| Upper 0.95
| Clozapine - <ref> 2.924388
| Typical - <ref> 1.651917
| Typical - Clozapine 0.526970Ld <- Ldiff(fit.lm, 'drug', ref = "Atypical")
wald(fit.lm, Ld)| numDF denDF F.value p.value
| 1 2 315 0.9501109 0.3878
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical 1.021073 0.967366 315 1.055519 0.29200
| Typical - Atypical -0.045073 0.862500 315 -0.052259 0.95836
| Typical - Clozapine -1.066146 0.809706 315 -1.316707 0.18889
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -0.882243 2.924388
| Typical - Atypical -1.742063 1.651917
| Typical - Clozapine -2.659262 0.526970# Which looks better: Clozapine or Typical? Remember that a higher response means
# more severe symptoms.
#
# LME model
#
library(nlme)|
| Attaching package: 'nlme'| The following object is masked from 'package:spida2':
|
| getDatafit <- lme( neg ~ drug, dd, random = ~ 1 | id )
summary(fit)| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1829.011 1847.773 -909.5053
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 5.240588 3.399672
|
| Fixed effects: neg ~ drug
| Value Std.Error DF t-value p-value
| (Intercept) 15.482816 0.8783362 263 17.627437 0.0000
| drugClozapine -1.300272 0.7980602 263 -1.629290 0.1044
| drugTypical 0.815021 0.6140648 263 1.327255 0.1856
| Correlation:
| (Intr) drgClz
| drugClozapine -0.444
| drugTypical -0.489 0.559
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.71174345 -0.59155383 -0.05454849 0.47584480 3.47145299
|
| Number of Observations: 318
| Number of Groups: 53wald(fit, -1)| numDF denDF F.value p.value
| 1 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value Lower 0.95
| drugClozapine -1.300272 0.798060 263 -1.629290 0.10445 -2.871672
| drugTypical 0.815021 0.614065 263 1.327255 0.18558 -0.394089
|
| Upper 0.95
| drugClozapine 0.271129
| drugTypical 2.024130# Note: overall significance although individual p-values were not
# Now, Clozapine looks best
Ld <- Ldiff ( fit, "drug") # hypothesis matrix to test differences between drugs
Ld| [,1] [,2] [,3]
| Clozapine - <ref> 0 1 0
| Typical - <ref> 0 0 1
| Typical - Clozapine 0 -1 1wald( fit, Ld )| numDF denDF F.value p.value
| 1 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value Lower 0.95
| Clozapine - <ref> -1.300272 0.798060 263 -1.629290 0.10445 -2.871672
| Typical - <ref> 0.815021 0.614065 263 1.327255 0.18558 -0.394089
| Typical - Clozapine 2.115292 0.682395 263 3.099806 0.00215 0.771639
|
| Upper 0.95
| Clozapine - <ref> 0.271129
| Typical - <ref> 2.024130
| Typical - Clozapine 3.458945Ld <- Ldiff ( fit, "drug", ref = "Atypical")
wald( fit, Ld )| numDF denDF F.value p.value
| 1 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.300272 0.798060 263 -1.629290 0.10445
| Typical - Atypical 0.815021 0.614065 263 1.327255 0.18558
| Typical - Clozapine 2.115292 0.682395 263 3.099806 0.00215
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.871672 0.271129
| Typical - Atypical -0.394089 2.024130
| Typical - Clozapine 0.771639 3.458945# Other ways
# QUESTIONS:
# 1) Can you explain why 'drug' is significant overall although
# neither p-value is significant in the summary output. Note that
# this illustrates the danger in dropping multiple terms on the
# basis of individual p-values.
#
# 2) How is clozapine doing ?
# Is it better than the others as expected?
# What explanation can you think of for the results?
#
#
# Hausman test: test whether contextual variables have a significant effect.
#
fit2 <- update( fit, . ~ . + cvar( drug, id))
summary( fit2 )| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1818.357 1844.581 -902.1785
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 4.984062 3.39656
|
| Fixed effects: neg ~ drug + cvar(drug, id)
| Value Std.Error DF t-value p-value
| (Intercept) 15.838657 2.121548 263 7.465614 0.0000
| drugClozapine -1.782102 0.834395 263 -2.135802 0.0336
| drugTypical 0.896756 0.627306 263 1.429535 0.1540
| cvar(drug, id)Clozapine 3.813141 2.872332 50 1.327542 0.1904
| cvar(drug, id)Typical -2.616211 2.891131 50 -0.904909 0.3699
| Correlation:
| (Intr) drgClz drgTyp c(,i)C
| drugClozapine 0.000
| drugTypical 0.000 0.552
| cvar(drug, id)Clozapine -0.771 -0.290 -0.160
| cvar(drug, id)Typical -0.884 -0.120 -0.217 0.668
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.78234220 -0.60736018 -0.06230846 0.49922986 3.32999988
|
| Number of Observations: 318
| Number of Groups: 53wald( fit2, 'cvar')| numDF denDF F.value p.value
| cvar 2 50 3.783136 0.02952
|
| Coefficients Estimate Std.Error DF t-value p-value
| cvar(drug, id)Clozapine 3.813141 2.872332 50 1.327542 0.19036
| cvar(drug, id)Typical -2.616211 2.891131 50 -0.904909 0.36985
|
| Coefficients Lower 0.95 Upper 0.95
| cvar(drug, id)Clozapine -1.956108 9.582391
| cvar(drug, id)Typical -8.423218 3.190796#
# QUESTIONS:
# 1) What is the interpretation of the contextual variable for a
# categorical effect:
head( cbind( dd['id'], getX( fit2) ), 18 )| id (Intercept) drugClozapine drugTypical cvar(drug, id)Clozapine
| 17 F36 1 0 1 0.6666667
| 72 F36 1 0 0 0.6666667
| 127 F36 1 1 0 0.6666667
| 182 F36 1 1 0 0.6666667
| 237 F36 1 1 0 0.6666667
| 292 F36 1 1 0 0.6666667
| 38 F41 1 0 1 0.6666667
| 93 F41 1 0 0 0.6666667
| 148 F41 1 1 0 0.6666667
| 203 F41 1 1 0 0.6666667
| 258 F41 1 1 0 0.6666667
| 313 F41 1 1 0 0.6666667
| 14 F2 1 0 1 0.0000000
| 69 F2 1 0 1 0.0000000
| 124 F2 1 0 1 0.0000000
| 179 F2 1 0 1 0.0000000
| 234 F2 1 0 1 0.0000000
| 289 F2 1 0 1 0.0000000
| cvar(drug, id)Typical
| 17 0.1666667
| 72 0.1666667
| 127 0.1666667
| 182 0.1666667
| 237 0.1666667
| 292 0.1666667
| 38 0.1666667
| 93 0.1666667
| 148 0.1666667
| 203 0.1666667
| 258 0.1666667
| 313 0.1666667
| 14 1.0000000
| 69 1.0000000
| 124 1.0000000
| 179 1.0000000
| 234 1.0000000
| 289 1.0000000#
# 2) Compare results with and without cvar's
#
# 3) Compare cvar effect with corresponding raw var: what does this suggest??
# E.g. 'Typical' and cvar(Typical)?
#
wald( fit2, -1)| numDF denDF F.value p.value
| 1 4 50 4.294531 0.0046
|
| Estimate Std.Error DF t-value p-value
| drugClozapine -1.782102 0.834395 263 -2.135802 0.03362
| drugTypical 0.896756 0.627306 263 1.429535 0.15404
| cvar(drug, id)Clozapine 3.813141 2.872332 50 1.327542 0.19036
| cvar(drug, id)Typical -2.616211 2.891131 50 -0.904909 0.36985
|
| Lower 0.95 Upper 0.95
| drugClozapine -3.425047 -0.139158
| drugTypical -0.338425 2.131937
| cvar(drug, id)Clozapine -1.956108 9.582391
| cvar(drug, id)Typical -8.423218 3.190796wald( fit , -1)| numDF denDF F.value p.value
| 1 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value Lower 0.95
| drugClozapine -1.300272 0.798060 263 -1.629290 0.10445 -2.871672
| drugTypical 0.815021 0.614065 263 1.327255 0.18558 -0.394089
|
| Upper 0.95
| drugClozapine 0.271129
| drugTypical 2.024130Ld <- Ldiff( fit2, "drug", ref = "Atypical")
wald( fit2, Ld ) # compare Clozapine with Typical| numDF denDF F.value p.value
| 1 2 263 7.174332 0.00093
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.782102 0.834395 263 -2.135802 0.03362
| Typical - Atypical 0.896756 0.627306 263 1.429535 0.15404
| Typical - Clozapine 2.678858 0.715431 263 3.744397 0.00022
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -3.425047 -0.139158
| Typical - Atypical -0.338425 2.131937
| Typical - Clozapine 1.270156 4.087560#
# QUESTIONS:
#
# 1) Would it ever be appropriate to include a contextual variable even
# though it isn't significant? e.g. In this case if the p-value were
# 0.06, what would be the consequence of including it or excluding it?
#
# 2) Compare the SEs for the raw variables and for the contextual variable.
# What could explain the difference?
#
# === EXERCISE: ===
#
# 1) Estimate the between-patient differences (compositional effects)
# in these three drugs. Note that there are at least two ways of
# doing this.
#
# 2) Carry out a similar analysis for general symptoms: 'gen'
#
#
# Taking time into account
#
# When the range of time is compact and similar for all subjects
# and when time is not expected to have a very different effect
# as time progresses, simple models for time are generally adequate
#
# In the next session, we will see much richer functions of
# time: splines, asymptotic models, Fourier analysis, etc.
#
fit2l <- update(fit2, . ~ . + year)
summary( fit2l )| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1810.39 1840.334 -897.1952
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 4.992041 3.325507
|
| Fixed effects: neg ~ drug + cvar(drug, id) + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.511006 2.1742271 262 8.053899 0.0000
| drugClozapine -1.153198 0.8362982 262 -1.378931 0.1691
| drugTypical -0.117037 0.6785156 262 -0.172489 0.8632
| cvar(drug, id)Clozapine 3.184237 2.8728859 50 1.108376 0.2730
| cvar(drug, id)Typical -1.602418 2.9026725 50 -0.552049 0.5834
| year -0.477814 0.1359168 262 -3.515489 0.0005
| Correlation:
| (Intr) drgClz drgTyp c(,i)C c(,i)T
| drugClozapine 0.047
| drugTypical -0.093 0.397
| cvar(drug, id)Clozapine -0.766 -0.291 -0.116
| cvar(drug, id)Typical -0.837 -0.093 -0.234 0.658
| year -0.219 -0.214 0.425 0.062 -0.099
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.50585771 -0.61471185 -0.09486059 0.56161475 3.26692262
|
| Number of Observations: 318
| Number of Groups: 53ww <-wald( fit2l )
wald( fit2l, 'cvar' )| numDF denDF F.value p.value
| cvar 2 50 2.061446 0.13795
|
| Coefficients Estimate Std.Error DF t-value p-value
| cvar(drug, id)Clozapine 3.184237 2.872886 50 1.108376 0.27300
| cvar(drug, id)Typical -1.602418 2.902672 50 -0.552049 0.58337
|
| Coefficients Lower 0.95 Upper 0.95
| cvar(drug, id)Clozapine -2.586125 8.954598
| cvar(drug, id)Typical -7.432608 4.227771wald( fit2l, 'drug' )| numDF denDF F.value p.value
| drug 4 50 1.224521 0.31232
|
| Coefficients Estimate Std.Error DF t-value p-value
| drugClozapine -1.153198 0.836298 262 -1.378931 0.16909
| drugTypical -0.117037 0.678516 262 -0.172489 0.86319
| cvar(drug, id)Clozapine 3.184237 2.872886 50 1.108376 0.27300
| cvar(drug, id)Typical -1.602418 2.902672 50 -0.552049 0.58337
|
| Coefficients Lower 0.95 Upper 0.95
| drugClozapine -2.799919 0.493524
| drugTypical -1.453074 1.219001
| cvar(drug, id)Clozapine -2.586125 8.954598
| cvar(drug, id)Typical -7.432608 4.227771Ld <- Ldiff( fit2l, "drug", ref = "Atypical")
wald ( fit2l, Ld ) # What does this say about Clozapine ?| numDF denDF F.value p.value
| 1 2 262 1.034312 0.35692
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.153198 0.836298 262 -1.378931 0.16909
| Typical - Atypical -0.117037 0.678516 262 -0.172489 0.86319
| Typical - Clozapine 1.036161 0.842019 262 1.230567 0.21959
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.799919 0.493524
| Typical - Atypical -1.453074 1.219001
| Typical - Clozapine -0.621825 2.694147# Q: How do you explain the differences in the estimation of the Typical - Clozapine
# comparison in the 4 analyses:
#
lapply(
list("pooled" = fit.lm, "no ctx" = fit,"ctx"= fit2,"ctx+year" = fit2l),
function( fit ) wald( fit , Ldiff( fit,'drug', ref = "Atypical"))
)| $pooled
| numDF denDF F.value p.value
| 1 2 315 0.9501109 0.3878
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical 1.021073 0.967366 315 1.055519 0.29200
| Typical - Atypical -0.045073 0.862500 315 -0.052259 0.95836
| Typical - Clozapine -1.066146 0.809706 315 -1.316707 0.18889
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -0.882243 2.924388
| Typical - Atypical -1.742063 1.651917
| Typical - Clozapine -2.659262 0.526970
|
|
| $`no ctx`
| numDF denDF F.value p.value
| 1 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.300272 0.798060 263 -1.629290 0.10445
| Typical - Atypical 0.815021 0.614065 263 1.327255 0.18558
| Typical - Clozapine 2.115292 0.682395 263 3.099806 0.00215
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.871672 0.271129
| Typical - Atypical -0.394089 2.024130
| Typical - Clozapine 0.771639 3.458945
|
|
| $ctx
| numDF denDF F.value p.value
| 1 2 263 7.174332 0.00093
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.782102 0.834395 263 -2.135802 0.03362
| Typical - Atypical 0.896756 0.627306 263 1.429535 0.15404
| Typical - Clozapine 2.678858 0.715431 263 3.744397 0.00022
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -3.425047 -0.139158
| Typical - Atypical -0.338425 2.131937
| Typical - Clozapine 1.270156 4.087560
|
|
| $`ctx+year`
| numDF denDF F.value p.value
| 1 2 262 1.034312 0.35692
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.153198 0.836298 262 -1.378931 0.16909
| Typical - Atypical -0.117037 0.678516 262 -0.172489 0.86319
| Typical - Clozapine 1.036161 0.842019 262 1.230567 0.21959
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.799919 0.493524
| Typical - Atypical -1.453074 1.219001
| Typical - Clozapine -0.621825 2.694147clist <- lapply(
list("pooled" = fit.lm, "no ctx" = fit,"ctx"= fit2,"ctx+year" = fit2l),
function( fit ) wald( fit , Ldiff( fit,'drug', ref = "Atypical"))
)
clist| $pooled
| numDF denDF F.value p.value
| 1 2 315 0.9501109 0.3878
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical 1.021073 0.967366 315 1.055519 0.29200
| Typical - Atypical -0.045073 0.862500 315 -0.052259 0.95836
| Typical - Clozapine -1.066146 0.809706 315 -1.316707 0.18889
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -0.882243 2.924388
| Typical - Atypical -1.742063 1.651917
| Typical - Clozapine -2.659262 0.526970
|
|
| $`no ctx`
| numDF denDF F.value p.value
| 1 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.300272 0.798060 263 -1.629290 0.10445
| Typical - Atypical 0.815021 0.614065 263 1.327255 0.18558
| Typical - Clozapine 2.115292 0.682395 263 3.099806 0.00215
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.871672 0.271129
| Typical - Atypical -0.394089 2.024130
| Typical - Clozapine 0.771639 3.458945
|
|
| $ctx
| numDF denDF F.value p.value
| 1 2 263 7.174332 0.00093
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.782102 0.834395 263 -2.135802 0.03362
| Typical - Atypical 0.896756 0.627306 263 1.429535 0.15404
| Typical - Clozapine 2.678858 0.715431 263 3.744397 0.00022
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -3.425047 -0.139158
| Typical - Atypical -0.338425 2.131937
| Typical - Clozapine 1.270156 4.087560
|
|
| $`ctx+year`
| numDF denDF F.value p.value
| 1 2 262 1.034312 0.35692
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.153198 0.836298 262 -1.378931 0.16909
| Typical - Atypical -0.117037 0.678516 262 -0.172489 0.86319
| Typical - Clozapine 1.036161 0.842019 262 1.230567 0.21959
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.799919 0.493524
| Typical - Atypical -1.453074 1.219001
| Typical - Clozapine -0.621825 2.694147do.call( rbind, lapply(clist, function(x) x[[1]][[2]][3,]))|
| Estimate Std.Error DF t-value p-value Lower 0.95
| pooled -1.066146 0.8097063 315 -1.316707 0.1888939025 -2.6592623
| no ctx 2.115292 0.6823951 263 3.099806 0.0021468358 0.7716393
| ctx 2.678858 0.7154312 263 3.744397 0.0002221581 1.2701563
| ctx+year 1.036161 0.8420192 262 1.230567 0.2195885901 -0.6218253
|
| Upper 0.95
| pooled 0.526970
| no ctx 3.458945
| ctx 4.087560
| ctx+year 2.694147#
# Can you think of explanations consistent with these results?
#
#
# EXERCISE:
# Is there evidence that there is curvature in the effect of time?
#
# Level 1 and Level 2 diagnostics
#
#
# Level 1
#
plot( fit2l )
plot( fit2l, id = .02)
plot( fit2l, resid(.) ~ fitted(.) | drug, id = .02)
plot( fit2l, resid(.) ~ fitted(.) | drug * Sex, id = .02)
#
# Diagnostic for heteroskedasticity
#
plot( fit2l, sqrt( abs( resid(.))) ~ fitted(.), id = .02) # exploratory version
# fancy version
plot( fit2l, sqrt( abs( resid(.))) ~ fitted(.), id = .02, sub = "scale-location plot",
ylab = expression( sqrt(abs(" resid "))),
panel = function(x, y, ...){
panel.xyplot( x, y, ...)
panel.loess(x ,y, ...)
})
# id didn't work so we do it manually
#
if(interactive) {
trellis.focus() # trellis = lattice ?!? (the archaeology should be less visible)
panel.identify(labels = dd$id)
trellis.unfocus()
}
# no strong pattern here so we wait until later to see what to do about it.
qqnorm( fit2l, id = .05 )
qqnorm( fit2l, ~ resid(.) | drug, id = .05 )
# We note 'M47', 'F36', 'F41' as suspicious cases
xyplot( neg ~ year | id, dd, groups = drug, auto.key = T) # not clear for lecture
show.settings()
gd( pch = 15:17, col = c('blue','green4','red'))
show.settings()
xyplot( neg ~ year | id, dd, groups = drug, auto.key = T) # look for M47, F36 and F41
#
# Variogram: diagnostics for autocorrelation
# New level 1 diagnostic for longitudinal data
#
vv <- Variogram( fit2l , form = ~ year | id, maxdist = .7)
vv # variance of differences between pairs as a function of distance in time| variog dist n.pairs
| 1 0.8641087 1 265
| 2 0.9580987 2 212
| 3 1.0372402 3 159
| 4 1.1749621 4 106
| 5 1.2156057 5 53str(vv)| Classes 'Variogram' and 'data.frame': 5 obs. of 3 variables:
| $ variog : num 0.864 0.958 1.037 1.175 1.216
| $ dist : num 1 2 3 4 5
| $ n.pairs: int 265 212 159 106 53
| - attr(*, "collapse")= logi TRUE# plot(vv) # produceds an error
xyplot(variog ~ dist, vv, pch = 16, cex =2)
# Shows that differences between pairs of residuals that are close have smaller
# variance than between those that are far apart
#
# Level 2 diagnostics
#
ranef(fit2l) # Level 2 residuals BLUPS for u_is| (Intercept)
| M1 -7.8659136
| M5 -6.4547653
| M3 -5.2322173
| F2 -5.0770287
| F4 -5.1886799
| M21 -5.5236336
| F11 -5.5917474
| F36 -6.2753069
| M7 -3.6803311
| M6 -5.0145303
| F10 -3.3699539
| M22 -3.3264165
| M14 -5.3877126
| M8 -3.0160392
| F25 -3.3945303
| M12 -3.6178327
| M13 -1.6628790
| M15 -1.6628790
| M16 -1.5076904
| F24 -3.0986038
| F32 -0.8869359
| M17 -2.0510039
| M19 -3.6419172
| M37 -2.0122834
| M29 -3.4480081
| F41 -2.3955914
| M26 -2.9824423
| F23 -1.7748369
| M30 -0.4555803
| M20 1.1740535
| M35 1.6396194
| M28 1.2611283
| F27 1.0378258
| M40 0.8294660
| M31 1.8287116
| M33 2.6142885
| M34 1.3772898
| M46 2.9633862
| M43 4.7433918
| M44 4.6317405
| M47 1.9835936
| M39 5.2524950
| M38 6.0965519
| F48 2.9147253
| M49 3.2251025
| M42 6.7608438
| M52 5.5919589
| M50 6.3142389
| F45 6.3727189
| F54 5.8633090
| M51 8.2349824
| F55 10.8731889
| M53 12.0126804class(ranef(fit2l))| [1] "ranef.lme" "data.frame"methods(class = 'ranef.lme')| [1] plot print
| see '?methods' for accessing help and source codeplot( ranef( fit2l, aug = T ), form = ~ Sex) 
# note ranef( fit2l, aug = T) has a 'bug' and doesn't work if there's a matrix in data frame
# To plot BLUPS
RE <- cbind( ranef( fit2l), up( dd, ~ id))
head( RE )| (Intercept) age1 yrsill1 status Sex Subj id
| M1 -7.865914 26.4 8 Single M M1 M1
| M5 -6.454765 31.8 13 Single M M5 M5
| M3 -5.232217 32.6 14 Single M M3 M3
| F2 -5.077029 29.9 10 Single F F2 F2
| F4 -5.188680 30.2 8 Single F F4 F4
| M21 -5.523634 31.9 12 Single M M21 M21xyplot( `(Intercept)` ~ age1 | Sex , RE)
qqnorm( RE[[1]])
#
# Acting on diagnostics
# Should we drop some observations
#
fit2l.d <- update( fit2l, data = subset(dd, !(id %in% c("M47","F36","F41")) ))
# or
fit2l.d <- update( fit2l, subset = !(id %in% c("M47","F36","F41")) )
summary (fit2l.d)| Linear mixed-effects model fit by REML
| Data: dd
| Subset: !(id %in% c("M47", "F36", "F41"))
| AIC BIC logLik
| 1669.854 1699.323 -826.9271
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 5.064185 3.075854
|
| Fixed effects: neg ~ drug + cvar(drug, id) + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.337154 2.1910795 247 7.912608 0.0000
| drugClozapine -1.402456 0.8165726 247 -1.717491 0.0871
| drugTypical 0.358380 0.6441458 247 0.556365 0.5785
| cvar(drug, id)Clozapine 3.708741 2.9771527 47 1.245734 0.2190
| cvar(drug, id)Typical -2.286751 2.9236816 47 -0.782148 0.4381
| year -0.375001 0.1289879 247 -2.907255 0.0040
| Correlation:
| (Intr) drgClz drgTyp c(,i)C c(,i)T
| drugClozapine 0.039
| drugTypical -0.090 0.406
| cvar(drug, id)Clozapine -0.738 -0.274 -0.111
| cvar(drug, id)Typical -0.844 -0.089 -0.220 0.629
| year -0.206 -0.189 0.436 0.052 -0.096
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -2.36776534 -0.63680435 -0.09251977 0.55370236 2.79884975
|
| Number of Observations: 300
| Number of Groups: 50#
# Note: a problem with this approach: a cluster or observation could have high
# influence even though it does not have a large residual. So we could miss
# some influential points.
#
# There was a package 'influence.ME' that was being developed to produce
# influence diagnostics for 'lme' object but the version eventually released
# on CRAN only work for 'lmer' objects.
#
#
#
# EXERCISE:
# Compare results with fit2l.d and fit2l. Any substantial differences?
#
# In the following analysis, we keep the original data but this does not
# imply that this is the best course. If some outliers prove to have
# uncorrectable measurement errors, or if they are not in scope (some other
# diagnosis) then it would be appropriate to drop them.
#
#
#
# Incorporation of possible autocorrelation
#
#
# Add a correlation structure to the R side.
# Most common:
# corAR1: autoregressive of order 1 for evenly spaced data
# corCAR1: continuous autoregressive of order 1 for
# unevenly spaced data
# corARMA: autoregressive moving average of variable order
# -- other classes are mainly for spatial correlation
# -- you can write your own but it's challenging
# See: ?corClasses for a complete list
#
# Auto-regressive process of order 1
fit2lc <- update( fit2l, correlation = corAR1( form = ~ year | id))
summary( fit2lc)| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1806.119 1839.806 -894.0595
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 4.863878 3.477444
|
| Correlation Structure: AR(1)
| Formula: ~year | id
| Parameter estimate(s):
| Phi
| 0.2045159
| Fixed effects: neg ~ drug + cvar(drug, id) + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.551586 2.1753872 262 8.068258 0.0000
| drugClozapine -1.198927 0.8940467 262 -1.341011 0.1811
| drugTypical -0.241801 0.6875520 262 -0.351684 0.7254
| cvar(drug, id)Clozapine 3.164622 2.8743291 50 1.100995 0.2762
| cvar(drug, id)Typical -1.494575 2.8857604 50 -0.517914 0.6068
| year -0.485736 0.1503674 262 -3.230327 0.0014
| Correlation:
| (Intr) drgClz drgTyp c(,i)C c(,i)T
| drugClozapine 0.050
| drugTypical -0.102 0.375
| cvar(drug, id)Clozapine -0.758 -0.309 -0.115
| cvar(drug, id)Typical -0.831 -0.088 -0.233 0.654
| year -0.244 -0.211 0.403 0.066 -0.094
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.2515034 -0.6100687 -0.1129732 0.5311779 3.0912181
|
| Number of Observations: 318
| Number of Groups: 53# do we need the autocorrelation?
intervals( fit2lc ) # look at CI for Phi (correlation)| Approximate 95% confidence intervals
|
| Fixed effects:
| lower est. upper
| (Intercept) 13.2681190 17.5515864 21.8350537
| drugClozapine -2.9593580 -1.1989267 0.5615045
| drugTypical -1.5956322 -0.2418011 1.1120299
| cvar(drug, id)Clozapine -2.6086379 3.1646220 8.9378819
| cvar(drug, id)Typical -7.2907957 -1.4945753 4.3016451
| year -0.7818181 -0.4857358 -0.1896535
| attr(,"label")
| [1] "Fixed effects:"
|
| Random Effects:
| Level: id
| lower est. upper
| sd((Intercept)) 3.897772 4.863878 6.069445
|
| Correlation structure:
| lower est. upper
| Phi 0.03372634 0.2045159 0.3636998
| attr(,"label")
| [1] "Correlation structure:"
|
| Within-group standard error:
| lower est. upper
| 3.125207 3.477444 3.869381# or
anova( fit2l, fit2lc ) # this test can be ok because there is no boundary at 0 for Phi| Model df AIC BIC logLik Test L.Ratio p-value
| fit2l 1 8 1810.390 1840.334 -897.1952
| fit2lc 2 9 1806.119 1839.806 -894.0595 1 vs 2 6.271422 0.0123# but there is when using corCAR1.
# plot( simulate( fit2l, nsim = 1000, m2 = fit2lc)) # out of luck with fancier model
wald( fit2l )| numDF denDF F.value p.value
| 6 50 84.67485 <.00001
|
| Coefficients Estimate Std.Error DF t-value p-value
| (Intercept) 17.511006 2.174227 262 8.053899 <.00001
| drugClozapine -1.153198 0.836298 262 -1.378931 0.16909
| drugTypical -0.117037 0.678516 262 -0.172489 0.86319
| cvar(drug, id)Clozapine 3.184237 2.872886 50 1.108376 0.27300
| cvar(drug, id)Typical -1.602418 2.902672 50 -0.552049 0.58337
| year -0.477814 0.135917 262 -3.515489 0.00052
|
| Coefficients Lower 0.95 Upper 0.95
| (Intercept) 13.229823 21.792190
| drugClozapine -2.799919 0.493524
| drugTypical -1.453074 1.219001
| cvar(drug, id)Clozapine -2.586125 8.954598
| cvar(drug, id)Typical -7.432608 4.227771
| year -0.745442 -0.210186wald( fit2lc )| numDF denDF F.value p.value
| 6 50 84.33179 <.00001
|
| Coefficients Estimate Std.Error DF t-value p-value
| (Intercept) 17.551586 2.175387 262 8.068258 <.00001
| drugClozapine -1.198927 0.894047 262 -1.341011 0.18108
| drugTypical -0.241801 0.687552 262 -0.351684 0.72536
| cvar(drug, id)Clozapine 3.164622 2.874329 50 1.100995 0.27617
| cvar(drug, id)Typical -1.494575 2.885760 50 -0.517914 0.60680
| year -0.485736 0.150367 262 -3.230327 0.00139
|
| Coefficients Lower 0.95 Upper 0.95
| (Intercept) 13.268119 21.835054
| drugClozapine -2.959358 0.561505
| drugTypical -1.595632 1.112030
| cvar(drug, id)Clozapine -2.608638 8.937882
| cvar(drug, id)Typical -7.290796 4.301645
| year -0.781818 -0.189654Ld| [,1] [,2] [,3] [,4] [,5] [,6]
| Clozapine - Atypical 0 1 0 0 0 0
| Typical - Atypical 0 0 1 0 0 0
| Typical - Clozapine 0 -1 1 0 0 0wald( fit2l, Ld )| numDF denDF F.value p.value
| 1 2 262 1.034312 0.35692
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.153198 0.836298 262 -1.378931 0.16909
| Typical - Atypical -0.117037 0.678516 262 -0.172489 0.86319
| Typical - Clozapine 1.036161 0.842019 262 1.230567 0.21959
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.799919 0.493524
| Typical - Atypical -1.453074 1.219001
| Typical - Clozapine -0.621825 2.694147wald( fit2lc, Ld )| numDF denDF F.value p.value
| 1 2 262 0.9123941 0.40283
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.198927 0.894047 262 -1.341011 0.18108
| Typical - Atypical -0.241801 0.687552 262 -0.351684 0.72536
| Typical - Clozapine 0.957126 0.900737 262 1.062602 0.28894
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.959358 0.561505
| Typical - Atypical -1.595632 1.112030
| Typical - Clozapine -0.816480 2.730731# QUESTION: When would you expect results to be affected by
# including AR in the model?
#
# Note: AR (with positive autocorrelation) implies we expect obs.
# close in time to be closer than expected under the assumption
# of independence and observations that are far in time
# to be farther.
#
# Therefore a treatment difference between
# adjoining times gets more weight
# than one between distant times.
#
# That explains why including autocorrelation can change estimated
# differences and p-values.
#
#
# EXERCISE:
# Look at diagnostics for models predicting 'gen' instead of 'neg'
#
#
#
# Heteroscedasticiy revisited
#
fitg <- lme( gen ~ drug + cvar(drug,id) , dd, random = ~ 1 | id)
summary( fitg )| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 2161.396 2187.62 -1073.698
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 5.771655 6.277779
|
| Fixed effects: gen ~ drug + cvar(drug, id)
| Value Std.Error DF t-value p-value
| (Intercept) 30.074931 2.589763 263 11.613004 0.0000
| drugClozapine -2.900413 1.542192 263 -1.880709 0.0611
| drugTypical 4.278485 1.159434 263 3.690148 0.0003
| cvar(drug, id)Clozapine 9.931334 3.692515 50 2.689585 0.0097
| cvar(drug, id)Typical -3.457844 3.634982 50 -0.951268 0.3460
| Correlation:
| (Intr) drgClz drgTyp c(,i)C
| drugClozapine 0.000
| drugTypical 0.000 0.552
| cvar(drug, id)Clozapine -0.732 -0.418 -0.231
| cvar(drug, id)Typical -0.858 -0.176 -0.319 0.658
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -2.18006804 -0.62712129 -0.08655213 0.52664934 4.25622647
|
| Number of Observations: 318
| Number of Groups: 53fitgl <- lme( gen ~ drug + cvar(drug,id) + year, dd, random = ~ 1 | id)
summary( fitgl )| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 2148.02 2177.964 -1066.01
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 5.803759 6.097606
|
| Fixed effects: gen ~ drug + cvar(drug, id) + year
| Value Std.Error DF t-value p-value
| (Intercept) 33.64712 2.732709 262 12.312733 0.0000
| drugClozapine -1.55705 1.533425 262 -1.015409 0.3108
| drugTypical 2.11299 1.244117 262 1.698387 0.0906
| cvar(drug, id)Clozapine 8.58797 3.688862 50 2.328082 0.0240
| cvar(drug, id)Typical -1.29235 3.662873 50 -0.352825 0.7257
| year -1.02063 0.249215 262 -4.095356 0.0001
| Correlation:
| (Intr) drgClz drgTyp c(,i)C c(,i)T
| drugClozapine 0.068
| drugTypical -0.136 0.397
| cvar(drug, id)Clozapine -0.723 -0.416 -0.165
| cvar(drug, id)Typical -0.761 -0.135 -0.340 0.636
| year -0.319 -0.214 0.425 0.089 -0.144
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -2.3078278 -0.5731915 -0.1464905 0.5223173 4.2237955
|
| Number of Observations: 318
| Number of Groups: 53wald( fitgl, Ld)| numDF denDF F.value p.value
| 1 2 262 3.137857 0.04501
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.557054 1.533425 262 -1.015409 0.31085
| Typical - Atypical 2.112993 1.244117 262 1.698387 0.09062
| Typical - Clozapine 3.670046 1.543915 262 2.377103 0.01817
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -4.576460 1.462353
| Typical - Atypical -0.336749 4.562734
| Typical - Clozapine 0.629985 6.710108# ... many diagnostics later ....
# Note: I believe there is an error in documentation and the default
# for resid is type = 'r' for 'raw' or 'response'. We prefer the
# standardized, 'pearson' residuals here.
plot( fitgl, sqrt( abs( resid(., type = 'p'))) ~ fitted(.), id = .03)
plot( fitgl, sqrt( abs( resid(., type = 'p'))) ~ fitted(.), id = .03,
panel = function(x, y, ...) {
panel.xyplot( x, y, ...)
panel.loess( x, y,...)
})
plot( fitgl, sqrt( abs( resid(., type = 'p'))) ~ fitted(.) | Sex, id = .03,
panel = function(x, y, ...) {
panel.xyplot( x, y, ...)
panel.loess( x, y,...)
})
plot( fitgl, sqrt( abs( resid(., type = 'p'))) ~ fitted(.) | drug, id = .03,
panel = function(x, y, ...) {
panel.xyplot( x, y, ...)
panel.loess( x, y,...)
})
#
# Perhaps: increased variability with increase in predicted value
# In practice, I wouldn't be too disturbed by this plot but
# for pedagogical purposes, let's see what we could do to address this.
#
#
# Accounting for heteroskedasticity
#
#
# We assume that the SD of residuals changes with the fitted value (variance covariate)
#
# Common assumptions are that:
#
# 1) variance is proportional to an unknown power of the variance covariate
#
# or
#
# 2) a unknown constant plus an unknown power of the covariate
#
# To see the full set of methods, available (you can also create your own):
?varClasses| starting httpd help server ... done# Here we illustrate varConstPower
fitglh <- update( fitgl, weights = varConstPower( form = ~fitted(.)| drug) )
summary(fitglh)| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 2131.661 2184.063 -1051.83
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 5.813634 0.03365691
|
| Variance function:
| Structure: Constant plus power of variance covariate, different strata
| Formula: ~fitted(.) | drug
| Parameter estimates:
| Typical Atypical Clozapine
| const 0.000895613 0.001125126 58.502605
| power 1.513026811 1.489908517 1.310149
| Fixed effects: gen ~ drug + cvar(drug, id) + year
| Value Std.Error DF t-value p-value
| (Intercept) 32.80306 2.662644 262 12.319732 0.0000
| drugClozapine -2.13808 1.508406 262 -1.417441 0.1575
| drugTypical 2.37355 1.142420 262 2.077649 0.0387
| cvar(drug, id)Clozapine 9.14154 3.712537 50 2.462343 0.0173
| cvar(drug, id)Typical -1.66796 3.581094 50 -0.465769 0.6434
| year -0.85033 0.215380 262 -3.948066 0.0001
| Correlation:
| (Intr) drgClz drgTyp c(,i)C c(,i)T
| drugClozapine 0.064
| drugTypical -0.091 0.349
| cvar(drug, id)Clozapine -0.709 -0.457 -0.163
| cvar(drug, id)Typical -0.785 -0.128 -0.328 0.625
| year -0.299 -0.174 0.420 0.081 -0.132
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -2.10072756 -0.55449512 -0.09397095 0.58046097 4.45778160
|
| Number of Observations: 318
| Number of Groups: 53anova( fitgl, fitglh )| Model df AIC BIC logLik Test L.Ratio p-value
| fitgl 1 8 2148.020 2177.964 -1066.01
| fitglh 2 14 2131.661 2184.063 -1051.83 1 vs 2 28.3594 1e-04plot( fitglh, sqrt( abs( resid(.,type = 'p'))) ~ fitted(.) | drug, id = .03,
panel = function(x, y, ...) {
panel.xyplot( x, y, ...)
panel.loess( x, y,...)
})
# note that 'pearson residuals' are plotted so that
# high variance residuals will be shrunk.
# with raw residuals we expect heteroskedasticity to look
# even worse because the high variance observations have
# less influence on the fit
# raw residuals:
plot( fitglh, sqrt( abs( resid(., type = 'r'))) ~ fitted(.) | drug, id = .03,
panel = function(x, y, ...) {
panel.xyplot( x, y, ...)
panel.loess( x, y,...)
})
plot( fitgl, sqrt( abs( resid(., type = 'r'))) ~ fitted(.) | drug, id = .03,
panel = function(x, y, ...) {
panel.xyplot( x, y, ...)
panel.loess( x, y,...)
})
wald( fitgl, Ld)| numDF denDF F.value p.value
| 1 2 262 3.137857 0.04501
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.557054 1.533425 262 -1.015409 0.31085
| Typical - Atypical 2.112993 1.244117 262 1.698387 0.09062
| Typical - Clozapine 3.670046 1.543915 262 2.377103 0.01817
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -4.576460 1.462353
| Typical - Atypical -0.336749 4.562734
| Typical - Clozapine 0.629985 6.710108wald( fitglh, Ld)| numDF denDF F.value p.value
| 1 2 262 4.774854 0.00919
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -2.138077 1.508406 262 -1.417441 0.15754
| Typical - Atypical 2.373549 1.142420 262 2.077649 0.03872
| Typical - Clozapine 4.511625 1.541482 262 2.926810 0.00373
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -5.108219 0.832065
| Typical - Atypical 0.124055 4.623042
| Typical - Clozapine 1.476355 7.546896#
# Note that the impact of accounting for heteroskedasticity is not negligible.
#
# As usual, seeking higher efficiency has an impact on what we're estimating.
# Differences between drugs in patients at lower levels of 'gen' receive
# more weight than those at higher levels. The estimate has lower variance
# and is more precise but it might estimate something slightly different if
# there are differences in drug effects at different levels of severity.
# If this is an important question to address we could reformulate the model
# to attempt to take it into account although we might not have much power
# to detect such an effect.
#
# EXERCISE:
# Refit 'fitglh' and 'fitgl' dropping a few residual outliers.
# Retest the difference between the two models.
# Does accounting for heteroskedastiticy still improve the fit?
# Might it be that heteroskedasticity was just accounting for a
# few outliers instead of capturing a general phenomenon?
#
# EXERCISE:
# Enlarge the RE model. How far can you go and what is the impact
# of doing so.
#
# Visualizing different models
#
#
# The following code was used to generate graphs for a lecture
# It shows how the different models above fit the data and
# attemps to explain how and why they differ.
#
gd( pch = 15:17, col = c('blue','green4','red'))
xyplot( neg ~ year | Subj, dd, groups = drug )
dd$id <- dd$Subj
dd$id <- reorder( dd$id, 1000*(dd$Sex == "M") + dd$neg )
gd( pch = 15:17, col = c('blue','green4','red'))
xyplot( neg ~ year | id, dd, groups = drug,
between = list( y = c(0,.5,0,0,0,0,0)) ,
skip = c( rep(F,15), T, rep(F,30)),
layout = c( 8,7),
auto.key = list( columns = 3))
# select a 'representative' sample of 3 from each sex
sel <- c( "F4","F25","F45","M12","M40","M51")
gd( pch = 15:17, col = c('blue','green4','red'), cex = 1.2)
xyplot( neg ~ year | id, dd, groups = drug,
between = list( y = c(0,.5,0,0,0,0,0)) ,
skip = c( rep(F,15), T, rep(F,30)),
layout = c( 3,2),
auto.key = list( columns = 3),
subset = id %in% sel)
fit.ols <- lm( neg ~ drug, dd )
summary( fit.ols )|
| Call:
| lm(formula = neg ~ drug, data = dd)
|
| Residuals:
| Min 1Q Median 3Q Max
| -9.299 -5.233 -1.233 3.767 18.767
|
| Coefficients:
| Estimate Std. Error t value Pr(>|t|)
| (Intercept) 15.27778 0.71557 21.351 <2e-16 ***
| drugClozapine 1.02107 0.96737 1.056 0.292
| drugTypical -0.04507 0.86250 -0.052 0.958
| ---
| Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
|
| Residual standard error: 6.072 on 315 degrees of freedom
| Multiple R-squared: 0.005996, Adjusted R-squared: -0.0003149
| F-statistic: 0.9501 on 2 and 315 DF, p-value: 0.3878L <- list(
'predicted' = cbind( 1, contrasts( dd$drug)),
'differences' = Ldiff ( fit.ols, 'drug', ref = 'Atypical'))
L| $predicted
| Clozapine Typical
| Atypical 1 0 0
| Clozapine 1 1 0
| Typical 1 0 1
|
| $differences
| [,1] [,2] [,3]
| Clozapine - Atypical 0 1 0
| Typical - Atypical 0 0 1
| Typical - Clozapine 0 -1 1wald( fit.ols, L)| numDF denDF F.value p.value
| predicted 3 315 694.4892 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| Atypical 15.27778 0.715570 315 21.35051 <.00001 13.86988 16.68568
| Clozapine 16.29885 0.650966 315 25.03793 <.00001 15.01806 17.57964
| Typical 15.23270 0.481526 315 31.63424 <.00001 14.28529 16.18012
|
| numDF denDF F.value p.value
| differences 2 315 0.9501109 0.3878
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical 1.021073 0.967366 315 1.055519 0.29200
| Typical - Atypical -0.045073 0.862500 315 -0.052259 0.95836
| Typical - Clozapine -1.066146 0.809706 315 -1.316707 0.18889
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -0.882243 2.924388
| Typical - Atypical -1.742063 1.651917
| Typical - Clozapine -2.659262 0.526970# Clozapine is worst but not significantly so
#
# Displaying data and fitted values
#
pred <- expand.grid( year = 1:6, drug = levels( dd$drug), id = levels(dd$id))
head( pred )| year drug id
| 1 1 Atypical F2
| 2 2 Atypical F2
| 3 3 Atypical F2
| 4 4 Atypical F2
| 5 5 Atypical F2
| 6 6 Atypical F2some( pred )| year drug id
| 38 2 Atypical F11
| 45 3 Clozapine F11
| 132 6 Atypical F32
| 336 6 Clozapine M21
| 365 5 Atypical M6
| 392 2 Typical M22
| 449 5 Typical M12
| 469 1 Atypical M15
| 475 1 Clozapine M15
| 735 3 Typical M34pred$neg <- predict( fit.ols, newdata = pred )
# Combine with data for plotting
dd$what <- factor("data")
pred$what <- factor("fit")
gd( pch = 15:17, col = c('blue','green4','red'), cex = 1.2,
lty = c(1,1,1,1,1,1,2,2,2), lwd = 2)
xyplot( neg ~ year | id, Rbind( dd, pred), groups = what:drug,
panel = panel.superpose.2,
type = c('p','p','p','l','l','l'),
subset = id %in% sel,
between = list( y =1 ),
auto.key = list( text = c("Atypical","Clozapine","Typical"),
points = T, lines = T, columns = 3))| [1] "age1" "yrsill1" "pos" "neg" "gen" "total" "year"
| [8] "drug" "status" "Sex" "Subj" "id" "what"
# overplotting problem
#
# make later lines less dense:
gd( pch = 15:17, col = c('blue','green4','red'), cex = 1.2,
lty = c(1,2,3,1,2,3,1,2,3),lwd = 2)
xyplot( neg ~ year | id, Rbind(dd, pred), groups = what:drug,
panel = panel.superpose.2,
type = c('p','p','p','l','l','l'),
subset = id %in% sel,
between = list( y =1 ),
auto.key = list( text = c("Atypical","Clozapine","Typical"),
points = T, lines = T, columns = 3))| [1] "age1" "yrsill1" "pos" "neg" "gen" "total" "year"
| [8] "drug" "status" "Sex" "Subj" "id" "what"
#
# Mixed model: first version
#
fit.mix <- lme( neg ~ drug, dd, random = ~ 1 | id )
summary(fit.mix)| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1829.011 1847.773 -909.5053
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 5.240588 3.399672
|
| Fixed effects: neg ~ drug
| Value Std.Error DF t-value p-value
| (Intercept) 15.482816 0.8783362 263 17.627437 0.0000
| drugClozapine -1.300272 0.7980602 263 -1.629290 0.1044
| drugTypical 0.815021 0.6140648 263 1.327255 0.1856
| Correlation:
| (Intr) drgClz
| drugClozapine -0.444
| drugTypical -0.489 0.559
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.71174345 -0.59155383 -0.05454849 0.47584480 3.47145299
|
| Number of Observations: 318
| Number of Groups: 53wald( fit.mix, L ) # same L as above because same FE model| numDF denDF F.value p.value
| predicted 3 263 148.3786 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| Atypical 15.48282 0.878336 263 17.62744 <.00001 13.75335 17.21228
| Clozapine 14.18254 0.886412 263 15.99995 <.00001 12.43718 15.92791
| Typical 16.29784 0.788330 263 20.67387 <.00001 14.74559 17.85008
|
| numDF denDF F.value p.value
| differences 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.300272 0.798060 263 -1.629290 0.10445
| Typical - Atypical 0.815021 0.614065 263 1.327255 0.18558
| Typical - Clozapine 2.115292 0.682395 263 3.099806 0.00215
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.871672 0.271129
| Typical - Atypical -0.394089 2.024130
| Typical - Clozapine 0.771639 3.458945# Clozapine is best and sig. diff from Typical
#
# Can use same pred for Popn pred and for BLUPS
#
pred$neg <- predict( fit.mix, pred, level = 0)
# Make a separate data frame for blupss
pred2 <- pred
pred2$neg <- predict( fit.mix, pred2, level = 1)
pred2$what <- factor("fit2")
# Cut and paste from last, then add 'pred2'
gd( pch = 15:17, col = c('blue','green4','red'), cex = 1.2,
lty = c(1,1,1,2,2,2,1,1,1),lwd = 2)
xyplot( neg ~ year | id, Rbind( dd, pred, pred2), groups = what:drug,
panel = panel.superpose.2,
lwd = 2,
type = c(rep('p',3),rep('l',6)),
subset = id %in% sel,
between = list( y =1 ),
auto.key = list( text = c("Atypical","Clozapine","Typical"),
points = T, lines = T, columns = 3))| [1] "age1" "yrsill1" "pos" "neg" "gen" "total" "year"
| [8] "drug" "status" "Sex" "Subj" "id" "what"
#
# Hausman test and mixed model with contextual categorical variable
#
dd$drug.m <- cvar( dd$drug, dd$id) # note that drug.m is matrix
dd$drug.m| Clozapine Typical
| [1,] 0.6666667 0.1666667
| [2,] 0.6666667 0.1666667
| [3,] 0.6666667 0.1666667
| [4,] 0.6666667 0.1666667
| [5,] 0.6666667 0.1666667
| [6,] 0.6666667 0.1666667
| [7,] 0.6666667 0.1666667
| [8,] 0.6666667 0.1666667
| [9,] 0.6666667 0.1666667
| [10,] 0.6666667 0.1666667
| [11,] 0.6666667 0.1666667
| [12,] 0.6666667 0.1666667
| [13,] 0.0000000 1.0000000
| [14,] 0.0000000 1.0000000
| [15,] 0.0000000 1.0000000
| [16,] 0.0000000 1.0000000
| [17,] 0.0000000 1.0000000
| [18,] 0.0000000 1.0000000
| [19,] 0.0000000 0.8333333
| [20,] 0.0000000 0.8333333
| [21,] 0.0000000 0.8333333
| [22,] 0.0000000 0.8333333
| [23,] 0.0000000 0.8333333
| [24,] 0.0000000 0.8333333
| [25,] 0.0000000 1.0000000
| [26,] 0.0000000 1.0000000
| [27,] 0.0000000 1.0000000
| [28,] 0.0000000 1.0000000
| [29,] 0.0000000 1.0000000
| [30,] 0.0000000 1.0000000
| [31,] 0.0000000 1.0000000
| [32,] 0.0000000 1.0000000
| [33,] 0.0000000 1.0000000
| [34,] 0.0000000 1.0000000
| [35,] 0.0000000 1.0000000
| [36,] 0.0000000 1.0000000
| [37,] 0.0000000 0.0000000
| [38,] 0.0000000 0.0000000
| [39,] 0.0000000 0.0000000
| [40,] 0.0000000 0.0000000
| [41,] 0.0000000 0.0000000
| [42,] 0.0000000 0.0000000
| [43,] 0.6666667 0.1666667
| [44,] 0.6666667 0.1666667
| [45,] 0.6666667 0.1666667
| [46,] 0.6666667 0.1666667
| [47,] 0.6666667 0.1666667
| [48,] 0.6666667 0.1666667
| [49,] 0.5000000 0.5000000
| [50,] 0.5000000 0.5000000
| [51,] 0.5000000 0.5000000
| [52,] 0.5000000 0.5000000
| [53,] 0.5000000 0.5000000
| [54,] 0.5000000 0.5000000
| [55,] 0.0000000 0.5000000
| [56,] 0.0000000 0.5000000
| [57,] 0.0000000 0.5000000
| [58,] 0.0000000 0.5000000
| [59,] 0.0000000 0.5000000
| [60,] 0.0000000 0.5000000
| [61,] 0.0000000 0.1666667
| [62,] 0.0000000 0.1666667
| [63,] 0.0000000 0.1666667
| [64,] 0.0000000 0.1666667
| [65,] 0.0000000 0.1666667
| [66,] 0.0000000 0.1666667
| [67,] 1.0000000 0.0000000
| [68,] 1.0000000 0.0000000
| [69,] 1.0000000 0.0000000
| [70,] 1.0000000 0.0000000
| [71,] 1.0000000 0.0000000
| [72,] 1.0000000 0.0000000
| [73,] 1.0000000 0.0000000
| [74,] 1.0000000 0.0000000
| [75,] 1.0000000 0.0000000
| [76,] 1.0000000 0.0000000
| [77,] 1.0000000 0.0000000
| [78,] 1.0000000 0.0000000
| [79,] 0.3333333 0.5000000
| [80,] 0.3333333 0.5000000
| [81,] 0.3333333 0.5000000
| [82,] 0.3333333 0.5000000
| [83,] 0.3333333 0.5000000
| [84,] 0.3333333 0.5000000
| [85,] 0.3333333 0.5000000
| [86,] 0.3333333 0.5000000
| [87,] 0.3333333 0.5000000
| [88,] 0.3333333 0.5000000
| [89,] 0.3333333 0.5000000
| [90,] 0.3333333 0.5000000
| [91,] 0.5000000 0.3333333
| [92,] 0.5000000 0.3333333
| [93,] 0.5000000 0.3333333
| [94,] 0.5000000 0.3333333
| [95,] 0.5000000 0.3333333
| [96,] 0.5000000 0.3333333
| [97,] 0.0000000 0.3333333
| [98,] 0.0000000 0.3333333
| [99,] 0.0000000 0.3333333
| [100,] 0.0000000 0.3333333
| [101,] 0.0000000 0.3333333
| [102,] 0.0000000 0.3333333
| [103,] 0.0000000 0.3333333
| [104,] 0.0000000 0.3333333
| [105,] 0.0000000 0.3333333
| [106,] 0.0000000 0.3333333
| [107,] 0.0000000 0.3333333
| [108,] 0.0000000 0.3333333
| [109,] 0.0000000 1.0000000
| [110,] 0.0000000 1.0000000
| [111,] 0.0000000 1.0000000
| [112,] 0.0000000 1.0000000
| [113,] 0.0000000 1.0000000
| [114,] 0.0000000 1.0000000
| [115,] 0.0000000 1.0000000
| [116,] 0.0000000 1.0000000
| [117,] 0.0000000 1.0000000
| [118,] 0.0000000 1.0000000
| [119,] 0.0000000 1.0000000
| [120,] 0.0000000 1.0000000
| [121,] 1.0000000 0.0000000
| [122,] 1.0000000 0.0000000
| [123,] 1.0000000 0.0000000
| [124,] 1.0000000 0.0000000
| [125,] 1.0000000 0.0000000
| [126,] 1.0000000 0.0000000
| [127,] 0.0000000 0.1666667
| [128,] 0.0000000 0.1666667
| [129,] 0.0000000 0.1666667
| [130,] 0.0000000 0.1666667
| [131,] 0.0000000 0.1666667
| [132,] 0.0000000 0.1666667
| [133,] 0.0000000 0.8333333
| [134,] 0.0000000 0.8333333
| [135,] 0.0000000 0.8333333
| [136,] 0.0000000 0.8333333
| [137,] 0.0000000 0.8333333
| [138,] 0.0000000 0.8333333
| [139,] 1.0000000 0.0000000
| [140,] 1.0000000 0.0000000
| [141,] 1.0000000 0.0000000
| [142,] 1.0000000 0.0000000
| [143,] 1.0000000 0.0000000
| [144,] 1.0000000 0.0000000
| [145,] 0.0000000 0.8333333
| [146,] 0.0000000 0.8333333
| [147,] 0.0000000 0.8333333
| [148,] 0.0000000 0.8333333
| [149,] 0.0000000 0.8333333
| [150,] 0.0000000 0.8333333
| [151,] 0.0000000 1.0000000
| [152,] 0.0000000 1.0000000
| [153,] 0.0000000 1.0000000
| [154,] 0.0000000 1.0000000
| [155,] 0.0000000 1.0000000
| [156,] 0.0000000 1.0000000
| [157,] 0.6666667 0.3333333
| [158,] 0.6666667 0.3333333
| [159,] 0.6666667 0.3333333
| [160,] 0.6666667 0.3333333
| [161,] 0.6666667 0.3333333
| [162,] 0.6666667 0.3333333
| [163,] 0.0000000 0.1666667
| [164,] 0.0000000 0.1666667
| [165,] 0.0000000 0.1666667
| [166,] 0.0000000 0.1666667
| [167,] 0.0000000 0.1666667
| [168,] 0.0000000 0.1666667
| [169,] 0.0000000 1.0000000
| [170,] 0.0000000 1.0000000
| [171,] 0.0000000 1.0000000
| [172,] 0.0000000 1.0000000
| [173,] 0.0000000 1.0000000
| [174,] 0.0000000 1.0000000
| [175,] 0.0000000 1.0000000
| [176,] 0.0000000 1.0000000
| [177,] 0.0000000 1.0000000
| [178,] 0.0000000 1.0000000
| [179,] 0.0000000 1.0000000
| [180,] 0.0000000 1.0000000
| [181,] 0.3333333 0.6666667
| [182,] 0.3333333 0.6666667
| [183,] 0.3333333 0.6666667
| [184,] 0.3333333 0.6666667
| [185,] 0.3333333 0.6666667
| [186,] 0.3333333 0.6666667
| [187,] 0.0000000 0.6666667
| [188,] 0.0000000 0.6666667
| [189,] 0.0000000 0.6666667
| [190,] 0.0000000 0.6666667
| [191,] 0.0000000 0.6666667
| [192,] 0.0000000 0.6666667
| [193,] 0.8333333 0.1666667
| [194,] 0.8333333 0.1666667
| [195,] 0.8333333 0.1666667
| [196,] 0.8333333 0.1666667
| [197,] 0.8333333 0.1666667
| [198,] 0.8333333 0.1666667
| [199,] 0.0000000 0.8333333
| [200,] 0.0000000 0.8333333
| [201,] 0.0000000 0.8333333
| [202,] 0.0000000 0.8333333
| [203,] 0.0000000 0.8333333
| [204,] 0.0000000 0.8333333
| [205,] 0.0000000 0.8333333
| [206,] 0.0000000 0.8333333
| [207,] 0.0000000 0.8333333
| [208,] 0.0000000 0.8333333
| [209,] 0.0000000 0.8333333
| [210,] 0.0000000 0.8333333
| [211,] 0.1666667 0.5000000
| [212,] 0.1666667 0.5000000
| [213,] 0.1666667 0.5000000
| [214,] 0.1666667 0.5000000
| [215,] 0.1666667 0.5000000
| [216,] 0.1666667 0.5000000
| [217,] 0.5000000 0.5000000
| [218,] 0.5000000 0.5000000
| [219,] 0.5000000 0.5000000
| [220,] 0.5000000 0.5000000
| [221,] 0.5000000 0.5000000
| [222,] 0.5000000 0.5000000
| [223,] 0.3333333 0.6666667
| [224,] 0.3333333 0.6666667
| [225,] 0.3333333 0.6666667
| [226,] 0.3333333 0.6666667
| [227,] 0.3333333 0.6666667
| [228,] 0.3333333 0.6666667
| [229,] 0.0000000 1.0000000
| [230,] 0.0000000 1.0000000
| [231,] 0.0000000 1.0000000
| [232,] 0.0000000 1.0000000
| [233,] 0.0000000 1.0000000
| [234,] 0.0000000 1.0000000
| [235,] 0.0000000 0.8333333
| [236,] 0.0000000 0.8333333
| [237,] 0.0000000 0.8333333
| [238,] 0.0000000 0.8333333
| [239,] 0.0000000 0.8333333
| [240,] 0.0000000 0.8333333
| [241,] 0.0000000 0.1666667
| [242,] 0.0000000 0.1666667
| [243,] 0.0000000 0.1666667
| [244,] 0.0000000 0.1666667
| [245,] 0.0000000 0.1666667
| [246,] 0.0000000 0.1666667
| [247,] 0.0000000 0.5000000
| [248,] 0.0000000 0.5000000
| [249,] 0.0000000 0.5000000
| [250,] 0.0000000 0.5000000
| [251,] 0.0000000 0.5000000
| [252,] 0.0000000 0.5000000
| [253,] 0.3333333 0.5000000
| [254,] 0.3333333 0.5000000
| [255,] 0.3333333 0.5000000
| [256,] 0.3333333 0.5000000
| [257,] 0.3333333 0.5000000
| [258,] 0.3333333 0.5000000
| [259,] 0.0000000 0.6666667
| [260,] 0.0000000 0.6666667
| [261,] 0.0000000 0.6666667
| [262,] 0.0000000 0.6666667
| [263,] 0.0000000 0.6666667
| [264,] 0.0000000 0.6666667
| [265,] 0.5000000 0.0000000
| [266,] 0.5000000 0.0000000
| [267,] 0.5000000 0.0000000
| [268,] 0.5000000 0.0000000
| [269,] 0.5000000 0.0000000
| [270,] 0.5000000 0.0000000
| [271,] 0.5000000 0.5000000
| [272,] 0.5000000 0.5000000
| [273,] 0.5000000 0.5000000
| [274,] 0.5000000 0.5000000
| [275,] 0.5000000 0.5000000
| [276,] 0.5000000 0.5000000
| [277,] 0.3333333 0.3333333
| [278,] 0.3333333 0.3333333
| [279,] 0.3333333 0.3333333
| [280,] 0.3333333 0.3333333
| [281,] 0.3333333 0.3333333
| [282,] 0.3333333 0.3333333
| [283,] 0.0000000 0.6666667
| [284,] 0.0000000 0.6666667
| [285,] 0.0000000 0.6666667
| [286,] 0.0000000 0.6666667
| [287,] 0.0000000 0.6666667
| [288,] 0.0000000 0.6666667
| [289,] 0.3333333 0.5000000
| [290,] 0.3333333 0.5000000
| [291,] 0.3333333 0.5000000
| [292,] 0.3333333 0.5000000
| [293,] 0.3333333 0.5000000
| [294,] 0.3333333 0.5000000
| [295,] 0.0000000 0.8333333
| [296,] 0.0000000 0.8333333
| [297,] 0.0000000 0.8333333
| [298,] 0.0000000 0.8333333
| [299,] 0.0000000 0.8333333
| [300,] 0.0000000 0.8333333
| [301,] 1.0000000 0.0000000
| [302,] 1.0000000 0.0000000
| [303,] 1.0000000 0.0000000
| [304,] 1.0000000 0.0000000
| [305,] 1.0000000 0.0000000
| [306,] 1.0000000 0.0000000
| [307,] 1.0000000 0.0000000
| [308,] 1.0000000 0.0000000
| [309,] 1.0000000 0.0000000
| [310,] 1.0000000 0.0000000
| [311,] 1.0000000 0.0000000
| [312,] 1.0000000 0.0000000
| [313,] 0.0000000 0.3333333
| [314,] 0.0000000 0.3333333
| [315,] 0.0000000 0.3333333
| [316,] 0.0000000 0.3333333
| [317,] 0.0000000 0.3333333
| [318,] 0.0000000 0.3333333head( dd, 18)| age1 yrsill1 pos neg gen total year drug status Sex
| 17 29.1 10 16 8 33 57 1 Typical Married F
| 72 29.1 10 21 23 60 104 2 Atypical Married F
| 127 29.1 10 9 7 20 36 3 Clozapine Married F
| 182 29.1 10 8 9 25 42 4 Clozapine Married F
| 237 29.1 10 8 7 25 40 5 Clozapine Married F
| 292 29.1 10 12 7 26 45 6 Clozapine Married F
| 38 36.4 16 12 8 29 49 1 Typical Separated/Divorced F
| 93 36.4 16 13 10 31 54 2 Atypical Separated/Divorced F
| 148 36.4 16 20 24 40 84 3 Clozapine Separated/Divorced F
| 203 36.4 16 15 16 35 66 4 Clozapine Separated/Divorced F
| 258 36.4 16 13 13 28 54 5 Clozapine Separated/Divorced F
| 313 36.4 16 11 15 22 48 6 Clozapine Separated/Divorced F
| 14 29.9 10 17 10 28 55 1 Typical Single F
| 69 29.9 10 14 10 25 49 2 Typical Single F
| 124 29.9 10 10 8 22 40 3 Typical Single F
| 179 29.9 10 15 10 28 53 4 Typical Single F
| 234 29.9 10 9 7 22 38 5 Typical Single F
| 289 29.9 10 11 7 20 38 6 Typical Single F
| Subj id what drug.m.Clozapine drug.m.Typical
| 17 F36 F36 data 0.6666667 0.1666667
| 72 F36 F36 data 0.6666667 0.1666667
| 127 F36 F36 data 0.6666667 0.1666667
| 182 F36 F36 data 0.6666667 0.1666667
| 237 F36 F36 data 0.6666667 0.1666667
| 292 F36 F36 data 0.6666667 0.1666667
| 38 F41 F41 data 0.6666667 0.1666667
| 93 F41 F41 data 0.6666667 0.1666667
| 148 F41 F41 data 0.6666667 0.1666667
| 203 F41 F41 data 0.6666667 0.1666667
| 258 F41 F41 data 0.6666667 0.1666667
| 313 F41 F41 data 0.6666667 0.1666667
| 14 F2 F2 data 0.0000000 1.0000000
| 69 F2 F2 data 0.0000000 1.0000000
| 124 F2 F2 data 0.0000000 1.0000000
| 179 F2 F2 data 0.0000000 1.0000000
| 234 F2 F2 data 0.0000000 1.0000000
| 289 F2 F2 data 0.0000000 1.0000000some( dd )| age1 yrsill1 pos neg gen total year drug status Sex Subj id what
| 59 36.6 20 16 23 40 79 2 Typical Single F F45 F45 data
| 139 32.6 14 23 7 26 56 3 Typical Single M M3 M3 data
| 231 34.0 12 12 9 23 44 5 Typical Single M M16 M16 data
| 319 28.5 11 17 9 25 51 6 Clozapine Single M M17 M17 data
| 171 32.8 10 13 19 29 61 4 Typical Single M M46 M46 data
| 186 37.3 16 14 19 30 63 4 Typical Single M M31 M31 data
| 12 25.0 8 18 21 52 91 1 Typical Single M M30 M30 data
| 67 25.0 8 12 13 29 54 2 Typical Single M M30 M30 data
| 251 32.7 14 22 17 41 80 5 Clozapine Single M M49 M49 data
| 277 38.8 24 10 29 24 63 6 Atypical Single M M53 M53 data
| drug.m.Clozapine drug.m.Typical
| 59 0.3333333 0.5000000
| 139 0.0000000 1.0000000
| 231 0.0000000 1.0000000
| 319 0.3333333 0.6666667
| 171 0.1666667 0.5000000
| 186 0.3333333 0.6666667
| 12 0.3333333 0.5000000
| 67 0.3333333 0.5000000
| 251 1.0000000 0.0000000
| 277 0.0000000 0.3333333# QUESTION:
# What is the interpretation of drug.m?
# Why not three columns, one for each drug?
fit.m <- lme( neg ~ drug + drug.m , dd, random = ~ 1 | id )
summary( fit.m )| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1818.357 1844.581 -902.1785
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 4.984062 3.39656
|
| Fixed effects: neg ~ drug + drug.m
| Value Std.Error DF t-value p-value
| (Intercept) 15.838657 2.121548 263 7.465614 0.0000
| drugClozapine -1.782102 0.834395 263 -2.135802 0.0336
| drugTypical 0.896756 0.627306 263 1.429535 0.1540
| drug.mClozapine 3.813141 2.872332 50 1.327542 0.1904
| drug.mTypical -2.616211 2.891131 50 -0.904909 0.3699
| Correlation:
| (Intr) drgClz drgTyp drg.mC
| drugClozapine 0.000
| drugTypical 0.000 0.552
| drug.mClozapine -0.771 -0.290 -0.160
| drug.mTypical -0.884 -0.120 -0.217 0.668
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.78234220 -0.60736018 -0.06230846 0.49922986 3.32999988
|
| Number of Observations: 318
| Number of Groups: 53wald( fit.m , 'drug.m') # Should we keep drug.m?| numDF denDF F.value p.value
| drug.m 2 50 3.783136 0.02952
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| drug.mClozapine 3.813141 2.872332 50 1.327542 0.19036 -1.956108
| drug.mTypical -2.616211 2.891131 50 -0.904909 0.36985 -8.423218
|
| Coefficients Upper 0.95
| drug.mClozapine 9.582391
| drug.mTypical 3.190796wald( fit.m , 'drug')| numDF denDF F.value p.value
| drug 4 50 4.294531 0.0046
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| drugClozapine -1.782102 0.834395 263 -2.135802 0.03362 -3.425047
| drugTypical 0.896756 0.627306 263 1.429535 0.15404 -0.338425
| drug.mClozapine 3.813141 2.872332 50 1.327542 0.19036 -1.956108
| drug.mTypical -2.616211 2.891131 50 -0.904909 0.36985 -8.423218
|
| Coefficients Upper 0.95
| drugClozapine -0.139158
| drugTypical 2.131937
| drug.mClozapine 9.582391
| drug.mTypical 3.190796L| $predicted
| Clozapine Typical
| Atypical 1 0 0
| Clozapine 1 1 0
| Typical 1 0 1
|
| $differences
| [,1] [,2] [,3]
| Clozapine - Atypical 0 1 0
| Typical - Atypical 0 0 1
| Typical - Clozapine 0 -1 1Lm <- lapply( L, function(x) cbind( x, 0, 0))
Lm| $predicted
| Clozapine Typical
| Atypical 1 0 0 0 0
| Clozapine 1 1 0 0 0
| Typical 1 0 1 0 0
|
| $differences
| [,1] [,2] [,3] [,4] [,5]
| Clozapine - Atypical 0 1 0 0 0
| Typical - Atypical 0 0 1 0 0
| Typical - Clozapine 0 -1 1 0 0wald( fit.m , Lm )| numDF denDF F.value p.value
| predicted 3 263 23.36135 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| Atypical 15.83866 2.121548 263 7.465614 <.00001 11.661277 20.01604
| Clozapine 14.05655 2.279732 263 6.165879 <.00001 9.567705 18.54540
| Typical 16.73541 2.212347 263 7.564553 <.00001 12.379247 21.09158
|
| numDF denDF F.value p.value
| differences 2 263 7.174332 0.00093
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.782102 0.834395 263 -2.135802 0.03362
| Typical - Atypical 0.896756 0.627306 263 1.429535 0.15404
| Typical - Clozapine 2.678858 0.715431 263 3.744397 0.00022
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -3.425047 -0.139158
| Typical - Atypical -0.338425 2.131937
| Typical - Clozapine 1.270156 4.087560# compare with
wald( fit.mix, L)| numDF denDF F.value p.value
| predicted 3 263 148.3786 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| Atypical 15.48282 0.878336 263 17.62744 <.00001 13.75335 17.21228
| Clozapine 14.18254 0.886412 263 15.99995 <.00001 12.43718 15.92791
| Typical 16.29784 0.788330 263 20.67387 <.00001 14.74559 17.85008
|
| numDF denDF F.value p.value
| differences 2 263 4.974805 0.00757
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.300272 0.798060 263 -1.629290 0.10445
| Typical - Atypical 0.815021 0.614065 263 1.327255 0.18558
| Typical - Clozapine 2.115292 0.682395 263 3.099806 0.00215
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.871672 0.271129
| Typical - Atypical -0.394089 2.024130
| Typical - Clozapine 0.771639 3.458945#
# QUESTION:
# With drug.m in the model, the difference between Typical and Clozapine
# is larger and more significant. Why?
#
#
# Plotting fit with Level 2 variables (here the contextual variables)
#
#
# If you want to show how the response depends on FE variables
# you need to create a prediction data frame with all the FE variables
# needed for prediction. (Level 1 or Level 2)
#
#
# If you want to show predictions for individual clusters ( popn average
# of BLUPS)
# AND
# you have Level 2 variables in your model
# you need to construct a data frame with
# ids and Level 1 variables
# then you need to merge the data.frame with the matching
# Level 2 variables.
#
#
pred <- expand.grid( year = 1:6, drug = levels( dd$drug), id = levels(dd$id))
head( pred)| year drug id
| 1 1 Atypical F2
| 2 2 Atypical F2
| 3 3 Atypical F2
| 4 4 Atypical F2
| 5 5 Atypical F2
| 6 6 Atypical F2dim( pred )| [1] 954 3# merge pred with specific Level 2 variable2 corresponding to id
ddu <- up( dd, ~ id)
head(ddu)| age1 yrsill1 status Sex Subj id what drug.m.Clozapine drug.m.Typical
| F2 29.9 10 Single F F2 F2 data 0.0000000 1.0000000
| F4 30.2 8 Single F F4 F4 data 0.0000000 0.8333333
| F11 31.8 11 Single F F11 F11 data 0.0000000 0.0000000
| F36 29.1 10 Married F F36 F36 data 0.6666667 0.1666667
| F10 26.2 9 Married F F10 F10 data 0.0000000 1.0000000
| F25 29.8 10 Single F F25 F25 data 0.0000000 0.5000000predc <- merge( pred, ddu, by = 'id')
dim( predc )| [1] 954 10head( predc )| id year drug age1 yrsill1 status Sex Subj what drug.m.Clozapine
| 1 F10 1 Atypical 26.2 9 Married F F10 data 0
| 2 F10 2 Atypical 26.2 9 Married F F10 data 0
| 3 F10 3 Atypical 26.2 9 Married F F10 data 0
| 4 F10 4 Atypical 26.2 9 Married F F10 data 0
| 5 F10 5 Atypical 26.2 9 Married F F10 data 0
| 6 F10 6 Atypical 26.2 9 Married F F10 data 0
| drug.m.Typical
| 1 1
| 2 1
| 3 1
| 4 1
| 5 1
| 6 1names( predc ) # note that drug.m matrix is preserved as a matrix| [1] "id" "year" "drug" "age1" "yrsill1" "status" "Sex"
| [8] "Subj" "what" "drug.m"# ready to predict popn average:
predc$neg <- predict( fit.m , predc, level = 0)
# predict BLUP
predc.blup <- predc
predc.blup$neg <- predict( fit.m , predc.blup, level = 1)
# Identify data frames so they can be combined and each plot differently
dd$what <- factor("data")
predc$what <- factor("fit0")
predc.blup$what <- factor('fit1')
comb <- Rbind(dd, predc, predc.blup)| [1] "age1" "yrsill1" "pos" "neg" "gen" "total" "year"
| [8] "drug" "status" "Sex" "Subj" "id" "what" "drug.m"# Note: this is all cut and paste except for the first line
gd( pch = 15:17, col = c('blue','green4','red'), cex = 1.2,
lty = c(1,1,1,2,2,2,1,1,1),lwd = 2)
xyplot( neg ~ year | id , comb, groups = what:drug,
panel = panel.superpose.2,
lwd = 2,
type = c(rep('p',3),rep('l',6)),
subset = id %in% sel,
between = list( y = 1 ),
#auto.key = list( text = c("Atypical","Clozapine","Typical"),
auto.key = list(
points = T, lines = T, columns = 3))
#
# Taking time into account
#
fit.my <- update( fit.m, . ~ . + year)
summary( fit.my ) # very significant drop with time| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1810.39 1840.334 -897.1952
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 4.992041 3.325507
|
| Fixed effects: neg ~ drug + drug.m + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.511006 2.1742271 262 8.053899 0.0000
| drugClozapine -1.153198 0.8362982 262 -1.378931 0.1691
| drugTypical -0.117037 0.6785156 262 -0.172489 0.8632
| drug.mClozapine 3.184237 2.8728859 50 1.108376 0.2730
| drug.mTypical -1.602418 2.9026725 50 -0.552049 0.5834
| year -0.477814 0.1359168 262 -3.515489 0.0005
| Correlation:
| (Intr) drgClz drgTyp drg.mC drg.mT
| drugClozapine 0.047
| drugTypical -0.093 0.397
| drug.mClozapine -0.766 -0.291 -0.116
| drug.mTypical -0.837 -0.093 -0.234 0.658
| year -0.219 -0.214 0.425 0.062 -0.099
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.50585771 -0.61471185 -0.09486059 0.56161475 3.26692262
|
| Number of Observations: 318
| Number of Groups: 53# Previous hypothesis matrix
Lm| $predicted
| Clozapine Typical
| Atypical 1 0 0 0 0
| Clozapine 1 1 0 0 0
| Typical 1 0 1 0 0
|
| $differences
| [,1] [,2] [,3] [,4] [,5]
| Clozapine - Atypical 0 1 0 0 0
| Typical - Atypical 0 0 1 0 0
| Typical - Clozapine 0 -1 1 0 0wald( fit.my )| numDF denDF F.value p.value
| 6 50 84.67485 <.00001
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| (Intercept) 17.511006 2.174227 262 8.053899 <.00001 13.229823
| drugClozapine -1.153198 0.836298 262 -1.378931 0.16909 -2.799919
| drugTypical -0.117037 0.678516 262 -0.172489 0.86319 -1.453074
| drug.mClozapine 3.184237 2.872886 50 1.108376 0.27300 -2.586125
| drug.mTypical -1.602418 2.902672 50 -0.552049 0.58337 -7.432608
| year -0.477814 0.135917 262 -3.515489 0.00052 -0.745442
|
| Coefficients Upper 0.95
| (Intercept) 21.792190
| drugClozapine 0.493524
| drugTypical 1.219001
| drug.mClozapine 8.954598
| drug.mTypical 4.227771
| year -0.210186# We only need to add year to Lm
Lmy <- Lm
Lmy <- lapply( Lm , function( x ) cbind(x, 0) )
# let's use the average year for predicted values
Lmy[[1]][,6] <- 3.5
Lmy| $predicted
| Clozapine Typical
| Atypical 1 0 0 0 0 3.5
| Clozapine 1 1 0 0 0 3.5
| Typical 1 0 1 0 0 3.5
|
| $differences
| [,1] [,2] [,3] [,4] [,5] [,6]
| Clozapine - Atypical 0 1 0 0 0 0
| Typical - Atypical 0 0 1 0 0 0
| Typical - Clozapine 0 -1 1 0 0 0#
# QUESTION: Should we do the same thing for the second matrix?
#
#
wald ( fit.my, Lmy )| numDF denDF F.value p.value
| predicted 3 262 19.26801 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| Atypical 15.83866 2.121548 262 7.465614 <.00001 11.66120 20.01611
| Clozapine 14.68546 2.280430 262 6.439777 <.00001 10.19516 19.17576
| Typical 15.72162 2.227408 262 7.058257 <.00001 11.33572 20.10752
|
| numDF denDF F.value p.value
| differences 2 262 1.034312 0.35692
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.153198 0.836298 262 -1.378931 0.16909
| Typical - Atypical -0.117037 0.678516 262 -0.172489 0.86319
| Typical - Clozapine 1.036161 0.842019 262 1.230567 0.21959
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.799919 0.493524
| Typical - Atypical -1.453074 1.219001
| Typical - Clozapine -0.621825 2.694147#
# QUESTION: What does this say about Clozapine?
#
# We don't need new pred data.frames
predc$neg <- predict( fit.my, predc, level = 0)
predc.blup$neg <- predict( fit.my, predc.blup, level = 1)
dd$what| [1] data data data data data data data data data data data data data data
| [15] data data data data data data data data data data data data data data
| [29] data data data data data data data data data data data data data data
| [43] data data data data data data data data data data data data data data
| [57] data data data data data data data data data data data data data data
| [71] data data data data data data data data data data data data data data
| [85] data data data data data data data data data data data data data data
| [99] data data data data data data data data data data data data data data
| [113] data data data data data data data data data data data data data data
| [127] data data data data data data data data data data data data data data
| [141] data data data data data data data data data data data data data data
| [155] data data data data data data data data data data data data data data
| [169] data data data data data data data data data data data data data data
| [183] data data data data data data data data data data data data data data
| [197] data data data data data data data data data data data data data data
| [211] data data data data data data data data data data data data data data
| [225] data data data data data data data data data data data data data data
| [239] data data data data data data data data data data data data data data
| [253] data data data data data data data data data data data data data data
| [267] data data data data data data data data data data data data data data
| [281] data data data data data data data data data data data data data data
| [295] data data data data data data data data data data data data data data
| [309] data data data data data data data data data data
| Levels: datapredc$what <- factor( 'fit0')
predc.blup$what <- factor( 'fit1')
rb <- Rbind( dd, predc, predc.blup)| [1] "age1" "yrsill1" "pos" "neg" "gen" "total" "year"
| [8] "drug" "status" "Sex" "Subj" "id" "what" "drug.m"# Note: this is identical to previous
xyplot( neg ~ year | id , Rbind( dd, predc, predc.blup ), groups = what:drug,
panel = panel.superpose.2,
lwd = 2,
type = c(rep('p',3),rep('l',6)),
subset = id %in% sel,
between = list( y =1 ),
auto.key = list( text = c("Atypical","Clozapine","Typical"),
points = T, lines = T, columns = 3))| [1] "age1" "yrsill1" "pos" "neg" "gen" "total" "year"
| [8] "drug" "status" "Sex" "Subj" "id" "what" "drug.m"
#
# Adding a possible autocorrelation
#
fit.myc <- update( fit.my, correlation = corAR1( form = ~ year | id ))
summary( fit.myc ) # very significant drop with time| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1806.119 1839.806 -894.0595
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 4.863878 3.477444
|
| Correlation Structure: AR(1)
| Formula: ~year | id
| Parameter estimate(s):
| Phi
| 0.2045159
| Fixed effects: neg ~ drug + drug.m + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.551586 2.1753872 262 8.068258 0.0000
| drugClozapine -1.198927 0.8940467 262 -1.341011 0.1811
| drugTypical -0.241801 0.6875520 262 -0.351684 0.7254
| drug.mClozapine 3.164622 2.8743291 50 1.100995 0.2762
| drug.mTypical -1.494575 2.8857604 50 -0.517914 0.6068
| year -0.485736 0.1503674 262 -3.230327 0.0014
| Correlation:
| (Intr) drgClz drgTyp drg.mC drg.mT
| drugClozapine 0.050
| drugTypical -0.102 0.375
| drug.mClozapine -0.758 -0.309 -0.115
| drug.mTypical -0.831 -0.088 -0.233 0.654
| year -0.244 -0.211 0.403 0.066 -0.094
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.2515034 -0.6100687 -0.1129732 0.5311779 3.0912181
|
| Number of Observations: 318
| Number of Groups: 53wald( fit.myc, 'drug.m') # what happened to the contextual effect? Why?| numDF denDF F.value p.value
| drug.m 2 50 1.94558 0.15357
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| drug.mClozapine 3.164622 2.874329 50 1.100995 0.27617 -2.608638
| drug.mTypical -1.494575 2.885760 50 -0.517914 0.60680 -7.290796
|
| Coefficients Upper 0.95
| drug.mClozapine 8.937882
| drug.mTypical 4.301645# Hypothesis matrix does not change
wald ( fit.myc, Lmy )| numDF denDF F.value p.value
| predicted 3 262 19.41893 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| Atypical 15.85151 2.109823 262 7.513194 <.00001 11.69714 20.00588
| Clozapine 14.65258 2.290340 262 6.397560 <.00001 10.14277 19.16240
| Typical 15.60971 2.216287 262 7.043182 <.00001 11.24571 19.97371
|
| numDF denDF F.value p.value
| differences 2 262 0.9123941 0.40283
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.198927 0.894047 262 -1.341011 0.18108
| Typical - Atypical -0.241801 0.687552 262 -0.351684 0.72536
| Typical - Clozapine 0.957126 0.900737 262 1.062602 0.28894
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -2.959358 0.561505
| Typical - Atypical -1.595632 1.112030
| Typical - Clozapine -0.816480 2.730731#
# Adding random Drug effects
#
fit.mycr <- update( fit.myc, random = ~ 1 + dvar( drug, id) | id)| Error in lme.formula(fixed = neg ~ drug + drug.m + year, data = dd, random = ~1 + : nlminb problem, convergence error code = 1
| message = iteration limit reached without convergence (10)fit.mycr <- update( fit.myc, random = ~ 1 + dvar( drug, id) | id,
control = list( msMaxIter = 200, msVerbose = TRUE ))| 0: 1347.6902: -0.439095 0.104470 0.163370 0.00405209 -0.313548 -1.05712 0.00000
| 1: 1345.0766: -0.437034 0.134923 0.188122 0.00338929 -0.317594 -1.06638 0.413664
| 2: 1344.8942: -0.374334 0.145914 0.194446 0.0116358 -0.313973 -1.07546 0.398566
| 3: 1344.7660: -0.352318 0.247302 0.237992 -0.0241443 -0.393254 -1.15584 0.446407
| 4: 1344.6270: -0.347783 0.282534 0.300490 0.0236578 -0.427254 -1.29728 0.427394
| 5: 1344.5999: -0.326152 0.287942 0.277751 0.0332323 -0.448758 -1.35109 0.413807
| 6: 1344.5505: -0.345930 0.301209 0.266764 0.0542005 -0.468121 -1.40492 0.429580
| 7: 1344.5190: -0.335589 0.328399 0.259043 0.0892275 -0.481017 -1.45338 0.428490
| 8: 1344.3794: -0.317595 0.423832 0.255534 0.218433 -0.710795 -1.86481 0.399681
| 9: 1344.2693: -0.332525 0.600546 0.256183 0.442419 -0.918780 -2.21538 0.435322
| 10: 1344.2435: -0.304302 0.607110 0.292488 0.476138 -0.993587 -2.32681 0.463354
| 11: 1344.2190: -0.286623 0.667102 0.313855 0.479519 -1.09430 -2.40992 0.437568
| 12: 1344.1978: -0.281306 0.687068 0.307426 0.542636 -1.18361 -2.50802 0.444023
| 13: 1344.1565: -0.254781 0.763180 0.336230 0.800035 -1.56308 -2.87658 0.453868
| 14: 1344.1309: -0.230403 0.860093 0.305654 1.04580 -1.94885 -3.24176 0.456444
| 15: 1344.1175: -0.200677 0.955116 0.328173 1.29080 -2.34311 -3.59899 0.458703
| 16: 1344.1083: -0.175381 1.01125 0.346014 1.55846 -2.77199 -3.90521 0.451133
| 17: 1344.1028: -0.150064 1.08022 0.348047 1.82266 -3.19249 -4.20386 0.452905
| 18: 1344.0977: -0.123453 1.12198 0.346923 2.10200 -3.64107 -4.44778 0.456100
| 19: 1344.0953: -0.0945663 1.17102 0.350206 2.37841 -4.08713 -4.69797 0.457902
| 20: 1344.0926: -0.0664422 1.20335 0.353877 2.66523 -4.55429 -4.89512 0.457830
| 21: 1344.0907: -0.0370896 1.23723 0.356406 2.95436 -5.02359 -5.08318 0.457725
| 22: 1344.0876: 0.0188193 1.27182 0.360981 3.47342 -5.87494 -5.32277 0.457498
| 23: 1344.0851: 0.0735086 1.31616 0.362033 3.99704 -6.72242 -5.56489 0.457951
| 24: 1344.0838: 0.122074 1.31827 0.358136 4.44943 -7.45683 -5.58626 0.458122
| 25: 1344.0799: 0.172725 1.33883 0.360442 4.89420 -8.18267 -5.72482 0.458141
| 26: 1344.0786: 0.222431 1.35804 0.362186 5.34107 -8.91033 -5.84653 0.458522
| 27: 1344.0775: 0.237173 1.35287 0.365458 5.44352 -9.07971 -5.82017 0.459272
| 28: 1344.0751: 0.277163 1.34640 0.368241 5.73547 -9.55361 -5.77370 0.460588
| 29: 1344.0711: 0.400471 1.33338 0.365608 6.66502 -11.0600 -5.68496 0.461927
| 30: 1344.0689: 0.452397 1.31052 0.360559 6.99693 -11.5934 -5.53261 0.462260
| 31: 1344.0680: 0.563252 1.29373 0.353045 7.81897 -12.9217 -5.41781 0.462363
| 32: 1344.0677: 0.540829 1.29078 0.352286 7.63601 -12.6272 -5.40902 0.461696
| 33: 1344.0677: 0.553374 1.28930 0.351942 7.73110 -12.7814 -5.39888 0.461842
| 34: 1344.0677: 0.558589 1.28841 0.351897 7.77049 -12.8456 -5.39307 0.461938
| 35: 1344.0677: 0.563458 1.28777 0.352000 7.80972 -12.9100 -5.38927 0.462038
| 36: 1344.0676: 0.576996 1.28677 0.352471 7.92628 -13.1026 -5.38376 0.462304
| 37: 1344.0676: 0.593571 1.28683 0.353216 8.08188 -13.3613 -5.38569 0.462553
| 38: 1344.0674: 0.619875 1.28946 0.354743 8.35504 -13.8185 -5.40543 0.462823
| 39: 1344.0672: 0.641469 1.29582 0.356721 8.62525 -14.2759 -5.44942 0.462876
| 40: 1344.0668: 0.660166 1.30852 0.358766 8.94454 -14.8245 -5.53708 0.462263
| 41: 1344.0668: 0.687453 1.31372 0.358835 9.27042 -15.3746 -5.57681 0.461721
| 42: 1344.0666: 0.690931 1.31619 0.358227 9.34672 -15.5071 -5.59529 0.461309
| 43: 1344.0666: 0.711359 1.31937 0.358658 9.59742 -15.9318 -5.61917 0.461155
| 44: 1344.0665: 0.757395 1.32118 0.359929 10.1220 -16.8187 -5.63908 0.461195
| 45: 1344.0663: 0.771001 1.31713 0.360429 10.2480 -17.0304 -5.61349 0.461610
| 46: 1344.0663: 0.816198 1.31473 0.361563 10.7405 -17.8618 -5.60476 0.461829
| 47: 1344.0662: 0.811773 1.31387 0.361082 10.6910 -17.7782 -5.59889 0.461688
| 48: 1344.0662: 0.828142 1.31402 0.361355 10.8888 -18.1133 -5.60218 0.461639
| 49: 1344.0661: 0.847053 1.31517 0.361849 11.1290 -18.5214 -5.61281 0.461592
| 50: 1344.0661: 0.864558 1.31738 0.362432 11.3688 -18.9296 -5.62999 0.461529
| 51: 1344.0661: 0.879582 1.32127 0.363119 11.6077 -19.3378 -5.65808 0.461407
| 52: 1344.0660: 0.897048 1.32241 0.363592 11.8477 -19.7461 -5.66923 0.461350
| 53: 1344.0660: 0.913854 1.32353 0.364124 12.0882 -20.1542 -5.67902 0.461369
| 54: 1344.0660: 0.931867 1.32412 0.364585 12.3289 -20.5621 -5.68499 0.461397
| 55: 1344.0659: 0.950631 1.32343 0.364785 12.5705 -20.9696 -5.68231 0.461421
| 56: 1344.0659: 0.969734 1.32068 0.364613 12.7998 -21.3528 -5.66493 0.461489
| 57: 1344.0659: 0.986369 1.32054 0.364771 13.0280 -21.7372 -5.66560 0.461473
| 58: 1344.0659: 1.02074 1.32130 0.365114 13.5286 -22.5789 -5.67343 0.461416
| 59: 1344.0658: 1.03767 1.32322 0.365306 13.8168 -23.0613 -5.68662 0.461357
| 60: 1344.0658: 1.05586 1.32478 0.365605 14.1044 -23.5442 -5.69778 0.461350
| 61: 1344.0658: 1.07374 1.32643 0.365897 14.3923 -24.0268 -5.70923 0.461360
| 62: 1344.0657: 1.07805 1.32824 0.366160 14.4863 -24.1838 -5.72040 0.461425
| 63: 1344.0657: 1.09577 1.33115 0.366455 14.8085 -24.7204 -5.73961 0.461444
| 64: 1344.0656: 1.12414 1.33426 0.366503 15.3161 -25.5599 -5.75899 0.461435
| 65: 1344.0656: 1.13290 1.33485 0.366524 15.4657 -25.8073 -5.76240 0.461443
| 66: 1344.0656: 1.15016 1.33573 0.366480 15.7658 -26.3015 -5.76684 0.461479
| 67: 1344.0655: 1.18517 1.33704 0.366238 16.3660 -27.2900 -5.77342 0.461506
| 68: 1344.0654: 1.21118 1.33576 0.365566 16.8064 -28.0065 -5.76184 0.461537
| 69: 1344.0654: 1.23582 1.33590 0.365122 17.2462 -28.7235 -5.76009 0.461554
| 70: 1344.0652: 1.28989 1.33551 0.363645 18.2666 -30.3719 -5.74922 0.461589
| 71: 1344.0652: 1.30736 1.33596 0.363328 18.5962 -30.9069 -5.75020 0.461595
| 72: 1344.0651: 1.32242 1.33836 0.363264 18.9261 -31.4415 -5.76309 0.461646
| 73: 1344.0651: 1.33386 1.34373 0.363616 19.2558 -31.9755 -5.79544 0.461716
| 74: 1344.0650: 1.33959 1.35226 0.364297 19.5552 -32.4591 -5.84869 0.461780
| 75: 1344.0650: 1.35258 1.35429 0.364004 19.8565 -32.9444 -5.85988 0.461788
| 76: 1344.0650: 1.36607 1.35537 0.363530 20.1589 -33.4291 -5.86477 0.461774
| 77: 1344.0649: 1.39888 1.35524 0.362066 20.8539 -34.5368 -5.85769 0.461741
| 78: 1344.0648: 1.41882 1.35237 0.360868 21.2479 -35.1566 -5.83455 0.461741
| 79: 1344.0648: 1.43641 1.35303 0.360420 21.6392 -35.7787 -5.83399 0.461800
| 80: 1344.0648: 1.45256 1.35468 0.360253 22.0310 -36.4005 -5.84003 0.461901
| 81: 1344.0647: 1.46675 1.35776 0.360459 22.4238 -37.0216 -5.85474 0.462069
| 82: 1344.0647: 1.48045 1.36071 0.360370 22.8171 -37.6424 -5.87008 0.462179
| 83: 1344.0647: 1.50115 1.36364 0.359496 23.3826 -38.5328 -5.88370 0.462204
| 84: 1344.0647: 1.51541 1.36592 0.358658 23.8140 -39.2049 -5.89390 0.462155
| 85: 1344.0646: 1.53162 1.36748 0.357970 24.2427 -39.8787 -5.89977 0.462142
| 86: 1344.0646: 1.54768 1.36903 0.357357 24.6719 -40.5523 -5.90535 0.462160
| 87: 1344.0646: 1.54497 1.36788 0.357602 24.5788 -40.4056 -5.89874 0.462206
| 88: 1344.0646: 1.54803 1.36767 0.357551 24.6498 -40.5164 -5.89642 0.462245
| 89: 1344.0646: 1.55426 1.36779 0.357369 24.8097 -40.7662 -5.89523 0.462293
| 90: 1344.0646: 1.55854 1.36820 0.357224 24.9297 -40.9538 -5.89654 0.462309
| 91: 1344.0646: 1.56344 1.36908 0.357037 25.0803 -41.1894 -5.90085 0.462308
| 92: 1344.0646: 1.56590 1.36966 0.356942 25.1599 -41.3140 -5.90394 0.462302
| 93: 1344.0646: 1.56683 1.36987 0.356910 25.1900 -41.3614 -5.90506 0.462301summary( fit.mycr )| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1809.73 1862.132 -890.865
|
| Random effects:
| Formula: ~1 + dvar(drug, id) | id
| Structure: General positive-definite, Log-Cholesky parametrization
| StdDev Corr
| (Intercept) 4.874986 (Intr) d(,i)C
| dvar(drug, id)Clozapine 3.636066 0.132
| dvar(drug, id)Typical 2.354937 0.359 0.972
| Residual 3.364987
|
| Correlation Structure: AR(1)
| Formula: ~year | id
| Parameter estimate(s):
| Phi
| 0.2271199
| Fixed effects: neg ~ drug + drug.m + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.222238 2.1614338 262 7.967969 0.0000
| drugClozapine -1.477185 1.0934029 262 -1.350998 0.1779
| drugTypical -0.443665 0.7984370 262 -0.555667 0.5789
| drug.mClozapine 3.784304 2.9985532 50 1.262043 0.2128
| drug.mTypical -0.970400 2.9348551 50 -0.330647 0.7423
| year -0.463628 0.1482115 262 -3.128150 0.0020
| Correlation:
| (Intr) drgClz drgTyp drg.mC drg.mT
| drugClozapine 0.107
| drugTypical 0.019 0.527
| drug.mClozapine -0.747 -0.428 -0.257
| drug.mTypical -0.838 -0.187 -0.318 0.668
| year -0.243 -0.145 0.342 0.059 -0.091
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.37154964 -0.55355114 -0.08730445 0.53241624 2.67344088
|
| Number of Observations: 318
| Number of Groups: 53summary( fit.myc ) # very significant drop with time| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1806.119 1839.806 -894.0595
|
| Random effects:
| Formula: ~1 | id
| (Intercept) Residual
| StdDev: 4.863878 3.477444
|
| Correlation Structure: AR(1)
| Formula: ~year | id
| Parameter estimate(s):
| Phi
| 0.2045159
| Fixed effects: neg ~ drug + drug.m + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.551586 2.1753872 262 8.068258 0.0000
| drugClozapine -1.198927 0.8940467 262 -1.341011 0.1811
| drugTypical -0.241801 0.6875520 262 -0.351684 0.7254
| drug.mClozapine 3.164622 2.8743291 50 1.100995 0.2762
| drug.mTypical -1.494575 2.8857604 50 -0.517914 0.6068
| year -0.485736 0.1503674 262 -3.230327 0.0014
| Correlation:
| (Intr) drgClz drgTyp drg.mC drg.mT
| drugClozapine 0.050
| drugTypical -0.102 0.375
| drug.mClozapine -0.758 -0.309 -0.115
| drug.mTypical -0.831 -0.088 -0.233 0.654
| year -0.244 -0.211 0.403 0.066 -0.094
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.2515034 -0.6100687 -0.1129732 0.5311779 3.0912181
|
| Number of Observations: 318
| Number of Groups: 53intervals( fit.myc ) # Beware of interpreting intervals for SDs| Approximate 95% confidence intervals
|
| Fixed effects:
| lower est. upper
| (Intercept) 13.2681190 17.5515864 21.8350537
| drugClozapine -2.9593580 -1.1989267 0.5615045
| drugTypical -1.5956322 -0.2418011 1.1120299
| drug.mClozapine -2.6086379 3.1646220 8.9378819
| drug.mTypical -7.2907957 -1.4945753 4.3016450
| year -0.7818181 -0.4857358 -0.1896535
| attr(,"label")
| [1] "Fixed effects:"
|
| Random Effects:
| Level: id
| lower est. upper
| sd((Intercept)) 3.897864 4.863878 6.069302
|
| Correlation structure:
| lower est. upper
| Phi 0.03389873 0.2045159 0.36355
| attr(,"label")
| [1] "Correlation structure:"
|
| Within-group standard error:
| lower est. upper
| 3.125290 3.477444 3.869279wald( fit.myc, 'drug.m') # what happened to the contextual effect? Why?| numDF denDF F.value p.value
| drug.m 2 50 1.94558 0.15357
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| drug.mClozapine 3.164622 2.874329 50 1.100995 0.27617 -2.608638
| drug.mTypical -1.494575 2.885760 50 -0.517914 0.60680 -7.290796
|
| Coefficients Upper 0.95
| drug.mClozapine 8.937882
| drug.mTypical 4.301645# Hypothesis matrix does not change
wald ( fit.mycr, Lmy )| numDF denDF F.value p.value
| predicted 3 262 19.72942 <.00001
|
| Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
| Atypical 15.59954 2.096422 262 7.441029 <.00001 11.471560 19.72752
| Clozapine 14.12236 2.435940 262 5.797498 <.00001 9.325844 18.91887
| Typical 15.15588 2.319962 262 6.532812 <.00001 10.587732 19.72402
|
| numDF denDF F.value p.value
| differences 2 262 0.9295739 0.39602
|
| Estimate Std.Error DF t-value p-value
| Clozapine - Atypical -1.477185 1.093403 262 -1.350998 0.17786
| Typical - Atypical -0.443665 0.798437 262 -0.555667 0.57891
| Typical - Clozapine 1.033520 0.955266 262 1.081919 0.28028
|
| Lower 0.95 Upper 0.95
| Clozapine - Atypical -3.630161 0.675790
| Typical - Atypical -2.015835 1.128505
| Typical - Clozapine -0.847456 2.914497getVarCov( fit.mycr)| Random effects variance covariance matrix
| (Intercept) dvar(drug, id)Clozapine
| (Intercept) 23.7650 2.3363
| dvar(drug, id)Clozapine 2.3363 13.2210
| dvar(drug, id)Typical 4.1188 8.3225
| dvar(drug, id)Typical
| (Intercept) 4.1188
| dvar(drug, id)Clozapine 8.3225
| dvar(drug, id)Typical 5.5457
| Standard Deviations: 4.875 3.6361 2.3549svd( getVarCov( fit.mycr) )| $d
| [1] 26.235036428 16.292475767 0.004679905
|
| $u
| [,1] [,2] [,3]
| [1,] -0.8756901 0.4731589 -0.09637141
| [2,] -0.3617502 -0.7750336 -0.51813096
| [3,] -0.3198494 -0.4188598 0.84985461
|
| $v
| [,1] [,2] [,3]
| [1,] -0.8756901 0.4731589 -0.09637141
| [2,] -0.3617502 -0.7750336 -0.51813096
| [3,] -0.3198494 -0.4188598 0.84985461coef( fit.mycr )| (Intercept) drugClozapine drugTypical drug.mClozapine drug.mTypical
| F2 12.258546 -1.477185 -0.4436649 3.784304 -0.9704
| F4 12.268543 -1.477185 -0.4436649 3.784304 -0.9704
| F11 12.103318 -1.477185 -0.4436649 3.784304 -0.9704
| F36 10.737796 -1.477185 -0.4436649 3.784304 -0.9704
| F10 13.916378 -1.477185 -0.4436649 3.784304 -0.9704
| F25 14.268775 -1.477185 -0.4436649 3.784304 -0.9704
| F24 14.069320 -1.477185 -0.4436649 3.784304 -0.9704
| F32 16.054955 -1.477185 -0.4436649 3.784304 -0.9704
| F41 14.616743 -1.477185 -0.4436649 3.784304 -0.9704
| F23 15.501026 -1.477185 -0.4436649 3.784304 -0.9704
| F27 18.491385 -1.477185 -0.4436649 3.784304 -0.9704
| F48 19.956942 -1.477185 -0.4436649 3.784304 -0.9704
| F45 23.301542 -1.477185 -0.4436649 3.784304 -0.9704
| F54 22.624043 -1.477185 -0.4436649 3.784304 -0.9704
| F55 27.632965 -1.477185 -0.4436649 3.784304 -0.9704
| M1 9.585229 -1.477185 -0.4436649 3.784304 -0.9704
| M5 11.066035 -1.477185 -0.4436649 3.784304 -0.9704
| M3 12.162397 -1.477185 -0.4436649 3.784304 -0.9704
| M21 11.970583 -1.477185 -0.4436649 3.784304 -0.9704
| M7 13.564028 -1.477185 -0.4436649 3.784304 -0.9704
| M6 12.676917 -1.477185 -0.4436649 3.784304 -0.9704
| M22 14.022968 -1.477185 -0.4436649 3.784304 -0.9704
| M14 12.240517 -1.477185 -0.4436649 3.784304 -0.9704
| M8 14.266775 -1.477185 -0.4436649 3.784304 -0.9704
| M12 14.067412 -1.477185 -0.4436649 3.784304 -0.9704
| M13 15.414158 -1.477185 -0.4436649 3.784304 -0.9704
| M15 15.614222 -1.477185 -0.4436649 3.784304 -0.9704
| M16 15.870422 -1.477185 -0.4436649 3.784304 -0.9704
| M17 15.208013 -1.477185 -0.4436649 3.784304 -0.9704
| M19 13.927730 -1.477185 -0.4436649 3.784304 -0.9704
| M37 15.511232 -1.477185 -0.4436649 3.784304 -0.9704
| M29 13.934164 -1.477185 -0.4436649 3.784304 -0.9704
| M26 14.222611 -1.477185 -0.4436649 3.784304 -0.9704
| M30 16.775732 -1.477185 -0.4436649 3.784304 -0.9704
| M20 18.361869 -1.477185 -0.4436649 3.784304 -0.9704
| M35 18.781482 -1.477185 -0.4436649 3.784304 -0.9704
| M28 18.628502 -1.477185 -0.4436649 3.784304 -0.9704
| M40 18.275107 -1.477185 -0.4436649 3.784304 -0.9704
| M31 18.875984 -1.477185 -0.4436649 3.784304 -0.9704
| M33 19.844552 -1.477185 -0.4436649 3.784304 -0.9704
| M34 18.620833 -1.477185 -0.4436649 3.784304 -0.9704
| M46 20.209872 -1.477185 -0.4436649 3.784304 -0.9704
| M43 21.772251 -1.477185 -0.4436649 3.784304 -0.9704
| M44 21.751034 -1.477185 -0.4436649 3.784304 -0.9704
| M47 18.819868 -1.477185 -0.4436649 3.784304 -0.9704
| M39 22.142500 -1.477185 -0.4436649 3.784304 -0.9704
| M38 22.982495 -1.477185 -0.4436649 3.784304 -0.9704
| M49 20.589381 -1.477185 -0.4436649 3.784304 -0.9704
| M42 23.790315 -1.477185 -0.4436649 3.784304 -0.9704
| M52 22.187680 -1.477185 -0.4436649 3.784304 -0.9704
| M50 23.238727 -1.477185 -0.4436649 3.784304 -0.9704
| M51 25.090523 -1.477185 -0.4436649 3.784304 -0.9704
| M53 28.912201 -1.477185 -0.4436649 3.784304 -0.9704
| year dvar(drug, id)Clozapine dvar(drug, id)Typical
| F2 -0.4636277 -0.48796613 -0.860257908
| F4 -0.4636277 -1.06065610 -1.208622225
| F11 -0.4636277 -0.50322608 -0.887160389
| F36 -0.4636277 -6.94564314 -4.971745231
| F10 -0.4636277 -0.32498950 -0.572938923
| F25 -0.4636277 1.68285053 0.692323493
| F24 -0.4636277 0.47050360 -0.071775539
| F32 -0.4636277 -0.11475218 -0.202301890
| F41 -0.4636277 5.44850075 3.021317777
| F23 -0.4636277 1.57650131 0.767187831
| F27 -0.4636277 0.87217103 0.676076885
| F48 -0.4636277 0.26884080 0.473951809
| F45 -0.4636277 -1.26353275 -0.083889194
| F54 -0.4636277 0.53103578 0.936187382
| F55 -0.4636277 4.25423594 3.775029411
| M1 -0.4636277 0.26067800 -0.708683863
| M5 -0.4636277 -1.80920904 -1.801707056
| M3 -0.4636277 -0.49741826 -0.876921517
| M21 -0.4636277 2.01387471 0.633916730
| M7 -0.4636277 -0.35962794 -0.634004617
| M6 -0.4636277 0.62707650 -0.132369524
| M22 -0.4636277 0.14756153 -0.272475455
| M14 -0.4636277 -1.96903346 -1.763056421
| M8 -0.4636277 -0.82553673 -0.838703711
| M12 -0.4636277 1.52453012 0.572887265
| M13 -0.4636277 -0.17774703 -0.313358413
| M15 -0.4636277 -0.15807933 -0.278685319
| M16 -0.4636277 -0.13289303 -0.234283235
| M17 -0.4636277 -1.49844263 -1.139977803
| M19 -0.4636277 -1.06218723 -1.019997218
| M37 -0.4636277 -0.73182373 -0.639315413
| M29 -0.4636277 -0.32324097 -0.569856352
| M26 -0.4636277 -0.29488455 -0.519865524
| M30 -0.4636277 -2.74000190 -1.722303004
| M20 -0.4636277 -0.42777035 -0.131918198
| M35 -0.4636277 -2.54597536 -1.377052354
| M28 -0.4636277 -0.97293608 -0.434404320
| M40 -0.4636277 -0.19580339 0.003260426
| M31 -0.4636277 0.30729869 0.374628588
| M33 -0.4636277 2.07024007 1.560561253
| M34 -0.4636277 0.55888782 0.499212330
| M46 -0.4636277 2.36614026 1.783792442
| M43 -0.4636277 -0.91140621 -0.040617374
| M44 -0.4636277 2.45098431 2.008953939
| M47 -0.4636277 0.15705841 0.276885491
| M39 -0.4636277 -1.30350745 -0.237951006
| M38 -0.4636277 0.56627417 0.998310752
| M49 -0.4636277 0.33101402 0.583559833
| M42 -0.4636277 0.98448521 1.345072401
| M52 -0.4636277 2.01118895 1.786846469
| M50 -0.4636277 -1.34961533 -0.141869333
| M51 -0.4636277 -0.51883631 0.578305862
| M53 -0.4636277 0.02480971 1.339799962zz <- ranef( fit.mycr )
names(zz) <- c("Int",'Cloz', "Typical")
names(zz)| [1] "Int" "Cloz" "Typical"if(interactive) {
library( p3d )
Init3d()
Plot3d(Int ~ Cloz +Typical , zz)
Id3d(pad = 2)
Ell3d()
}
#
# some of the potential RE variance is accounted for by within-cluster
# variance. This is analogous to having an F at or below 1 for blocks
# in a blocked randomized design.
#
fit.mycrd <- update( fit.mycr, subset = !(id %in% c("F41", "F55", "F36", "M53")))| 0: 1218.2848: -0.350809 0.409399 0.419278 0.0606497 0.197817 -0.929162 0.00000
| 1: 1214.6465: -0.348493 0.484975 0.456475 0.0664595 0.202482 -0.926218 0.481076
| 2: 1214.2631: -0.249104 0.531671 0.493657 0.0648044 0.189951 -0.928716 0.456148
| 3: 1214.0069: -0.289863 0.611199 0.556716 0.0786657 0.188032 -0.935074 0.501201
| 4: 1213.6081: -0.230851 0.784148 0.693169 0.108721 0.177443 -0.957063 0.442872
| 5: 1212.9669: -0.245281 1.13922 0.945018 0.227100 0.176893 -1.09821 0.506170
| 6: 1212.8109: -0.175260 1.43672 1.16301 0.400409 0.213672 -1.28302 0.602055
| 7: 1212.6017: -0.225219 1.75511 1.39737 0.555898 0.224998 -1.45494 0.547093
| 8: 1212.5510: -0.225556 1.98982 1.57317 0.705819 0.260614 -1.59721 0.534801
| 9: 1212.5222: -0.223161 2.21568 1.75927 0.855509 0.281178 -1.74432 0.538837
| 10: 1212.5011: -0.222257 2.52465 2.01948 1.06369 0.306900 -1.94779 0.543573
| 11: 1212.4913: -0.223239 2.80790 2.26214 1.25738 0.332138 -2.13482 0.544418
| 12: 1212.4861: -0.224315 3.08885 2.50696 1.45215 0.358484 -2.32121 0.544446
| 13: 1212.4833: -0.225112 3.36799 2.75349 1.64763 0.385423 -2.50725 0.544529
| 14: 1212.4814: -0.225793 3.71203 3.05996 1.89022 0.419445 -2.73723 0.544650
| 15: 1212.4805: -0.226192 4.04391 3.35714 2.12544 0.453042 -2.95951 0.544709
| 16: 1212.4801: -0.226454 4.37528 3.65469 2.36094 0.486998 -3.18171 0.544736
| 17: 1212.4798: -0.226606 4.70636 3.95242 2.59664 0.521234 -3.40386 0.544746
| 18: 1212.4797: -0.226701 5.03727 4.25024 2.83245 0.555608 -3.62598 0.544753
| 19: 1212.4797: -0.226753 5.36809 4.54811 3.06834 0.590088 -3.84809 0.544756
| 20: 1212.4796: -0.226785 5.69886 4.84600 3.30426 0.624616 -4.07019 0.544759
| 21: 1212.4796: -0.226803 6.02961 5.14389 3.54021 0.659182 -4.29229 0.544761
| 22: 1212.4796: -0.226813 6.36033 5.44180 3.77618 0.693764 -4.51438 0.544762
| 23: 1212.4796: -0.226819 6.69105 5.73971 4.01215 0.728360 -4.73648 0.544762
| 24: 1212.4796: -0.226822 7.02177 6.03761 4.24813 0.762960 -4.95857 0.544762
| 25: 1212.4796: -0.226824 7.35248 6.33552 4.48411 0.797565 -5.18067 0.544763
| 26: 1212.4796: -0.226825 7.68319 6.63343 4.72009 0.832171 -5.40276 0.544763
| 27: 1212.4796: -0.226826 8.01390 6.93134 4.95607 0.866779 -5.62485 0.544763
| 28: 1212.4796: -0.226826 8.34461 7.22925 5.19205 0.901388 -5.84695 0.544763
| 29: 1212.4796: -0.226826 8.67532 7.52716 5.42803 0.935997 -6.06904 0.544763summary(fit.mycrd)| Linear mixed-effects model fit by REML
| Data: dd
| Subset: !(id %in% c("F41", "F55", "F36", "M53"))
| AIC BIC logLik
| 1639.335 1690.617 -805.6676
|
| Random effects:
| Formula: ~1 + dvar(drug, id) | id
| Structure: General positive-definite, Log-Cholesky parametrization
| StdDev Corr
| (Intercept) 4.2162336556 (Intr) d(,i)C
| dvar(drug, id)Clozapine 0.0005738158 -0.001
| dvar(drug, id)Typical 0.0018088928 -0.001 0.003
| Residual 3.3605866597
|
| Correlation Structure: AR(1)
| Formula: ~year | id
| Parameter estimate(s):
| Phi
| 0.2658394
| Fixed effects: neg ~ drug + drug.m + year
| Value Std.Error DF t-value p-value
| (Intercept) 16.057880 2.0000393 242 8.028782 0.0000
| drugClozapine -1.729732 0.9491360 242 -1.822428 0.0696
| drugTypical -0.069781 0.6817590 242 -0.102354 0.9186
| drug.mClozapine 5.463603 2.6625298 46 2.052034 0.0459
| drug.mTypical -0.577079 2.6176684 46 -0.220455 0.8265
| year -0.437601 0.1504172 242 -2.909247 0.0040
| Correlation:
| (Intr) drgClz drgTyp drg.mC drg.mT
| drugClozapine 0.052
| drugTypical -0.104 0.387
| drug.mClozapine -0.738 -0.355 -0.135
| drug.mTypical -0.821 -0.098 -0.253 0.647
| year -0.265 -0.205 0.374 0.074 -0.094
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -3.3738536 -0.5858574 -0.0880750 0.5438053 2.7528831
|
| Number of Observations: 294
| Number of Groups: 49wald( fit.mycrd, "drug")| numDF denDF F.value p.value
| drug 4 46 2.250855 0.07808
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| drugClozapine -1.729732 0.949136 242 -1.822428 0.06962 -3.599354
| drugTypical -0.069781 0.681759 242 -0.102354 0.91856 -1.412720
| drug.mClozapine 5.463603 2.662530 46 2.052034 0.04588 0.104208
| drug.mTypical -0.577079 2.617668 46 -0.220455 0.82649 -5.846172
|
| Coefficients Upper 0.95
| drugClozapine 0.139891
| drugTypical 1.273158
| drug.mClozapine 10.822997
| drug.mTypical 4.692014wald( fit.mycrd, Ldiff( fit.mycrd, "drug[^\\.]", ref ="Atypical"))| numDF denDF F.value p.value
| 1 2 242 1.874021 0.15573
|
| Estimate Std.Error DF t-value p-value Lower 0.95
| pine - Atypical -1.729732 0.949136 242 -1.822428 0.06962 -3.599354
| al - Atypical -0.069781 0.681759 242 -0.102354 0.91856 -1.412720
| al - pine 1.659951 0.930129 242 1.784645 0.07557 -0.172232
|
| Upper 0.95
| pine - Atypical 0.139891
| al - Atypical 1.273158
| al - pine 3.492134# Finally: random year
fit.mycry <- update( fit.myc, random = ~ 1 + dvar( drug, id) + year | id)
fit.mycry <- update( fit.myc, random = ~ 1 + dvar( drug, id) + year| id,
control = list( msMaxIter = 200, msVerbose = TRUE ))| 0: 1342.1088: -0.529535 0.00935858 0.0825324 1.27208 0.0368797 -0.333924 -0.985909 2.61981 4.27663 -0.941670 0.00000
| 1: 1342.0089: -0.524417 0.0437405 0.128175 1.28043 0.0395296 -0.343160 -0.997026 2.61923 4.28079 -0.937800 0.0893559
| 2: 1341.9801: -0.503277 0.0485775 0.138368 1.29903 0.0375529 -0.345711 -1.00736 2.61997 4.28312 -0.936681 0.0710878
| 3: 1341.9793: -0.527049 0.0541661 0.148864 1.31127 0.0217401 -0.356637 -1.01796 2.62250 4.28602 -0.934734 0.0757131
| 4: 1341.9560: -0.511187 0.0547405 0.153031 1.31227 0.0292068 -0.354844 -1.02110 2.62178 4.28709 -0.934336 0.0737904
| 5: 1341.9415: -0.496887 0.0618125 0.174333 1.32275 0.0323529 -0.362520 -1.04153 2.62286 4.29388 -0.930675 0.0773681
| 6: 1341.9090: -0.504847 0.0759527 0.208867 1.33028 0.0358415 -0.388951 -1.09302 2.62833 4.31501 -0.919726 0.0838823
| 7: 1341.9055: -0.482569 0.0924420 0.213671 1.32621 0.00899266 -0.412967 -1.13851 2.64435 4.34355 -0.903922 0.0798157
| 8: 1341.8769: -0.490251 0.107285 0.213359 1.32996 0.0399562 -0.428768 -1.17719 2.66022 4.38531 -0.883410 0.0897290
| 9: 1341.8745: -0.497464 0.105226 0.217439 1.33475 0.0415688 -0.429223 -1.18034 2.66262 4.38995 -0.881218 0.0898676
| 10: 1341.8722: -0.491577 0.104456 0.220498 1.33588 0.0456205 -0.428093 -1.18423 2.66577 4.39629 -0.878179 0.0882957
| 11: 1341.8693: -0.495693 0.105681 0.222656 1.33654 0.0446148 -0.429741 -1.18825 2.67032 4.40404 -0.874328 0.0895832
| 12: 1341.8230: -0.496765 0.154759 0.258690 1.34041 0.0773498 -0.464767 -1.28307 2.78241 4.61297 -0.774099 0.0943019
| 13: 1341.8082: -0.479411 0.158903 0.278163 1.38281 0.141519 -0.582627 -1.50039 3.00503 5.05766 -0.569463 0.132950
| 14: 1341.7612: -0.473119 0.230124 0.275044 1.41117 0.143274 -0.648404 -1.54112 3.08186 5.21337 -0.505585 0.0908765
| 15: 1341.7277: -0.456599 0.242358 0.295407 1.37159 0.140966 -0.653899 -1.52159 3.09499 5.23952 -0.501358 0.110911
| 16: 1341.7170: -0.465162 0.235911 0.302909 1.36527 0.152375 -0.666087 -1.53320 3.11916 5.28997 -0.484084 0.104294
| 17: 1341.7128: -0.462066 0.224836 0.303253 1.36551 0.160734 -0.663829 -1.53523 3.14851 5.34278 -0.468767 0.103290
| 18: 1341.7051: -0.466706 0.231688 0.303983 1.37010 0.162482 -0.674659 -1.54246 3.17633 5.39516 -0.451391 0.107088
| 19: 1341.6609: -0.454341 0.291094 0.334134 1.37699 0.229970 -0.808177 -1.66300 3.48851 6.01441 -0.260108 0.110423
| 20: 1341.6356: -0.448586 0.363129 0.355820 1.38311 0.276903 -0.859557 -1.79418 3.82871 6.64659 -0.134038 0.104673
| 21: 1341.6235: -0.427820 0.423566 0.351099 1.42758 0.345425 -0.957342 -1.91266 4.15504 7.28551 -0.0394686 0.137098
| 22: 1341.6058: -0.419046 0.470867 0.336768 1.39255 0.402861 -1.07553 -2.00142 4.46101 7.90583 -0.0597390 0.123320
| 23: 1341.5893: -0.421636 0.488806 0.350097 1.39505 0.417211 -1.07745 -2.01891 4.66159 8.26768 -0.184556 0.118362
| 24: 1341.5673: -0.431931 0.553695 0.296058 1.39937 0.464165 -1.04450 -2.14274 5.17778 9.18236 -0.509412 0.124909
| 25: 1341.5598: -0.423394 0.585179 0.282630 1.40252 0.509422 -1.10760 -2.22493 5.46262 9.70782 -0.634050 0.131681
| 26: 1341.5573: -0.418975 0.615892 0.277159 1.40447 0.537862 -1.15130 -2.31883 5.62534 10.0004 -0.676252 0.136697
| 27: 1341.5551: -0.417407 0.637870 0.277960 1.40267 0.565757 -1.19428 -2.43355 5.79586 10.2809 -0.725693 0.141115
| 28: 1341.5540: -0.415143 0.658884 0.290562 1.40318 0.586297 -1.24494 -2.51225 5.90881 10.4692 -0.735306 0.140997
| 29: 1341.5537: -0.415918 0.660152 0.292820 1.40362 0.590137 -1.25233 -2.52591 5.93908 10.5010 -0.735399 0.139546
| 30: 1341.5535: -0.415200 0.664032 0.295138 1.40435 0.598462 -1.26741 -2.54350 6.00137 10.5751 -0.734569 0.138924
| 31: 1341.5533: -0.414202 0.668831 0.295930 1.40503 0.607536 -1.28158 -2.55928 6.06359 10.6498 -0.734587 0.138883
| 32: 1341.5528: -0.409593 0.688851 0.297709 1.40738 0.646310 -1.34001 -2.62032 6.31836 10.9430 -0.732310 0.139470
| 33: 1341.5523: -0.404595 0.709217 0.298163 1.40912 0.687911 -1.40056 -2.68034 6.58308 11.2267 -0.730936 0.140609
| 34: 1341.5517: -0.397108 0.735600 0.297849 1.41069 0.749151 -1.48840 -2.75890 6.96368 11.5860 -0.729128 0.142538
| 35: 1341.5512: -0.390327 0.754248 0.297247 1.41126 0.803250 -1.56578 -2.81790 7.29289 11.8294 -0.729137 0.144172
| 36: 1341.5506: -0.381999 0.768810 0.297221 1.41125 0.864460 -1.65584 -2.87138 7.65833 11.9995 -0.731168 0.145786
| 37: 1341.5500: -0.373075 0.773676 0.298451 1.41074 0.920234 -1.74357 -2.90346 7.98465 12.0177 -0.737061 0.146950
| 38: 1341.5495: -0.366772 0.764795 0.300738 1.40999 0.945558 -1.79263 -2.89591 8.12597 11.8501 -0.745772 0.147172
| 39: 1341.5493: -0.365842 0.751851 0.302633 1.40956 0.933905 -1.78577 -2.86723 8.05165 11.6522 -0.754077 0.146770
| 40: 1341.5493: -0.368244 0.740565 0.303238 1.40926 0.908924 -1.75279 -2.83540 7.90140 11.5091 -0.756761 0.146145
| 41: 1341.5493: -0.370239 0.737860 0.303239 1.40923 0.895327 -1.73230 -2.82628 7.82196 11.4733 -0.759036 0.145947
| 42: 1341.5493: -0.370687 0.737399 0.303046 1.40917 0.893826 -1.72899 -2.82397 7.81324 11.4729 -0.758157 0.145864
| 43: 1341.5493: -0.370650 0.737463 0.303040 1.40915 0.894183 -1.72938 -2.82405 7.81536 11.4733 -0.758124 0.145864summary(fit.mycry)| Linear mixed-effects model fit by REML
| Data: dd
| AIC BIC logLik
| 1812.699 1880.073 -888.3496
|
| Random effects:
| Formula: ~1 + dvar(drug, id) + year | id
| Structure: General positive-definite, Log-Cholesky parametrization
| StdDev Corr
| (Intercept) 5.7116322 (Intr) d(,i)C d(,i)T
| dvar(drug, id)Clozapine 4.7831944 0.435
| dvar(drug, id)Typical 2.2451318 0.199 0.481
| year 0.7303589 -0.485 -0.722 0.182
| Residual 2.9889859
|
| Correlation Structure: AR(1)
| Formula: ~year | id
| Parameter estimate(s):
| Phi
| 0.07280286
| Fixed effects: neg ~ drug + drug.m + year
| Value Std.Error DF t-value p-value
| (Intercept) 17.373201 2.1886369 262 7.937909 0.0000
| drugClozapine -1.561402 1.1228580 262 -1.390560 0.1655
| drugTypical -0.582994 0.7728203 262 -0.754372 0.4513
| drug.mClozapine 3.858564 3.0086968 50 1.282470 0.2056
| drug.mTypical -0.827937 2.9391945 50 -0.281688 0.7793
| year -0.504260 0.1670087 262 -3.019364 0.0028
| Correlation:
| (Intr) drgClz drgTyp drg.mC drg.mT
| drugClozapine 0.153
| drugTypical 0.006 0.436
| drug.mClozapine -0.759 -0.426 -0.227
| drug.mTypical -0.828 -0.153 -0.312 0.657
| year -0.276 -0.312 0.342 0.121 -0.091
|
| Standardized Within-Group Residuals:
| Min Q1 Med Q3 Max
| -2.8579153 -0.4981676 -0.1084193 0.4554007 2.6856555
|
| Number of Observations: 318
| Number of Groups: 53wald( fit.mycry )| numDF denDF F.value p.value
| 6 50 83.96209 <.00001
|
| Coefficients Estimate Std.Error DF t-value p-value Lower 0.95
| (Intercept) 17.373201 2.188637 262 7.937909 <.00001 13.063644
| drugClozapine -1.561402 1.122858 262 -1.390560 0.16554 -3.772376
| drugTypical -0.582994 0.772820 262 -0.754372 0.45130 -2.104723
| drug.mClozapine 3.858564 3.008697 50 1.282470 0.20559 -2.184582
| drug.mTypical -0.827937 2.939195 50 -0.281688 0.77935 -6.731483
| year -0.504260 0.167009 262 -3.019364 0.00278 -0.833110
|
| Coefficients Upper 0.95
| (Intercept) 21.682758
| drugClozapine 0.649573
| drugTypical 0.938736
| drug.mClozapine 9.901709
| drug.mTypical 5.075609
| year -0.175410#
# Please contact Georges Monette <georges@yorku.ca>
# with any questions or comments.
#