3
$\begingroup$

I'm doing partial least squares regression (PLSR), using the df below, to investigate how to predictors (catchment characteristics) influence the dependent (nitrogen in the river). In this data.frame, all predictors are areal percentage (percent of the variable in the study area) except hydrology and morphology predictors.

I'm following a number of studies (example 1 and example 2) that used:

  • The regression coefficient: to show the direction of the relationship (+ve or -ve); and
  • The variable influence on projection (VIP): to infer the most influential variables (influential variables have VIP>1)

Yet, none of these studies mentioned for which component (1 or 2 or 3..) they reported the VIP and regression coefficients.

I got the plot that shows the VIP and regression coefficients for each of the 4 components I used in my analysis.

enter image description here

As you can see in this plot, for many predictors the VIP is larger than 1 (i.e. influential) for all components. The same for many variables that don't have significant influence on the dependent for each component. Similarly, many variables have either +ve or -ve regression coefficients across all components. However, there are some variables that have different VIP values across different components. Also, some variables have different signs for regression coefficients (for example c_C3 has negative regression coefficient for component 1 but positive regression coefficients for all other components; the vice versa for c_C4).

In other words, if I decide to choose component1 to discuss the magnitude (VIP) and direction (regression coefficient sign) between predictors and dependent:

  • c_C3 doesnt influence the dependent (VIP <1)

  • c_C3 has negative relationship with the dependent

Conversely, according to component2"

  • c_C3 is influential (VIP > 1) and positively correlated with the dependent.

My question:

As none of these studies reported which component they used to show the VIP and regression coefficients for the predictors, I wonder which component should be chosen to show the magnitude (VIP) and direction (regression coefficient) between predictors and dependent?

CODE

library(pls)
library(tidyverse)
#devtools::install_github("baptiste/egg")
library(egg)

y <- as.matrix(df[,1])  
x <- as.matrix(df[,2:36])
df.pls <- mvr(y  ~ x , 
              ncomp = 4, 
              method = "oscorespls" , 
              scale = T)

summary(df.pls)

coef(df.pls, ncomp = 1:4) 


df_coef <- as.data.frame(coef(df.pls, ncomp = 1:4))

names(to_plot_df_coef)

to_plot_df_coef <- df_coef %>% 
  dplyr::mutate(variables = rownames(.)) %>% 
  dplyr::rename("comp1" = "y.1 comps",
                "comp2" = "y.2 comps",
                "comp3" = "y.3 comps",
                "comp4" = "y.4 comps") %>% 
  tidyr::gather("comp_id", "value", 1:4) %>% 
  dplyr::mutate(comp_id = factor(comp_id)) %>% 
  dplyr::mutate(category = "regression_coefficient")



ggplot(to_plot_df_coef , aes(x = variables, y = value, color = comp_id, group = comp_id))+
  geom_point()+
  geom_line()+
  geom_hline(yintercept=0, size = 0.3, linetype = 3) +
  theme(axis.text.x = element_text(angle=65,
                                   #vjust=1, 
                                   hjust=1,
                                   size = 8))+
  labs(y = "regression_coefficients")





VIP <- function(object) {
  if (object$method != "oscorespls")
    stop("Only implemented for orthogonal scores algorithm.  Refit with 'method = \"oscorespls\"'")
  if (nrow(object$Yloadings) > 1)
    stop("Only implemented for single-response models")

  SS <- c(object$Yloadings)^2 * colSums(object$scores^2)
  Wnorm2 <- colSums(object$loading.weights^2)
  SSW <- sweep(object$loading.weights^2, 2, SS / Wnorm2, "*")
  sqrt(nrow(SSW) * apply(SSW, 1, cumsum) / cumsum(SS))
}


df_vip <- VIP(df.pls)

df_vip <- as.data.frame(df_vip) %>% 
  tibble::rownames_to_column() %>% 
  dplyr::rename(comp_id = rowname) %>% 
  dplyr::mutate(comp_id  = factor(comp_id, levels=c("Comp 1", "Comp 2", "Comp 3", "Comp 4"), 
                            labels =  c("comp1", "comp2", "comp3", "comp4"))) %>% 
  tidyr::gather("variables", "value", 2:36) %>%
  dplyr::mutate(category = "VIP")


