[R-sig-ME] MCMCglmm predictions with fixed effects Splines
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Mar 26 18:38:58 CET 2013
Hi Antonio,
The penalised bit of a penalised spline is achieved by having the
spline coefficients as random effects rather than fixed.
Cheers,
Jarrod
Quoting "Antonio P. Ramos" <ramos.grad.student at gmail.com> on Tue, 26
Mar 2013 10:06:57 -0700:
> Hi Jarrod,
>
> Thanks for your reply.
>
> I think I need a spline as a fixed effect only - time doesn't vary by
> CASEID and all observations are subjects for the same time trends. I would
> be nice to have a cubic spline though. Thanks
>
>
> On Tue, Mar 26, 2013 at 2:38 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk>wrote:
>
>> Hi,
>>
>> Check out the examples in ?spl: you haven't fitted a penalised spline. You
>> probably want something like:
>>
>> glm.MC.2 <- MCMCglmm(mortality.under.2 ~ maternal_age_c +
>> I(maternal_age_c^2) + birth_year + residence + maternal_educ +
>> sex + wealth, nitt=20000, thin=10, burnin=1000,random=
>> ~idv(spl(birdth_year))+CASEID,
>> prior=prior.2,data=rwanda2,**family='categorical',
>> pr=TRUE)
>>
>> Note that I have saved the random effects (pr=TRUE) because the first k
>> random effects are the spline coefficients. You will need to associate
>> these with the relevant columns of Z (rather than X) to get predictions.
>> Remember to include the fixed birth_year effect too.
>>
>> Cheers,
>>
>> Jarrod
>>
>>
>>
>> Quoting "Antonio P. Ramos" <ramos.grad.student at gmail.com> on Mon, 25 Mar
>> 2013 19:57:51 -0700:
>>
>> I think I get what is going on a little better now. Does it make sense?
>>>
>>> # getting the predictions
>>> # where the nots are?
>>> k <- 10
>>> x <- quantile(rwanda2$birth_year, 1:k/(k + 1), na.rm = T)
>>> y <- unique(rwanda2$birth_year)
>>>
>>> # creating new data for the poor
>>> pred.data <- data.frame(maternal_age_c=rep(**0,28),wealth=rep("Lowest
>>> quintile",28),
>>> sex=rep("Female",28),**residence=rep("Rural",28),
>>> "spl(birth_year)1"=ifelse(y<**1976,1,0),
>>> "spl(birth_year)2"=ifelse(y>**1975&y<1979,1,0),
>>> "spl(birth_year)3"=ifelse(y>**1978&y<1981,1,0),
>>> "spl(birth_year)4"=ifelse(y>**1980&y<1983,1,0),
>>> "spl(birth_year)5"=ifelse(y>**1982&y<1985,1,0),
>>> "spl(birth_year)6"=ifelse(y>=**1984&y<1988,1,0),
>>> "spl(birth_year)7"=ifelse(y>=**1987&y<=1990,1,0),
>>> "spl(birth_year)8"=ifelse(y>**1989&y<1993,1,0),
>>> "spl(birth_year)9"=ifelse(y>**1993&y<1996,1,0),
>>> "spl(birth_year)10"=ifelse(y>**1995,1,0),
>>> birth_order=rep(1,28),
>>> maternal_educ=rep("No education",28))
>>>
>>> pred.data$wealth <- factor(pred.data$wealth,
>>> levels=c("Lowest quintile", "Second
>>> quintile","Middle quintile","Fourth quintile","Highest quintile"))
>>> pred.data$sex <- factor(pred.data$sex, levels=c("Male","Female"))
>>> pred.data$residence <- factor(pred.data$residence,**
>>> levels=c("Rural","Urban"))
>>> pred.data$maternal_educ <- factor(pred.data$maternal_**educ,
>>> levels=c("No education", "Primary",
>>> "Secondary","Higher"))
>>>
>>>
>>> # design matrix
>>> X <- model.matrix(~ maternal_age_c + I(maternal_age_c^2) + residence +
>>> + maternal_educ +
>>> birth_order + wealth +
>>> spl.birth_year.1 +
>>> spl.birth_year.2 +
>>> spl.birth_year.3 +
>>> spl.birth_year.4 +
>>> spl.birth_year.5 +
>>> spl.birth_year.6 +
>>> spl.birth_year.7 +
>>> spl.birth_year.8 +
>>> spl.birth_year.9 +
>>> spl.birth_year.10, data=pred.data)
>>>
>>>
>>> V <- rowSums(glm.MC.2$VCV) # marginalizing over random effects
>>> beta <- glm.MC.2$Sol # fixed effects
>>> c2 <- (16*sqrt(3)/(15*pi))^2 # alterative parametrization for the logistic
>>> distribution
>>> pred<- t(plogis(t(beta%*%t(X)/sqrt(1+**c2*V))))
>>> pred <- as.data.frame(pred)
>>> colnames(pred) <- 1970:1997 # predictions for the poor for every year
>>> colSums(pred)
>>> pred.poor <- pred
>>>
>>>
>>>
>>> On Mon, Mar 25, 2013 at 5:14 PM, Antonio P. Ramos <
>>> ramos.grad.student at gmail.com> wrote:
>>>
>>> maybe the model's summary would also help:
>>>>
>>>> > summary(glm.MC.2)
>>>>
>>>> Iterations = 1001:19991
>>>> Thinning interval = 10
>>>> Sample size = 1900
>>>>
>>>> DIC: 23202.78
>>>>
>>>> G-structure: ~CASEID
>>>>
>>>> post.mean l-95% CI u-95% CI eff.samp
>>>> CASEID 1.008 0.8508 1.139 73.88
>>>>
>>>> R-structure: ~units
>>>>
>>>> post.mean l-95% CI u-95% CI eff.samp
>>>> units 1 1 1 0
>>>>
>>>> Location effects: mortality.under.2 ~ maternal_age_c +
>>>> I(maternal_age_c^2) + spl(birth_year) + residence + maternal_educ + sex +
>>>> wealth
>>>>
>>>> post.mean l-95% CI u-95% CI eff.samp
>>>> pMCMC
>>>> (Intercept) -2.2844882 -4.0378822 -0.5243228 270.8
>>>> 0.00947 **
>>>> maternal_age_c -0.0278874 -0.0396679 -0.0169772 409.5 <
>>>> 5e-04 ***
>>>> I(maternal_age_c^2) 0.0067512 0.0040534 0.0096366 369.2 <
>>>> 5e-04 ***
>>>> spl(birth_year)1 -0.0069841 -0.0172375 0.0024631 387.8
>>>> 0.15789
>>>> spl(birth_year)2 0.0257588 0.0003497 0.0511228 425.9
>>>> 0.05474 .
>>>> spl(birth_year)3 -0.0251424 -0.0871379 0.0381827 376.8
>>>> 0.41368
>>>> spl(birth_year)4 -0.0451816 -0.1068456 0.0188489 315.1
>>>> 0.17895
>>>> spl(birth_year)5 -0.0506369 -0.1256510 0.0256118 325.8
>>>> 0.20737
>>>> spl(birth_year)6 0.0571108 -0.0450529 0.1564138 229.7
>>>> 0.25368
>>>> spl(birth_year)7 -0.0993668 -0.1830175 -0.0135886 275.6
>>>> 0.01474 *
>>>> spl(birth_year)8 -0.0812237 -0.1417047 -0.0322527 273.6
>>>> 0.00316 **
>>>> spl(birth_year)9 0.0626604 0.0448635 0.0797381 300.8 <
>>>> 5e-04 ***
>>>> spl(birth_year)10 -0.0148629 -0.0234434 -0.0054245 387.5
>>>> 0.00105 **
>>>> residenceUrban -0.2141484 -0.3716601 -0.0511390 280.1
>>>> 0.00947 **
>>>> maternal_educNo education 1.0334332 0.1228220 2.0524150 207.8
>>>> 0.03263 *
>>>> maternal_educPrimary 0.7954471 -0.1554284 1.7616937 206.3
>>>> 0.10316
>>>> maternal_educSecondary -0.0083354 -0.9622503 1.0182426 195.0
>>>> 0.98000
>>>> sexMale 0.1893753 0.1072094 0.2710264 275.2 <
>>>> 5e-04 ***
>>>> wealthSecond quintile 0.1773074 0.0444755 0.3283123 398.4
>>>> 0.01368 *
>>>> wealthMiddle quintile -0.0667648 -0.2060366 0.0676636 376.9
>>>> 0.34316
>>>> wealthFourth quintile 0.0485797 -0.0913787 0.1942995 416.2
>>>> 0.47579
>>>> wealthHighest quintile -0.1339028 -0.3017483 0.0077947 338.2
>>>> 0.08000 .
>>>> ---
>>>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> On Mon, Mar 25, 2013 at 5:12 PM, Antonio P. Ramos <
>>>> ramos.grad.student at gmail.com> wrote:
>>>>
>>>> Hi all,
>>>>>
>>>>> I am trying to get some predictions from a MCMCglmm model but it is not
>>>>> working. I guess I don't really following what the model is doing with
>>>>> tje
>>>>> spl() command. Here is an example of the issue.
>>>>>
>>>>> Thanks a bunch
>>>>>
>>>>>
>>>>>
>>>>> # inve.wishart(V=1,nu=4) is equivalent to inv-gamma(shape=2,scale=2) for
>>>>> mothers random effects
>>>>> prior.2 <- list(R = list(V = 1, fix = 1), G = list(G1 = list(V = 1,nu =
>>>>> 4)))
>>>>> glm.MC.2 <- MCMCglmm(mortality.under.2 ~ maternal_age_c +
>>>>> I(maternal_age_c^2) +
>>>>> spl(birth_year) + residence + maternal_educ +
>>>>> sex + wealth,
>>>>> nitt=20000, thin=10, burnin=1000,
>>>>> random= ~CASEID, prior=prior.2,data=rwanda2,
>>>>> family='categorical')
>>>>>
>>>>>
>>>>> > # creating new data for the poor
>>>>> > pred.data <- data.frame(maternal_age_c=rep(**
>>>>> 18,25),wealth=rep("Lowest
>>>>> quintile",25),
>>>>> + sex=rep("Female",25),**
>>>>> residence=rep("Rural",25),
>>>>> + "spl(birth_year)1"=1971:1995,
>>>>> + "spl(birth_year)2"=1971:1995,
>>>>> + "spl(birth_year)3"=1971:1995,
>>>>> + "spl(birth_year)4"=1971:1995,
>>>>> + "spl(birth_year)5"=1971:1995,
>>>>> + "spl(birth_year)6"=1971:1995,
>>>>> + "spl(birth_year)7"=1971:1995,
>>>>> + "spl(birth_year)8"=1971:1995,
>>>>> + "spl(birth_year)9"=1971:1995,
>>>>> + "spl(birth_year)10"=1971:1995,
>>>>> + birth_order=rep(1,25),
>>>>> + maternal_educ=rep("No education",25))
>>>>> >
>>>>> > pred.data$wealth <- factor(pred.data$wealth,
>>>>> + levels=c("Lowest quintile", "Second
>>>>> quintile","Middle quintile","Fourth quintile","Highest quintile"))
>>>>> > pred.data$sex <- factor(pred.data$sex, levels=c("Male","Female"))
>>>>> > pred.data$residence <-
>>>>> factor(pred.data$residence,**levels=c("Rural","Urban"))
>>>>> > # pred.data$birth_year <- factor(pred.data$birth_year,
>>>>> levels=1970:1997)
>>>>> > pred.data$maternal_educ <- factor(pred.data$maternal_**educ,
>>>>> + levels=c("No education", "Primary",
>>>>> "Secondary","Higher"))
>>>>> >
>>>>> >
>>>>> > # design matrix
>>>>> > X <- model.matrix(~ maternal_age_c + I(maternal_age_c^2) + residence +
>>>>> + + maternal_educ +
>>>>> + birth_order + wealth +
>>>>> + spl.birth_year.1 +
>>>>> + spl.birth_year.2 +
>>>>> + spl.birth_year.3 +
>>>>> + spl.birth_year.4 +
>>>>> + spl.birth_year.5 +
>>>>> + spl.birth_year.6 +
>>>>> + spl.birth_year.7 +
>>>>> + spl.birth_year.8 +
>>>>> + spl.birth_year.9 +
>>>>> + spl.birth_year.10, data=pred.data)
>>>>> >
>>>>> >
>>>>> > V <- rowSums(glm.MC.2$VCV) # marginalizing over random effects
>>>>> > beta <- glm.MC.2$Sol # fixed effects
>>>>> > c2 <- (16*sqrt(3)/(15*pi))^2 # alterative parametrization for the
>>>>> logistic distribution
>>>>> > pred<- t(plogis(t(beta%*%t(X)/sqrt(1+**c2*V))))
>>>>> > pred <- as.data.frame(pred)
>>>>> > colnames(pred) <- 1971:1995 # predictions for the poor for every year
>>>>> > colSums(pred)
>>>>> 1971 1972 1973 1974 1975 1976 1977 1978
>>>>> 1979 1980 1981 1982 1983
>>>>> 1677.094 1677.093 1677.093 1677.093 1677.093 1677.093 1677.093 1677.093
>>>>> 1677.093 1677.093 1677.093 1677.093 1677.093
>>>>> 1984 1985 1986 1987 1988 1989 1990 1991
>>>>> 1992 1993 1994 1995
>>>>> 1677.092 1677.092 1677.092 1677.092 1677.092 1677.092 1677.092 1677.092
>>>>> 1677.092 1677.092 1677.092 1677.092
>>>>> >
>>>>>
>>>>>
>>>>
>>>>
>>> [[alternative HTML version deleted]]
>>>
>>>
>>>
>>
>>
>> --
>> The University of Edinburgh is a charitable body, registered in
>> Scotland, with registration number SC005336.
>>
>>
>>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list