[R] glmmPQL question!

[Simon Wood]

Will this do? best, Simon

## simulate some data...
set.seed(1)
joint <- c(rep(1,20),rep(2,20),rep(3,20))
time <- runif(60)+1
subject <- factor(rep(1:12,rep(5,12)))
mu <- time*joint
joint <- factor(joint)
y <- rgamma(mu,mu)

## fit model
b <- glmmPQL(y~joint*time,random=~1|subject,family=Gamma(link="identity"))

## extract fixed effect parameter estimates and covariance matrix
fix.b <- b$coefficients$fixed
V.b <- b$varFix

## Create a `prediction matrix'
pd <- data.frame(time = rep(seq(1,2,length=100),3),
                       joint=factor(c(rep(1,100),rep(2,100),rep(3,100))))

Xp <- model.matrix(~joint*time,pd)

## use it to get predictions and associated standard errors
mu <- Xp %*% fix.b
mu.se <- diag(Xp%*%V.b%*%t(Xp))^.5 ## inefficient for readability
## check this is done right
range(mu - predict(b,pd,level=0))

## produce plot
plot(pd$time[1:100],mu[1:100],main="joint==1",type="l")
lines(pd$time[1:100],mu[1:100]+2*mu.se[1:100],lty=2)
lines(pd$time[1:100],mu[1:100]-2*mu.se[1:100],lty=3)



> m-krutky at northwestern.edu wrote:
> > Hello Folks-
> >
> > Is there a way to create confidence bands with 'glmmPQL' ???
> >
> > I am performing a stroke study for Northwestern University in Chicago,
> > Illinois.  I am trying to decide a way to best plot the model which we
> > created with the glmmPQL function in R.   I would like to plot my actual
> > averaged data points within 95 % confidence intervals from the model.
> > Plotting the model is easy, but determining confidence bands is not.
> >
> > Here is my model:
> >
> > ratiomodel<-glmmPQL(ratio~as.factor(joint)*time, random = ~ 1 | subject,
> > family = Gamma(link = "identity"),alldata3)
> >
> > I am used to seeing confidence intervals from models that increase,
> > “flair out” in the y direction, at the beginning and ending time points
> > (x values) of the simulated data.  If I use 'lm' and pass the command
> > 'int = "c" ' 'to create this model I can easily find and plot this type
> > of confidence band for 'ratio~time'.  But I need to take into account
> > 'as.factor(joint)', and in fact I can produce confidence bands with 'glm'
> > by passing in 'se.fit = TRUE', but the problem is I need to make subject
> > a random variable, and take into account my ratio with the Gamma
> > distribution.
> >
> > Is there a way to create confidence bands with 'glmmPQL' ??? '
> > as.factor(joint)' has 3 levels, so I would like to produce this linear
> > model with three levels and confidence bands for comparison of the levels
> > of 'joint'.
> >
> > Any Help at all with my problem would be greatly appreciated!!
> > LJ
> >
> >
> > ------------------------------------------------------------------------
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch 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.
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.

-- 
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603  www.maths.bath.ac.uk/~sw283