[R] Covariance structure for lme

Charles Determan Jr deter088 at umn.edu
Thu May 17 18:03:44 CEST 2012


Greetings again R users,

Some of you will likely recognize me but I hope you can help me once
more.  I have tried the mixed model mailing list for this question but have
yet to find a solution.  As such I hope someone will have another idea.

I have previously attempted to replicate the UN, CS, and AR(1)
covariance structures used in SAS PROC MIXED.  However, my efforts
have fallen short on replicating the Variance Components (VC)
structure.  I have read that it is also known as a diagonal structure.
 Below I have copied over all the models I have tried and their output
with no success.  Perhaps someone here will see my error or something
I have overlooked.  I have attached the data for this particular
model.  Thanks to all, I certainly cannot thank this help list enough.
 I you need any further information/clarification, please ask.
Cheers, Charles

dat=read.table("C:/subset.csv",sep=",",header=TRUE, na.strings=".")
attach(dat)

dat34=dat[Group %in% c("3", "4"),]
attach(dat34)
liver34=within(dat34, {
       Group=factor(Group)
       Event_name=factor(Event_name)
       Died=factor(Died)
       ID=factor(ID)
})
attach(liver34)

contrasts(liver34$Old_Event_name)=contr.sum(n=6)
contrasts(liver34$Pig_group)=contr.sum(n=2)
contrasts(liver34$Died)=contr.sum(n=2)

##What is should be from SAS
#CV var
#Type 3 Tests of Fixed Effects
#Effect          NumDF DenDF F Value Pr > F
#Group                   1      22      0.25    0.6244
#Died                    1      22      6.55    0.0179
#Group*Died              1      22      4.43    0.0470

fit.1=lme(var~Group*Died,
       random=~1|ID,
       data=dat34)
anova(fit.1, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 227.58700  <.0001
#Group          1    22   0.18320  0.6728
#Died           1    22   3.63388  0.0698
#Group:Died     1    22   3.04103  0.0951

fit.2=lme(var~Group*Died,
       data=dat34,
       random=~1|ID/Died)
anova(fit.2, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 77.99004  <.0001
#Group          1    22  1.46275  0.2393
#Died           1    22  5.84535  0.0243
#Group:Died     1    22  3.04103  0.0951

fit.3=lme(var~Group*Died,
       random=list(ID=pdSymm(~Event_name)),
       data=dat34)
anova(fit.3, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 273.10918  <.0001
#Group          1    22   0.69692  0.4128
#Died           1    22   1.43316  0.2440
#Group:Died     1    22   5.74399  0.0255

fit.4=lme(var~Group*Died,
       random=list(ID=pdSymm(~Group)),
       data=dat34)
anova(fit.4, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 235.13889  <.0001
#Group          1    22   0.15878  0.6941
#Died           1    22   3.83253  0.0631
#Group:Died     1    22   3.01222  0.0966

fit.5=lme(var~Group*Died,
       random=list(ID=pdSymm(~Group)),
       data=dat34,
       weights=varIdent(form=~1|Event_name))
anova(fit.5, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 277.16705  <.0001
#Group          1    22   0.23901  0.6298
#Died           1    22   3.99283  0.0582
#Group:Died     1    22   3.23135  0.0860

fit.6=lme(var~Group*Died,
       random=list(ID=pdSymm(~Group)),
       data=dat34,
       weights=varIdent(form=~1|Event_name))
anova(fit.6, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 277.16705  <.0001
#Group          1    22   0.23901  0.6298
#Died           1    22   3.99283  0.0582
#Group:Died     1    22   3.23135  0.0860

fit.7=lme(var~(Group*Died),
       random=list(ID=pdCompSymm(~Died)),
       data=dat34)
anova(fit.7, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 85.83799  <.0001
#Group          1    22  1.60624  0.2183
#Died           1    22  4.71795  0.0409
#Group:Died     1    22  2.65379  0.1175

fit.8=lme(var~(Group*Died),
       data=dat34,
       random=~1|ID,
       corr=corSymm())
anova(fit.8, type="marginal", adjustSigma=F)
#               numDF denDF   F-value p-value
#(Intercept)    1   101 119.54403  <.0001
#Group          1    22   4.58972  0.0435
#Died           1    22   8.01715  0.0097
#Group:Died     1    22   5.27470  0.0315

fit.9=lme(var~(Group*Died),
       data=dat34,
       random=list(ID=pdDiag(~Group*Died)),
       corr=corSymm(, ~1|ID))
#  Error in lme.formula(var ~ (Group * Died), data = dat34, random =
list(ID = pdDiag(~Group *  :
#  nlminb problem, convergence error code = 1
#  message = iteration limit reached without convergence (9)

fit.10=lme(var~(Group*Died),
       data=dat34,
       random=list(ID=pdDiag(~Group*Died)),
       corr=NULL)
anova(fit.10, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 93.90211  <.0001
#Group          1    22  1.75311  0.1991
#Died           1    22  6.84379  0.0158
#Group:Died     1    22  3.11458  0.0915

fit.11=lme(var~Group*Died,
       data=dat34,
       random=list(ID=pdDiag(~Group*Died)))
anova(fit.11, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 93.90211  <.0001
#Group          1    22  1.75311  0.1991
#Died           1    22  6.84379  0.0158
#Group:Died     1    22  3.11458  0.0915

fit.12=lme(var~Group*Died,
       data=dat34,
       random=list(ID=pdDiag(~Event_name)))
anova(fit.12, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 87.33040  <.0001
#Group          1    22  1.09661  0.3064
#Died           1    22  5.46329  0.0289
#Group:Died     1    22  2.94589  0.1001
summary(fit.12)

fit.13=lme(var~Group*Died,
       data=dat34,
       random=list(ID=pdDiag(~Group)))
anova(fit.13, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 77.99004  <.0001
#Group          1    22  1.46275  0.2393
#Died           1    22  5.84535  0.0243
#Group:Died     1    22  3.04103  0.0951

fit.14=lme(var~Group*Died,
       data=dat34,
       random=list(ID=pdDiag(~Died)))
anova(fit.14, type="marginal", adjustSigma=F)
#               numDF denDF  F-value p-value
#(Intercept)    1   101 85.83800  <.0001
#Group          1    22  1.60624  0.2183
#Died           1    22  4.71795  0.0409
#Group:Died     1    22  2.65379  0.1175

fit.15=lme(var~Group*Died,
       data=dat34,
       random=~1|ID,
       corr=corCompSymm())
anova(fit.15, type="marginal", adjustSigma=F)
#same as fit.13

fit.16=lme(var~Group*Died,
       data=dat34,
       random=~1|ID/Event_name)
anova(fit.16, type="marginal", adjustSigma=F)
#same as fit.13
#######################


More information about the R-help mailing list