combine_vip_coef <- df_vip %>% 
  dplyr::full_join(., to_plot_df_coef) %>% 
  dplyr::mutate(variables1 = variables) %>% 
  tidyr::separate(variables1, into = c("group", "variable_id"), sep = "_", remove = FALSE) %>% 
  dplyr::mutate(group= as.factor(group)) %>% 
  dplyr::mutate(group1 = ifelse(group == "c", "carbon", 
                                   ifelse(group == "d", "drainage", 
                                          ifelse(group == "h" , "hydrology", 
                                                 ifelse(group == "l", "landuse", 
                                                        ifelse(group == "m", "morphology", 
                                                               ifelse(group == "Mr" , "mainrocks", 
                                                                      ifelse(group == "s", "soils", 
                                                                             ifelse(group == "Sr", "Subrock", NA)))))))))




p1_vip <- 
  ggplot(combine_vip_coef[combine_vip_coef$category=="VIP", ] , aes(x = variables, y=value, fill  =  comp_id))+
  geom_bar(stat="identity", position = "dodge") +
  geom_hline(yintercept = 1, size = 0.8, linetype = 3) +
  facet_grid(~group1 , scales = "free_x", space = "free_x")+
  theme_bw()+
  theme(axis.text.x = element_text(angle=65,
                                   #vjust=1, 
                                   hjust=1,
                                   size = 8))+
  labs(y = "VIP") +
  labs(x= "")

p2_coef <-
  ggplot(combine_vip_coef[combine_vip_coef$category=="regression_coefficient", ] , aes(x = variables, y = value, color = comp_id, group = comp_id))+
  geom_point()+
  geom_line()+
  geom_hline(yintercept=0, size = 0.8, linetype = 3) +
  facet_grid(~group1 , scales = "free_x", space = "free_x")+
  theme_bw()+
  theme(axis.text.x = element_text(angle=65,
                                   #vjust=1, 
                                   hjust=1,
                                   size = 8))+
  labs(y = "regression_coefficients") +
  labs(x= "")

DATA

In this .csv file or below

