HLM Model

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

HLM Model

Belle
Hi

I am trying to convert SAS codes to R, but some of the result are quite different from SAS.

When I ran proc mixed, I have an option ddfm=bw followed by the model. How can I show this method in R (I am thinking that this maybe the reason that I can't get the similar results)

below is my SAS codes:

proc mixed data=test covtest empirical;
class pair grade team school;
        model score = trt pair grade school/ solution covb ddfm=bw ;
        random int /  sub=team solution type=un;
run;

I have tried both lmer and hglm, but non of them works.

Could anyone tell me how can I covert this SAS codes to R? Thanks
Reply | Threaded
Open this post in threaded view
|

Re: HLM Model

Doran, Harold
I think it should be

fm <- lmer(score ~ trt + pair + grade + school + (1|team), test)

The unstructured covariance matrix you use in proc mixed is not available in Rs function for mixed models

> -----Original Message-----
> From: [hidden email] [mailto:[hidden email]] On
> Behalf Of Belle
> Sent: Thursday, January 27, 2011 2:44 PM
> To: [hidden email]
> Subject: [R] HLM Model
>
>
> Hi
>
> I am trying to convert SAS codes to R, but some of the result are quite
> different from SAS.
>
> When I ran proc mixed, I have an option ddfm=bw followed by the model. How
> can I show this method in R (I am thinking that this maybe the reason that I
> can't get the similar results)
>
> below is my SAS codes:
>
> proc mixed data=test covtest empirical;
> class pair grade team school;
> model score = trt pair grade school/ solution covb ddfm=bw ;
> random int /  sub=team solution type=un;
> run;
>
> I have tried both lmer and hglm, but non of them works.
>
> Could anyone tell me how can I covert this SAS codes to R? Thanks
> --
> View this message in context: http://r.789695.n4.nabble.com/HLM-Model-
> tp3242999p3242999.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: HLM Model

Belle
Hi Harold:

Yes, this was the R code that I tried, and got different result from SAS.

Is that mean I cannot actually use R to run unstructured covariance matrix? How can I solve this problem if I need an unstructured covariance matrix method?

Thanks for the help.
Reply | Threaded
Open this post in threaded view
|

Re: HLM Model

Doran, Harold
Belle:

Before I provide any more help, you'll need to follow some posting guide rules and describe what you have done with some R code, describe what the problem is, and describe what you are hoping to accomplish.

The fact that the results you get from R and SAS are different doesn't mean there is a problem. If you can provide more details on what you are doing, what hasn't worked (specifically), then you will also find others on this list will provide help too.

> -----Original Message-----
> From: [hidden email] [mailto:[hidden email]] On
> Behalf Of Belle
> Sent: Thursday, January 27, 2011 3:54 PM
> To: [hidden email]
> Subject: Re: [R] HLM Model
>
>
> Hi Harold:
>
> Yes, this was the R code that I tried, and got different result from SAS.
>
> Is that mean I cannot actually use R to run unstructured covariance matrix?
> How can I solve this problem if I need an unstructured covariance matrix
> method?
>
> Thanks for the help.
> --
> View this message in context: http://r.789695.n4.nabble.com/HLM-Model-
> tp3242999p3243158.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: HLM Model

Belle
Hi Harold:

I know the outputs are different between SAS and R, but the results that I got have big difference.

Here is part of the result based on the SAS code I provided earlier:

                                Cov Parm     Subject    Estimate       Error     Value      Pr > Z

                                UN(1,1)      team        177.53      273.66        0.65       0.2583
                                Residual                    2161.15     67.1438      32.19      <.0001

                                                   Solution for Fixed Effects

                                                                                     Standard
 Effect          pairs    grade      school        Estimate       Error      DF    t Value    Pr > |t|

 Intercept                                             638.82       4.6127       5     138.49      <.0001
 trt                                                      -0.2955      3.4800       5      -0.08      0.9356
 pairs              1                                    0.1899       7.1651       5       0.03      0.9799
 pairs              2                                    31.1293      6.0636      5       5.13      0.0037
   .
   .
   .

In R:
library(lme4)
mixed<- lmer(Pre~trt+pairs+grade+school+(1|team), test)

result:

Random effects:
 Groups                 Name        Variance Std.Dev.
 team               (Intercept)     568.61  23.846  
 Residual                               2161.21  46.489

Fixed effects:
                          Estimate        Std. Error    t value
(Intercept)                540.402     43.029       12.559
trt                           7.291        13.084        0.557
pairs                       -3.535         6.150        -0.575

In random effect, the variance of team in SAS is 177.53, but it is 568.61 in R. Also I have negative estimate for trt in SAS but positive estimate for trt in R. I am wondering how this happened, and how can I solve this problem so that I can get similar result from both software.

Also does R provides result for fixed effect of each level? For example, the result of pair1, pair2,pair3,..., and grade1, grade2, grade3,...
Reply | Threaded
Open this post in threaded view
|

Re: HLM Model

Jeremy Miles-2
The empirical statement on the proc mixed line gives you robust
standard errors, I don't think you get them in R.

In SAS you specify that the predictors are to be dummy coded using the
class .  Are they factors in R?  I can't tell from the SAS output,
because the formatting has been lost.  However, it appears that in R
you did not dummy code them.  It also appears you haven't given use
all of the SAS output.

Jeremy





On 27 January 2011 15:52, Belle <[hidden email]> wrote:

>
> Hi Harold:
>
> I know the outputs are different between SAS and R, but the results that I
> got have big difference.
>
> Here is part of the result based on the SAS code I provided earlier:
>
>                                Cov Parm     Subject    Estimate       Error
> Value      Pr > Z
>
>                                UN(1,1)      team        177.53      273.66
> 0.65       0.2583
>                                Residual                    2161.15
> 67.1438      32.19      <.0001
>
>                                                   Solution for Fixed
> Effects
>
>
> Standard
>  Effect          pairs    grade      school        Estimate       Error
> DF    t Value    Pr > |t|
>
>  Intercept                                             638.82       4.6127
> 5     138.49      <.0001
>  trt                                                      -0.2955
> 3.4800       5      -0.08      0.9356
>  pairs              1                                    0.1899       7.1651
> 5       0.03      0.9799
>  pairs              2                                    31.1293      6.0636
> 5       5.13      0.0037
>   .
>   .
>   .
>
> In R:
> library(lme4)
> mixed<- lmer(Pre~trt+pairs+grade+school+(1|team), test)
>
> result:
>
> Random effects:
>  Groups                 Name        Variance Std.Dev.
>  team               (Intercept)     568.61  23.846
>  Residual                               2161.21  46.489
>
> Fixed effects:
>                          Estimate        Std. Error    t value
> (Intercept)                540.402     43.029       12.559
> trt                           7.291        13.084        0.557
> pairs                       -3.535         6.150        -0.575
>
> In random effect, the variance of team in SAS is 177.53, but it is 568.61 in
> R. Also I have negative estimate for trt in SAS but positive estimate for
> trt in R. I am wondering how this happened, and how can I solve this problem
> so that I can get similar result from both software.
>
> Also does R provides result for fixed effect of each level? For example, the
> result of pair1, pair2,pair3,..., and grade1, grade2, grade3,...
>
> --
> View this message in context: http://r.789695.n4.nabble.com/HLM-Model-tp3242999p3243475.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



--
Jeremy Miles
Psychology Research Methods Wiki: www.researchmethodsinpsychology.com

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: HLM Model

Silvano-7
In reply to this post by Belle
Hi Belle,

try this:

SAS:
proc mixed data=test noclprint noinfo covtest noitprint
method=reml;
class pair grade team school;
model score = trt pair grade school / solution ddfm=bw
notest;
random int /  sub=team solution type=un r;
run;

R:
require(nlme)
unstruct <- gls(score~trt+pair+grade+school, test,
                correlation=corSymm(form = ~ 1 |id),
                weights=varIdent(form = ~ 1|team),
method="REML")
summary(unstruct)

--------------------------------------
Silvano Cesar da Costa
Departamento de Estatística
Universidade Estadual de Londrina
Fone: 3371-4346
--------------------------------------
----- Original Message -----
From: "Belle" <[hidden email]>
To: <[hidden email]>
Sent: Thursday, January 27, 2011 5:43 PM
Subject: [R] HLM Model


>
> Hi
>
> I am trying to convert SAS codes to R, but some of the
> result are quite
> different from SAS.
>
> When I ran proc mixed, I have an option ddfm=bw followed
> by the model. How
> can I show this method in R (I am thinking that this maybe
> the reason that I
> can't get the similar results)
>
> below is my SAS codes:
>
> proc mixed data=test covtest empirical;
> class pair grade team school;
> model score = trt pair grade school/ solution covb ddfm=bw
> ;
> random int /  sub=team solution type=un;
> run;
>
> I have tried both lmer and hglm, but non of them works.
>
> Could anyone tell me how can I covert this SAS codes to R?
> Thanks
> --
> View this message in context:
> http://r.789695.n4.nabble.com/HLM-Model-tp3242999p3242999.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

survreg 3-way interaction

"Dénes TÓTH"

Dear All,

I was wondering why survreg (in survival package) can not handle three-way
interactions. I have an AFT (accelerated failure time) model with a
frailty term. If the model is restricted to two-way interactions, survreg
and gamlss.cens (with random() term) give very similar results, except for
the intercept (I do not know why it differs, anyway, if other coefficients
are pretty the same). The aftreg function in the eha package could handle
three-way interactions, but no frailty. Streg in Stata works fine with
three-way interactions and frailty. However, 1) in the two-way interaction
setting, Stata's results differ from those of gamlss and survreg, 2) I
prefer R over Stata.
So it seems that I can use gamlss.cens to analyze my data, but it would be
more comfortable to use a dedicated library. Is there a theoretical reason
that survreg does not provide the possibility of higher-order
interactions, or is merely a programmatic issue?

Best,
  Denes

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: survreg 3-way interaction

Therneau, Terry M., Ph.D.
 > I was wondering why survreg (in survival package) can not handle
> three-way interactions. I have an AFT .....

  You have given us no data to diagnose your problem. What do you mean
by "cannot handle" -- does the package print a message "no 3 way
interactions", gives wrong answers, your laptop catches on fire when you
run it, ....?
  Also, make sure you read help(frailty).  The following is a key line
"Use of gamma/ml or gaussian/reml with survreg does not lead to valid
results."

Terry Therneau

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: survreg 3-way interaction

"Dénes TÓTH"

The fire devoured my laptop, but the PC is still working... :)
Sorry for the lazyness, you are totally right.

Here goes a reproducible example, which resembles the main features of our
dataset:

# three experimental factors, 30 subjects
expgrid <- expand.grid(F1=0:1,F2=0:1,F3=0:1,id=1:30)

# set up the dataset (make it unbalanced)
inpdata <- expgrid[sample(nrow(expgrid),nrow(expgrid)*3,replace=T),]
inpdata <- inpdata[order(inpdata$id),]

# define coefficients
# fixed effects
coefs <- c(5,rnorm(7))
# random effects
subj.coefs <- rnorm(30)

# compute fixed effects
modmat <- model.matrix(logtimes~F1*F2*F3,
        data=data.frame(logtimes=0,inpdata))
fix.eff <- rep(modmat%*%coefs)

# compute logtimes (fixed effects + random effects + noise)
logtimes <- fix.eff+rep(subj.coefs,table(inpdata$id))+
        rnorm(nrow(inpdata),0,0.5)

# random censoring
cens <- runif(nrow(inpdata))<0.2
logtimes[cens] <- logtimes[cens]-rnorm(sum(cens))

# finalize dataset
inpdata$status <- !cens
inpdata$times <- exp(logtimes)
inpdata$id <- factor(inpdata$id)

# run survreg with fixed df in frailty
# (method="aic" leads to error, don't know why)
library(survival)
survreg3 <- survreg(Surv(times,status) ~
        F1*F2*F3 + frailty.gamma(id,df=29), data=inpdata,
        dist="lognormal")
# Error in survpenal.fit(X, Y, weights, offset, init = init, controlvals =
control,  :
#  Invalid pcols or pattr arg

# note, that without the frailty term or without the tree-way interaction
term, it runs smoothly


# the same with gamlss
library(gamlss.cens)
gamlss3 <- gamlss(Surv(times,status) ~
        F1*F2*F3 + random(id,df=29), data=inpdata,
        family=cens(LOGNO))






>  > I was wondering why survreg (in survival package) can not handle
>> three-way interactions. I have an AFT .....
>
>   You have given us no data to diagnose your problem. What do you mean
> by "cannot handle" -- does the package print a message "no 3 way
> interactions", gives wrong answers, your laptop catches on fire when you
> run it, ....?
>   Also, make sure you read help(frailty).  The following is a key line
> "Use of gamma/ml or gaussian/reml with survreg does not lead to valid
> results."
>
> Terry Therneau
>
>
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: HLM Model

Belle
In reply to this post by Silvano-7
Hi Silvano:

Could you tell me what "correlation=corSymm(form = ~ 1 |id)" represents?  In our case, team is random effect, trt, pairs, grade, school are fixed effect, and each team is within school.

I still got the different results from both SAS and R.
> unstruct <- gls(score~trt+pairs+grade+school, test,correlation=corSymm(form = ~ 1 |StudentID),
+          weights=varIdent(form=~1|team), method="REML")

I tried school instead of StudentID, but I got error.
> unstruct <- gls(score~trt+pairs+grade+school, test,correlation=corSymm(form = ~ 1 |school),
+          weights=varIdent(form=~1|team), method="REML")
Error in vector("double", length) : vector size specified is too large

Thanks for the help