[R-sig-ME] MCMCglmm predictions with fixed effects Splines
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Mar 26 10:38:24 CET 2013
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.
More information about the R-sig-mixed-models
mailing list