library(nlme)
Dialyzer[1:3,]
plot(Dialyzer, outer= ~QB,layout=c(2,1))
m1Dial.lis <- nlsList(rate~SSasympOff(pressure,th1,th2,th3) | QB, data=Dialyzer)
m1Dial.lis
plot(intervals(m1Dial.lis))
t(model.matrix(~QB,data=Dialyzer))
# the model below uses dummy variable coding of the second type mentioned
# in the notes. That is, th1.(Intercept) is the th1-value for QB=200
# (the 1st level of QB) and th1.(Intercept)+th1.QB300 is the th1-value
# for QB=300 (the 2nd level of QB)
m2Dial.gnls <- gnls( rate ~ SSasympOff(pressure, th1,th2,th3),
data = Dialyzer, params = list(th1 ~ QB, th2 ~ QB, th3 ~QB),
start=c(44.99,62.22-44.99, .76,.25-.76, .22,0))
summary(m2Dial.gnls)
plot(m2Dial.gnls,grid=F,main="Residuals vs Fitteds")
plot(m2Dial.gnls, resid(.) ~ pressure, abline=0,grid=F,
main="Residuals vs x=Transmembrane Pressure")
# the model below uses dummy variable coding of the first type mentioned
# in the notes. That is, th1.QB200 is the th1-value for QB=200
# (the 1st level of QB) and th1.QB300 is the th1-value
# for QB=300 (the 2nd level of QB)
m2altDial.gnls <- gnls( rate ~ SSasympOff(pressure, th1,th2,th3),
data = Dialyzer, params = list(th1 ~ QB-1, th2 ~ QB-1, th3 ~QB-1),
start=c(44.99,62.22, .76,.25, .22,.22))
summary(m2altDial.gnls)
# Note that m2Dial.gnls and m2altDial.gnls are equivalent models differing
# only in their parameterization. Their equivalence can be seen from the fact
# that they have the same loglikelihood and AIC values
m3Dial.gnls <- update(m2Dial.gnls,weights=varPower())
anova(m2Dial.gnls,m3Dial.gnls)
m4Dial.gnls <- update(m2Dial.gnls,weights=varPower(form = ~ pressure))
anova(m2Dial.gnls,m4Dial.gnls)
anova(m3Dial.gnls,m4Dial.gnls)
plot(ACF(m4Dial.gnls, form=~1|Subject),alpha=.05)
m5Dial.gnls <- update(m4Dial.gnls,corr=corAR1(.7,form=~1|Subject))
anova(m4Dial.gnls,m5Dial.gnls)
plot(ACF(m5Dial.gnls, resType="n",form=~1|Subject),alpha=.05)
summary(m5Dial.gnls)
# The t-test of th3.QB300 in the above summary of model m5Dial.gnls
# is the Wald test of H_0:gamma_3=0 (that is, of the hypothesis
# that the offset parameter does not differ across QB levels
# Alternatively, we can get an equivalent F test via the anova function
# by specifying the hypothesis in the form A*theta=0 where A is given
# by Lmat below. Note that this Wald test is the F test given as (**) at
# the top of p. 174 of the class notes.
Lmat=matrix(c(0,0,0,0,0,1),nrow=1,byrow=T)
Lmat
anova(m5Dial.gnls,type="marginal",L=Lmat,adjustSigma=F)
# Instead of the Wald test of H_0:gamma_3=0, it is better to use a LR test.
# There are two versions of the LR test: an F version (given in two different,
# but equivalent forms on p.171 and toward the top of p.173) and a chi-square
# version (given as 2*log(lambda) in the middle of p.173). The F version is
# preferable, but they are equivalent for large samples.
# To do the LR tests, we need to fit the model under the null hypothesis
# that theta_3=0. When we do that we need to keep the variance covariance
# parameter fixed at the values from the model where theta_3 is not assumed
# equal to 0 (model m5Dial.gnls)
# The following model call should work but doesn't due to some unknown bug
# in the update.gnls function.
#m6Dial.gnls <- update(m5Dial.gnls, params=list(th1 ~ QB, th2 ~ QB, th3 ~ 1),
# weights=varPower(form = ~ pressure,fixed=.5849288),
# corr=corAR1(value=.7432391,fixed=T,form=~1|Subject),
# start=c(55.11,7.77,.37,-.16,.22))
# So, we'll try to fit this model without using the update function
m6Dial.gnls <- gnls( rate ~ SSasympOff(pressure, th1,th2,th3),
data = Dialyzer, params=list(th1 ~ QB, th2 ~ QB, th3 ~ 1),
weights=varPower(form = ~ pressure,fixed=.5849288),
corr=corAR1(value=.7432391,fixed=T,form=~1|Subject),
start=c(47,15,.54,-.33,.22))
# The following statement gives the chi-square version of the LR test,
# but note that it gets the difference in df of the two models wrong, and
# hence also gives the wrong p-value
anova(m5Dial.gnls,m6Dial.gnls) # gives right test stat, wrong df and p-val
# This can be fixed by computing the p-value from the chi-square
# distribution and correctly specifying the df as 1, not 3:
1-pchisq(2*(logLik(m5Dial.gnls)-logLik(m6Dial.gnls)),1) #correct p-value based on 1 df
#A way to get around the above problem that the anova function gets the wrong
# df and p-value for this test is as follows: refit model m5 fixing the
# covariance parameters
# at their estimated values. Call this model m5a (will be exactly the same as
# model m5, but the software will no longer regard the covariance parameters
# as free parameters and won't count them in the df). Then test model m6 vs
# m5a with the anova function.
m5aDial.gnls <- update(m4Dial.gnls,weights=varPower(form=~pressure, fixed=.5849288),
corr=corAR1(value=.7432391,fixed=T,form=~1|Subject))
logLik(m5Dial.gnls)
logLik(m5aDial.gnls) # should be the same because models are identical
# this will give the correct chi-square version of the LRT
anova(m6Dial.gnls,m5aDial.gnls)
# It is a bit harder to get the F version of the LR test, but here it is
# (cf. the formula on the top of p.173):
lam <- as.numeric(exp(logLik(m5Dial.gnls)-logLik(m6Dial.gnls)))
F <- (lam^(2/nrow(Dialyzer))-1)*( nrow(Dialyzer)-length(coef(m5Dial.gnls)))/1
1-pf( F,1,nrow(Dialyzer)-length(coef(m5Dial.gnls)) )