Convert row names to columns and split them using 'stringr' package

[Patrick C.]

Hello All,

We’ve all been there when certain R functions return outputs that are
not very useable. I can think of a good example: when we use the
ranef() function to return the random effects of an nlme model, and
the effects are labeled as row names instead of columns. In the cases
where the predict() function will not work for whatever reason, it is
handy to get these random effects and merge them back to the dataframe
we’re using in order to predict manually. To do this, I have always
copy/pasted the output directly from R into a text file, and then
opened and saved it as .csv to get the columns I want. Then I would
have to read that .csv file back into R and merge it to my dataframe.
This seems like a waste of time, especially if you have to make
changes later, you’ll be doing this over and over again. The package
‘stringr’ by Hadley Wickham has offered a way around this. There are
several useful functions in this package. Below is an example of how
to make these row names into useful columns. This is done first by
converting the row names into a column and renaming it, then by using
the ‘str_split_fixed’ function to split the column of nested random
effects names by the ‘/’ symbol to make it ready to merge.

##install 'stringr' package if you haven't already, un-comment next
line if not
#install.packages("stringr")

##install 'nlme' package if you haven't already, un-comment next line
if not
#install.packages("nlme")

##install 'systemfit' package for the ppine dataset, un-comment next
line if not
#install.packages("systemfit")

##load 'stringr', 'nlme' and 'systemfit' libraries
library(stringr); library(nlme); library(systemfit)

##load 'ppine' dataset
data(ppine)

##change name of 'ppine' dataset for alteration
p.pine=ppine

##add 'Stand' column to 'p.pine' (based on the elevation measurement)
##This will create 4 different stands
p.pine$Stand=c(rep("S-One",41),rep("S-Two",39),
rep("S-Three",42),rep("S-Four",44))

##add 'Plot' column to 'p.pine'
##(I just split the number of trees in each stand into two plots)
p.pine$Plot=c(rep(1,21),rep(2,20),rep(1,18),rep(2,21),
rep(1,19),rep(2,23),rep(1,20),rep(2,24))

## NOTE: I'm adding these columns to use as random effects in the
model below.
## This is just to illustrate how to split the effects from the
ranef() output.

##fit nlme model to height and dbh with random effects 'Stand' and
'Plot'
ppine.ht=nlme(tht~4.5+exp(b0+b1*dbh^b2),data=p.pine,
fixed=b0+b1+b2~1,random=b0~1|Stand/Plot,
start=c(b0=4.3,b1=-2.5,b2=-0.6),
control=nlmeControl(minScale=1e-10,returnObject=T))

##make separate table of 'Stand' level random effects
stand.ranef=ranef(ppine.ht)$Stand

##make separate table of 'Plot' level random effects
plot.ranef=ranef(ppine.ht)$Plot

##Site random effects. This one is easy, just convert the rownames
## and change the column names
standid.ranef<-
data.frame(as.character(rownames(stand.ranef)),stand.ranef)
colnames(standid.ranef)[1]="Stand"
colnames(standid.ranef)[2]="b0.Stand"

##Plot random effects. This one is more tricky, convert to a dataframe
as before
## then change names, as before
plotid<-data.frame(as.character(rownames(plot.ranef)),plot.ranef)
colnames(plotid)[1]="effect"
colnames(plotid)[2]="b0.Plot"

##now we break out stringr to split the 'effect' column into two
##str_split_fixed(dataframe column, symbol to split by, no. of pieces
to return)
## then change names
stand.plot=data.frame(str_split_fixed(plotid$effect, "/", 2))
colnames(stand.plot) <- c("Stand","Plot")

##bind the stand and plot columns to the random effects column
## then subset the dataframe to have only the column we want
s.plot=cbind(plotid,stand.plot)
plotid.ranef=subset(s.plot,,select=c(Stand,Plot,b0.Plot))

##we can now merge the tables back to the p.pine dataset
p.pine=merge(p.pine,standid.ranef,by=c("Stand") , all=T)
p.pine=merge(p.pine,plotid.ranef,by=c("Stand","Plot"),all=T)

##next we can predict manually
p.pine$pHT1=4.5+exp((as.numeric(fixef(ppine.ht)[1])+p.pine$b0.Stand
+p.pine$b0.Plot)
+as.numeric(fixef(ppine.ht)[2])*p.pine
$dbh^as.numeric(fixef(ppine.ht)[3]))

##we get the same result as with the predict function
p.pine$pHT2=predict(ppine.ht,newdata=p.pine,na.action=na.omit)

head(p.pine)

##obviously the predict function does work in this case, but in my
experience
##it will fail sometimes. This is a quick and easy way to get the
ranefs,
##make the output useful, and predict without having to leave R