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

