[R-sig-ME] MCMCglmm predictions with fixed effects Splines
Joshua Wiley
jwiley.psych at gmail.com
Tue Mar 26 02:30:53 CET 2013
Realized it may not be clear how to extract the design matrix:
model$X
it is stored in the X element of the model
On Mon, Mar 25, 2013 at 6:25 PM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
> Hi Antonio,
>
> Have you looked at the documentation for spl? The default is quantile
> smooths. If you read Simon Wood's book, it provides much more detail
> on how smooths are constructed from basis functions.
>
> In particular, from the docs, that is not functioning as a smooth term
> in gam (which you recently posted about on Rhelp so I assume you are
> familiar with), where the approximate df are estimated using thin
> plate splines, but it is fixed at 10, (the default k parameter).
> Anyway, you cannot simply feed the same birth year values to every one
> of the smooth parameters and use those. You need to create the proper
> model matrix based on the data to be predicted from. Extract the
> design matrix from your MCMCglmm object, and you may find that
> instructive.
>
> After you have tried that, if you are still struggling, email back or
> since you are at UCLA, you can come to the UCLA IDRE statistical
> consulting group where I or Neal can help you set it up :)
>
> Cheers,
>
> Josh
>
>
> 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]]
>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>
>
> --
> Joshua Wiley
> Ph.D. Student, Health Psychology
> University of California, Los Angeles
> http://joshuawiley.com/
> Senior Analyst - Elkhart Group Ltd.
> http://elkhartgroup.com
--
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://joshuawiley.com/
Senior Analyst - Elkhart Group Ltd.
http://elkhartgroup.com
More information about the R-sig-mixed-models
mailing list