> dput(df)
structure(list(dependent = c(0.86397454211987, 0.787954497352421, 
0.659691072486949, 0.946075761583277, 0.728654822779142, 0.62950913750375, 
0.220547032431762, 0.644444993765386, 0.0932051430418795, 0.770186377283592, 
0.649755817096116, 1.2620621832137, 0.813883734861209, 0.756789828278448, 
0.59333732933648), h_f = c(1041.41975308642, 15773.6534246575, 
9657.94383561644, 63956.4219178082, 197778.257534247, 35966.0917808219, 
36205.2424657534, 1846.36849315069, 13306.0657534247, 43568.8849315069, 
10483.9588477366, 4790.59726027397, 2604.7397260274, 8224.63561643836, 
39813.4506849315), h_QuFl = c(326.84540048392, 6557.27843422791, 
3261.10693883351, 24068.3321704268, 65529.9097068129, 15416.1394401651, 
15774.1168214205, 807.867897832351, 5522.27019290237, 16081.754959384, 
4612.86443524532, 1528.19683548948, 872.149503125132, 2895.91238144059, 
16903.29346385), h_BF = c(714.20985306257, 9213.90344104575, 
6396.83689678293, 39897.6782605484, 132174.271350296, 20549.9523406568, 
20431.1256443329, 1038.50059531833, 7778.31571189528, 27484.0175382703, 
5871.0944124913, 3262.40042478449, 1732.59022290227, 5326.69048213606, 
22910.1572210815), h_BFIh = c(0.686, 0.584, 0.662, 0.624, 0.668, 
0.571, 0.564, 0.562, 0.585, 0.631, 0.56, 0.681, 0.665, 0.648, 
0.575), h_Ra = c(6.17674897119342, 8.23551369863014, 7.08599315068493, 
5.4904187128457, 5.67950786402841, 5.41065841802916, 10.8220662100457, 
6.46143835616438, 18.9247260273973, 7.5814924297044, 5.31914951989026, 
7.76993150684932, 6.46958904109589, 7.10945205479452, 5.00932876712329
), h_PET = c(1.81097393689986, 1.65354452054795, 1.73058219178082, 
1.60063065391574, 1.80907762557078, 1.5343526292532, 1.92253424657534, 
1.80904109589041, 1.70986301369863, 1.7879956741168, 1.52829903978052, 
1.71349315068493, 1.59897260273973, 1.75561643835616, 1.62924200913242
), h_ER = c(5.40727023319616, 7.30272260273973, 6.2747602739726, 
4.80876195399328, 4.79026281075596, 4.77028281042863, 9.85849315068493, 
5.66623287671233, 17.9095205479452, 6.6363590483057, 4.69740740740741, 
6.97095890410959, 5.73369863013699, 6.27739726027397, 4.28383561643836
), m_a = c(11.77, 163.19, 135.39, 1261.88, 3191.83, 714.57, 278.69, 
26.58, 57.35, 401.28, 220.41, 53.36, 42.31, 77.35, 744.65), m_MeEl = c(549.25, 
328.67, 451.43, 343.41, 316.31, 362.37, 470.26, 280.56, 521.56, 
308.44, 362.23, 385.66, 312.29, 466.72, 288.97), m_MeSlPe = c(40.55, 
19.53, 28.77, 19.65, 22.82, 19.67, 38.77, 19.97, 39.85, 19.1, 
22.87, 22.55, 17.87, 29.46, 25.77), m_ElRa = c(0.32, 0.36, 0.45, 
0.39, 0.5, 0.42, 0.26, 0.41, 0.45, 0.34, 0.43, 0.47, 0.32, 0.38, 
0.43), m_DrDe = c(1.7, 1.73, 1.57, 1.68, 1.64, 1.65, 1.53, 1.74, 
1.62, 1.7, 1.63, 1.96, 1.81, 1.55, 1.53), c_C2 = c(6.38465773020042, 
9.78985858901589, 1.50356582433208, 7.65339989412654, 5.49752073990677, 
3.90059314867797, 5.9409185818413, 17.4685189240362, 1.3046321051485, 
5.8492626793834, 0, 22.6041793153936, 32.247470774266, 14.1062362886855, 
0.843516522957992), c_C3 = c(54.6485735221962, 47.7105460942601, 
87.9747102844327, 51.8527177042445, 52.9410050371843, 59.9933071920204, 
59.2002531727095, 6.30611865292691, 75.3578079057062, 47.5217195679834, 
73.9175281998161, 33.0000057897949, 17.7895691750417, 50.7102968695948, 
71.2739582196542), c_C4 = c(38.9667687476034, 42.499595316724, 
10.5217238912352, 40.3249893635017, 41.4876095819823, 35.8078469269088, 
34.8588282454492, 76.2253624230368, 23.3375599891453, 46.6290177526332, 
26.0824718001839, 44.3958148948115, 49.9629600506923, 35.1834668417197, 
27.8825252573878), d_D2 = c(6.38465773020042, 2.68799095824262, 
3.83422657879301, 3.10357361725961, 5.72101253575559, 0.582584426668564, 
0.310429629241765, 23.5236369297241, 0, 7.20751236564974, 0, 
16.9040756921581, 11.776860021365, 1.14463243233264, 5.41431731847442
), d_D3 = c(0, 6.00876004990636, 4.60414912318469, 21.2125301554554, 
20.5553248044751, 23.1934190213699, 16.0656101595563, 8.76631382526795, 
0.991353257873467, 2.82722923028608, 33.1101411654781, 1.4597129328763, 
27.0133941457253, 7.18460363621184, 34.038045603972), d_D4 = c(29.812474631149, 
16.9997855063674, 16.9484733602449, 12.2962240035803, 19.4487078537099, 
13.3162195315644, 39.7673906398948, 0, 71.0194979738029, 21.4164343190374, 
23.4607809727846, 5.70504485786588, 5.4888828491522, 11.8151517616597, 
25.9226181415465), d_D5 = c(63.8028676386506, 74.3034634854837, 
74.6131509377774, 63.3876722237047, 54.2749548060594, 62.9077770203972, 
43.8565695713072, 67.7100492450079, 27.9891487683237, 68.5488240850268, 
43.4290778617374, 75.9311665170997, 55.7208629837575, 79.8556121697958, 
34.6250189360071), l_Da = c(22.2987231970424, 37.302281222996, 
0, 15.6888076579869, 15.3143227640784, 7.49552886039105, 9.21111717074486, 
15.6483412661417, 1.09281136009292, 34.6450027695119, 0.180628929489219, 
43.0778604483166, 35.4442243795419, 22.3103188721025, 1.69075281628809
), l_NaCo = c(65.1901939669953, 16.9440437901174, 16.8592486381856, 
10.0217003815836, 16.9178216725931, 7.27861971152794, 66.3473082501395, 
26.9484833920213, 64.9511228483726, 19.8836148068671, 4.13435918931213, 
32.9788279487656, 22.9070497052634, 35.5649600755622, 11.3521194977445
), l_ShBe = c(10.7608030609987, 42.8900057987031, 80.663402529751, 
69.9119944821054, 63.7520715780007, 81.8269148029647, 23.1058374218561, 
52.3253889206384, 32.7228187389004, 43.0244133098626, 90.9194088664656, 
22.7834238897341, 39.6658979608829, 39.4195030936271, 81.3124392127133
), Mr_Co = c(0, 0, 55.1128728093047, 6.51090929661756, 7.75041439460692, 
7.95893058086494, 4.14643495765247, 0, 0, 2.54938917987072, 20.7580709729255, 
0, 0.471612101548274, 0, 11.324878993859), Mr_GRAv = c(35.3896058163273, 
30.7432333308418, 5.10851905044963, 39.847735861398, 30.6826416725748, 
31.7526597888325, 19.4881703291849, 32.2650320691635, 12.3518678893064, 
39.099244961592, 7.05710760397929, 58.8850326871133, 50.9556529985623, 
41.3555074493459, 11.2694297654175), Mr_GREy = c(9.88107227563547, 
27.8058034065098, 7.06352142365074, 6.11279891418493, 13.8435635061977, 
3.3341124261228, 68.1839854780714, 27.5842838556726, 77.3296388036071, 
25.7336688950093, 10.8093669952665, 14.9315696591391, 19.6728001455373, 
3.66489820552796, 1.35155474262874), Mr_Mu = c(0, 27.9788730214221, 
22.4789874821358, 10.9120473953551, 28.3402024501966, 10.44521479672, 
8.18140923509129, 0, 10.1563537553917, 26.509122106797, 13.9141292330255, 
0, 15.2986691974923, 0, 66.5277988300636), Mr_Sa = c(54.7293219080373, 
0, 0, 34.630486683079, 16.1537362551322, 43.2860322394116, 0, 
40.1506840751639, 0, 0, 45.0619898493789, 26.1833976537477, 13.6012655568597, 
54.7388963413616, 4.60393247565709), s_SaLo = c(0, 3.62885344924762, 
49.7156992962877, 22.8742965465896, 16.1037497470103, 26.4108812066677, 
0, 8.76631382526795, 0, 2.28989887729153, 41.204133609971, 1.4597129328763, 
24.7165698755149, 0, 22.1203660572375), s_SiLo = c(68.3924324838086, 
46.7825172630918, 40.2918029150206, 50.9952195120652, 52.4616283965265, 
49.8656534812195, 29.2832265709446, 73.9920650538172, 17.3120134354451, 
47.8438174044157, 35.5069030606913, 69.5515837880942, 62.4089302957283, 
70.6441053507101, 61.3184871234248), s_St = c(9.34306648617202, 
16.9407115250502, 0, 3.02375098621997, 7.49142545780987, 0.739987935905543, 
47.9639272491623, 0, 45.1441255322013, 15.8949013302606, 2.37812587269736, 
4.51166909866183, 1.22573367622539, 1.41402832539501, 0.0517684668983939
), s_sSiLo = c(0, 10.2631744532951, 3.69657384646686, 3.51661128643063, 
5.38382097572785, 4.32919720682223, 17.0908901169596, 0, 26.3981021895912, 
10.2162013196148, 1.64750150735275, 0, 0.3886506275644, 11.6553395492559, 
0.672088258465062), Sr_Li = c(0, 0, 1.70442568792041, 5.71976303452223, 
5.0036595147534, 6.56182424061139, 4.14643495765247, 0, 2.17838680427836, 
2.86071234634266, 3.2464151051141, 0, 0.471612101548274, 0, 0.724170451915613
), Sr_SaLoSi = c(31.4915152861977, 13.0626789594875, 4.58852637794785, 
8.07133532113299, 12.3784892438516, 1.41827358332991, 6.77359973391309, 
28.7019160437988, 0.147389283434159, 20.7639813245479, 4.55174136965654, 
40.9469057406407, 34.0553180610046, 0, 10.1596460774595), Sr_SaCoCoTe = c(0, 
27.9788730214221, 21.8693540906945, 5.7796242808183, 21.9582513650753, 
1.38173157278083, 8.18140923509129, 0, 6.30678806617604, 25.9258663994387, 
4.4796460797524, 0, 15.2986691974923, 0, 48.1840998259004), Sr_SaMu = c(9.88107227563547, 
15.9350047162996, 0, 1.62740092981472, 10.4709347286845, 0, 68.1839854780714, 
27.5842838556726, 77.3296388036071, 20.9061094414991, 0, 14.9315696591391, 
19.6728001455373, 3.66489820552796, 0)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -15L), .Names = c("dependent", 
"h_f", "h_QuFl", "h_BF", "h_BFIh", "h_Ra", "h_PET", "h_ER", "m_a", 
"m_MeEl", "m_MeSlPe", "m_ElRa", "m_DrDe", "c_C2", "c_C3", "c_C4", 
"d_D2", "d_D3", "d_D4", "d_D5", "l_Da", "l_NaCo", "l_ShBe", "Mr_Co", 
"Mr_GRAv", "Mr_GREy", "Mr_Mu", "Mr_Sa", "s_SaLo", "s_SiLo", "s_St", 
"s_sSiLo", "Sr_Li", "Sr_SaLoSi", "Sr_SaCoCoTe", "Sr_SaMu"))
$\endgroup$
3
+50
$\begingroup$

