[R-sig-ME] Questions on the results from glmmPQL(MASS)

[Ken Beath]

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
>