library( BMA )
library( lattice )
library( latticeExtra )

setwd( "~/Egyetem/Kutatas/SulyokHIVFibrosis/" )

SI <- read.csv( "SI.csv", sep = ";" )

Labels <- c( "CD4%", "CD8%", "Age (years)", "Liver Stiffness (kPa)", "BMI (kg/m^2)", "CD4/8 ratio", "Triglyceride (mmol/L)", "Cholesterol (mmol/L)", "Sex", "Diabetes",
             "Hypertension", "Lipodystrophy", "Kaletra", "Known HIV positivity (years)", "Darunavir", "Atanazavir", "Raltegravir", "Etravirine",
             "Nevirapin", "Efavirenz", "Tenofovir", "Abacavir", "Zidovudine", "Lamivudine", "Lopinavir", "ART ever taken", "CAP (dB/m)" )

range( SI$Liver.Stiffness )
median( SI$Liver.Stiffness )
IQR( SI$Liver.Stiffness )
table( cut( SI$Liver.Stiffness, c( -Inf, 7.2, 14.6, Inf ), right = FALSE ) )
table( cut( SI$Liver.Stiffness, c( -Inf, 5.3, Inf ), right = FALSE ) )

categstat <- function( var, ref ) {
  varnotna <- var[ !is.na( var ) ]
  return( c( length( varnotna ), table( varnotna )[ ref ], prop.table( table( varnotna ) )[ ref ] ) )  
}

categosszefuzo <- function( data ) {
  return( paste0( "n=", data[ 1 ], ", ", round( data[ 2 ], 0 ), " (", round( data[ 3 ]*100, 1 ), "%)" ) )
}

contstat <- function( var, IQR = TRUE ) {
  varnotna <- var[ !is.na( var ) ]
  if ( IQR )
    return( c( length( varnotna ), round( c( mean( varnotna ), median( varnotna ), sd( varnotna ),
                                             IQR( varnotna ), range( varnotna ) ), 1 ) ) )
  else
    return( c( length( varnotna ), round( c( mean( varnotna ), median( varnotna ), sd( varnotna ), quantile( varnotna ),
                                             range( varnotna ) ), 1 ) ) )
}

contosszefuzo <- function( data, IQR = TRUE ) {
  if ( is.vector( data ) )
    data <- matrix( data, nr = 1 )
  if( IQR )
    return( paste0( "n=", data[ , 1 ], ", ", round( data[ , 2 ], 1 ), " (", round( data[ , 3 ], 1 ), ") ± ", round( data[ , 4 ], 1 ),
                    " (", round( data[ , 5 ], 1 ), ") [", round( data[ , 6 ], 1 ), "-", round( data[ , 7 ], 1 ), "]" ) )
  else
    return( paste0( "n=", data[ , 1 ], ", ", round( data[ , 2 ], 1 ), " (", round( data[ , 3 ], 1 ), ") ± ", round( data[ , 4 ], 1 ),
                    " (", round( data[ , 5 ], 1 ), "-", round( data[ , 6 ], 1 ), ") [", round( data[ , 7 ], 1 ), "-",
                    round( data[ , 8 ], 1 ), "]" ) )
}

categstatossze <- function( a, b, ref ) {
  data.frame( lapply( tapply( a, b, categstat, ref = ref ), function( x ) if(!is.null(x)) categosszefuzo(x) ),
              p = fisher.test( table( a, b ) )$p.value
  )
}

contstatossze <- function( a, b ) {
  data.frame( lapply( tapply( a, b, contstat ), function( x ) if(!is.null(x)) contosszefuzo(x) ),
              p = wilcox.test( a~b )$p.value
  )
}

ContVars <- c( "cd4", "cd8", "cd48", "age", "BMI", "triglicerid", "cholesterin", "betegsegfennalla", "Liver.Stiffness", "cap" )
CatVars <- c( names( SI )[ 15:26 ], "sex", "Diabetes", "Hypertension", "Lipodystrophy")

write.csv2( rbind( data.frame( Var = ContVars,  Stat = sapply( ContVars, function( x ) contosszefuzo( contstat( SI[[ x ]] ) ) ) ),
                   data.frame( Var = CatVars,  Stat = sapply( CatVars, function( x )
                     categosszefuzo( categstat( SI[[ x ]], ref = 2 ) ) ) ) ),
            "Table1.csv", row.names = FALSE )

cairo_pdf( "Fig2.pdf" )
xyplot( Liver.Stiffness ~ cap + cd4 + cd8 + age + BMI + cd48 + triglicerid + cholesterin + betegsegfennalla, outer = TRUE, data = SI,
        scales = list( relation = "free" ), ylab= "Liver Stiffness (kPa)", xlab = "", strip = strip.custom( factor.levels = Labels[ c( 27, 1:3, 5:8, 14 )], par.strip.text=list(cex=0.8) ),
        #par.settings = standard.theme( "pdf", color = FALSE ),
        panel = function( x, y, ... ) {
          panel.xyplot( x, y, ... )
          panel.abline( lm( y ~ x ), lwd = 1.5 )
          panel.smoother( x, y, se = FALSE, col = "red" )
        } )
dev.off()



UnivariateCors <- data.frame( Var = ContVars[ -9 ], do.call( rbind, lapply( ContVars[ -9 ], function( x ) {
  return( c( corPearson = cor.test( SI[[ x ]], SI$Liver.Stiffness )$estimate,
             pPearson = cor.test( SI[[ x ]],  SI$Liver.Stiffness )$p.value,
             corKendall = cor.test( SI[[ x ]],  SI$Liver.Stiffness, method = "kendall" )$estimate,
             pKendall = cor.test( SI[[ x ]],  SI$Liver.Stiffness, method = "kendall")$p.value ) )
} ) ), stringsAsFactors = FALSE )
write.csv2( UnivariateCors, "Table2.csv" )

UnivariateMixed <- data.frame( Var = CatVars, do.call( rbind, lapply( CatVars, function( x )
  contstatossze( SI$Liver.Stiffness, SI[[ x ]] ) ) ), stringsAsFactors = FALSE )
write.csv2(  UnivariateMixed, "Table3.csv", row.names = FALSE )

data.frame( Var = c( rep( UnivariateCors$Var, 2 ), UnivariateMixed$Var ),
            p.adj = p.adjust( c( UnivariateCors$pPearson, UnivariateCors$pKendall, UnivariateMixed$p ), "holm" ) )




SI <- na.omit( SI )


y <- SI$Liver.Stiffness
x <- data.frame( SI[,-c(4,13,26) ] )
#col 13 is a duplicated of col "lopinavir"; col 26 was only intended to be used in the descriptive part
colnames( x ) <- Labels[ -c( 4, 13, 26 ) ]


fit <- bicreg( x, y )
fit$namesx <- Labels[ -c( 4, 13, 26 ) ]
fit
summary( fit )
write.csv2( temp[1:(nrow(temp)-4),1:3], "Table4.csv" )
summary( fit, n.models = 10 )
plot( fit )

cairo_pdf( "Fig3.pdf", width = 15, height = 13 )
imageplot.bma( fit, ylab = "" )
dev.off()

fit <- lm( Liver.Stiffness ~ age + cd48, data = SI )
summary( fit )
confint( fit )

library("rms", lib.loc="D:/R/R-3.3.1/library")

fit<-ols( as.formula(paste0("y ~",paste(colnames(x),collapse="+"))),data=x,x=T,y=T )
fit
validate(fit) # one may have to run it more than one times given the small sample size
p<-pentrace(fit,seq(0,500,0.1))
plot(p)
fit<-update(fit,penalty=p$penalty)
fit