Short answer: If you choose h number of components (Latent Variables, LVs) for your modelling, then the VIP scores and regression coefficients to look at are for h LVs.

Long Answer: A full scenario: Let's say you have a data having 100 samples and 500 variables. Using a method such as leave one out cross validation, calculate RMSEP or PRESS value for each number of component. Plotting the results may look something like this: enter image description here

At this point, let's say you choose 14 LVs. You can now use VIP scores of this PLSR model with 14 LVs. For using regression coefficients, again, you calculate the regression coefficients belonging to 14 LVs.

There are few things to be careful though. You usually need to use VIP scores only once, because even after you remove some variables and rebuild a model there will always be some variables whose VIP values are under 1 due to the nature of its calculation (see the reference below). Regression coefficients can be slightly tricky too. If your regression coefficients account for scaling of the data (which means they can be used on the data directly) then they are probably misleading. One needs to check that does not account for scaling. One work around is to center and scale your data prior to the PLSR. This way you can be sure that no matter which regression coefficients your R package reports, you obtain regression coefficients that are independent of scale of variables, hence comparable.

Very good paper on the subject: Performance of some variable selection methods when multicollinearity is present (PLS). It also provides calculations of VIP scores, regression coefficients, and other variable selection methods for PLS.

$\endgroup$
  • $\begingroup$ Thanks for your time and help. Just to confirm, this means there is no overall VIP value or regression coefficients for the PLSR with many components. This also means the only way to have one VIP or regression coefficient value for each predictor is to have 1 component only? $\endgroup$ – aelwan Nov 7 '17 at 9:01
  • $\begingroup$ The papers I have shared in the question, they have used more than one component but only presented one VIP value. It appeared to me that there is a way to get a whole VIP value for all components. $\endgroup$ – aelwan Nov 7 '17 at 9:03
  • 1
    $\begingroup$ The VIP scores corresponding to first component is just for the PLS model with 1 component. The VIP scores for the 2 components is for the 1st AND the 2nd component combined and so on. VIP scores are not for each individual component, it is for H number of components retained, in a cumulative way. Is it clear? $\endgroup$ – theGD Nov 7 '17 at 11:10
  • 1
    $\begingroup$ In that document, they are probably using the model with 3 components. By 3 components they are not using the 3rd alone, it actually is 1st, 2nd and 3rd component together. $\endgroup$ – theGD Nov 7 '17 at 11:12
  • 1
    $\begingroup$ If you want that information then you can use X loadings for any spesific component alone. $\endgroup$ – theGD Nov 7 '17 at 11:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.