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]]

