[R-sig-ME] Questions on the results from glmmPQL(MASS)
Ken Beath
ken at kjbeath.com.au
Sun Dec 7 04:14:19 CET 2008
On 07/12/2008, at 12:34 PM, zhijie zhang wrote:
> Dear Rusers,
> Sorry for re-post my question in this list. Some person has
> recommend
> that this list is more specific for me.
> I have used R,S-PLUS and SAS to analyze the sample data "bacteria"
> in
> MASS package. Their results are listed below.
> I have three questions, anybody can give me possible answers?
> Q1:From the results, we see that R get 'NAs'for AIC,BIC and logLik,
> while
> S-PLUS8.0 gave the exact values for them. Why? I had thought that R
> should
> give the same results as SPLUS here.
>
PQL is not maximum likelihood (it is an approximation which uses lme
internally and this is what generates the results) so the results
should be NA. The R and S-Plus versions have obviously diverged.
SAS nlmixed uses maximum likelihood so a log likelihood is available.
> Q2: The model to analyse the data is logity=b0+u+b1*trt
> +b2*I(week>2), but
> the results for Random effects in R/SPLUS confused me. SAS may be
> clearer.
> Random effects:
> Formula: ~1 | ID
> (Intercept) Residual
> StdDev: 1.410637 0.7800511
> Which is the random effect 'sigma'? I think it is "1.410637", but
> what
> does "0.7800511" mean? That is, i want ot know how to explain/use
> the above
> two data for Random effects.
>
I wonder if in PQL these have any real meaning, as they are obtained
from the linear mixed effects step. Try using lmer in the lme4 package.
> Q3:In SAS and other softwares, we can get p-values for the random
> effect
> 'sigma', but i donot see the p-values in the results of R/SPLUS. I
> have used
> attributes() to look for them, but no p values. Anybody knows how to
> get
> p-values for the random effect 'sigma',.
The standard answer is of the form "Just because SAS has it, doesn't
mean it is a good idea". There was a recent discussion on this list
about significance of random effects.
Ken
>
> Any suggestions or help are greatly appreciated.
> #R Results:MASS' version 7.2-44; R version 2.7.2
> library(MASS)
> summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,family =
> binomial,
> data = bacteria))
> Linear mixed-effects model fit by maximum likelihood
> Data: bacteria
> AIC BIC logLik
> NA NA NA
> Random effects:
> Formula: ~1 | ID
> (Intercept) Residual
> StdDev: 1.410637 0.7800511
> Variance function:
> Structure: fixed weights
> Formula: ~invwt
> Fixed effects: y ~ trt + I(week > 2)
> Value Std.Error DF t-value
> p-value
> (Intercept) 3.412014 0.5185033 169 6.580506 0.0000
> trtdrug -1.247355 0.6440635 47 -1.936696 0.0588
> trtdrug+ -0.754327 0.6453978 47 -1.168779 0.2484
> I(week > 2)TRUE -1.607257 0.3583379 169 -4.485311 0.0000
> Correlation:
> (Intr) trtdrg trtdr+
> trtdrug -0.598
> trtdrug+ -0.571 0.460
> I(week > 2)TRUE -0.537 0.047 -0.001
> #S-PLUS8.0: The results are the same as R except the followings:
> AIC BIC logLik
> 1113.622 1133.984 -550.8111
> #SAS9.1.3
> proc nlmixed data=b;
> parms b0=-1 b1=1 b2=1 sigma=0.4;
> yy=b0+u+b1*trt+b2*week;
> p=1/(1+exp(-yy));
> Model response~binary(p);
> Random u~normal(0,sigma) subject=id;
> Run;
> -2 Log Likelihood = 192.2
> AIC (smaller is better)=200.2
> AICC (smaller is better) =200.3
> BIC (smaller is better)= 207.8
>
> Parameter Estimates
> Standard
> Parameter Estimate Error DF t Value Pr > |t| Alpha
> Lower Upper Gradient
> b0 3.4966 0.6512 49 5.37 <.0001 0.05
> 2.1880 4.8052 -4.69E-6
> trt -0.6763 0.3352 49 -2.02 0.0491 0.05
> -1.3500 -0.00266 -0.00001
> I(week>2) -1.6132 0.4785 49 -3.37 0.0015 0.05
> -2.5747
> -0.6516 -9.35E-7
> sigma 1.5301 0.9632 49 1.59 0.1186 0.05
> -0.4054 3.4656 -2.42E-6
>
> --
> With Kind Regards,
>
> oooO:::::::::
> (..):::::::::
> :\.(:::Oooo::
> ::\_)::(..)::
> :::::::)./:::
> ::::::(_/::::
> :::::::::::::
> [***********************************************************************]
> ZhiJie Zhang ,PhD
> Dept.of Epidemiology, School of Public Health,Fudan University
> Office:Room 443, Building 8
> Office Tel./Fax.:+86-21-54237410
> Address:No. 138 Yi Xue Yuan Road,Shanghai,China
> Postcode:200032
> Email:epistat at gmail.com <Email%3Aepistat at gmail.com>
> Website: www.statABC.com
> [***********************************************************************]
> oooO:::::::::
> (..):::::::::
> :\.(:::Oooo::
> ::\_)::(..)::
> :::::::)./:::
> ::::::(_/::::
> :::::::::::::
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list