title: geepack: parameter estimates change sign depending on correlation structure

text: I'm working on a dataset with the following variables:               Y:  a boolean telling whether the subject has experienced a seizure     ID:  the id of the subject     Sess:  the subject's session number (a subject has been observed multiple times)     X3:  numeric measurements of the subject's behavior during the session      X4:  "     X6:  "      The idea is to see if the behavior measurements(X3,X4,X6) can be used predict the Seizure status (Y) in the population. If there were just one session per subject I'd model the data with a logistic regression, but since each subject has multiple sessions the observations cannot be assumed independent. It seems a GEE logistic model makes the most sense for this dataset, but I'm having trouble understanding the results when compared to a regular glm. When introducing correlation structure to the GEE equation the signs of the coefficients are reversed, and the predicted values are opposite of what they would be with an independent correlation structure.               library("geepack")     #regular logistic model     mod0 = glm(Y~X3+X4+X6, family=binomial("logit"), data=mice)     #gee logistic model, independent correlation structure     mod1 = geeglm(Y~X3+X4+X6, id=ID, family=binomial("logit"), corstr="indep", scale.fix=T, waves=Sess, data=mice)     #gee logistic model, exchangeable correlation structure     mod2 = geeglm(Y~X3+X4+X6, id=ID, family=binomial("logit"), corstr="exchangeable", scale.fix=T, waves=Sess, data=mice)      As expected, the parameter estimates of the glm model(mod0) and the independent correlation gee model(mod1) are the same. But when an exchangeable correlation structure(mod2) is introduced, the estimates are completely different and change sign.               >      > summary(mod1)                              Call:     geeglm(formula = Y ~ X3 + X4 + X6, family = binomial("logit"),          data = mice, id = ID, waves = Sess, corstr = "indep", scale.fix = T)           Coefficients:                 Estimate Std.err  Wald Pr(>|W|)         (Intercept)   -1.494   0.459 10.59  0.00114 **      X3            -2.377   0.626 14.42  0.00015 ***     X4             1.090   0.398  7.50  0.00619 **      X6             1.233   0.403  9.36  0.00222 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Scale is fixed.          Correlation: Structure = independenceNumber of clusters:   28   Maximum cluster size: 4      > summary(mod2)          Call:     geeglm(formula = Y ~ X3 + X4 + X6, family = binomial("logit"),          data = mice, id = ID, waves = Sess, corstr = "exchangeable",          scale.fix = T)           Coefficients:                 Estimate Std.err Wald Pr(>|W|)        (Intercept)  -0.8928  0.4255 4.40   0.0359 *      X3            0.1059  0.0383 7.64   0.0057 **     X4           -0.0427  0.0299 2.04   0.1535        X6           -0.0528  0.0213 6.16   0.0131 *      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Scale is fixed.          Correlation: Structure = exchangeable  Link = identity           Estimated Correlation Parameters:           Estimate Std.err     alpha     1.09   0.133     Number of clusters:   28   Maximum cluster size: 4       My biggest concern is that the predicted Y's of the two GEE models show opposite trends, i.e. behaviors that would predict a positive seizure status in mod1, predict a negative seizure status in mod2.               plot(mod1$fitted, mod2$fitted)      ![scatter plot of fitted values from mod1, mod2](http://i.stack.imgur.com/nAogH.png) Also, what's with the estimated correlation parameter of alpha = 1.09 in mod2? Unless I'm interpreting this wrong, shouldn't this always fall between -1 and 1?               Correlation: Structure = exchangeable  Link = identity           Estimated Correlation Parameters:           Estimate Std.err     alpha     1.09   0.133      It seems odd to me that the results are completely flipped depending on whether or not the observations are assumed to have dependence structure. Can anyone else offer insight on this behavior? Here is the data:               > dput(mice)     structure(list(Y = c(0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L,      0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L     ), ID = c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L,      5L, 5L, 6L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L,      11L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 15L,      15L, 15L, 15L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 18L, 18L, 18L,      19L, 19L, 19L, 20L, 20L, 20L, 21L, 21L, 21L, 22L, 22L, 22L, 23L,      23L, 23L, 23L, 24L, 24L, 24L, 24L, 25L, 25L, 25L, 25L, 26L, 26L,      26L, 26L, 27L, 27L, 27L, 27L, 28L, 28L, 28L, 28L), Sess = c(2L,      3L, 4L, 2L, 3L, 4L, 2L, 3L, 4L, 2L, 3L, 4L, 2L, 3L, 4L, 2L, 3L,      4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 1L, 3L, 4L, 1L, 3L, 4L, 1L,      2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 2L,      3L, 4L, 2L, 3L, 4L, 2L, 3L, 4L, 2L, 3L, 4L, 1L, 2L, 3L, 1L, 2L,      3L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L,      4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), X3 = c(0.511015060794264,      0.898356533693696, 0.798280430157052, 1.31144372617517, 0.829923189452201,      0.289089506643144, -0.763028944257538, -0.944459588217789, -1.16474609928919,      -0.182524267014845, -0.338967193889711, -0.896037509887988, 0.00426073081308205,      0.0576592603165749, -1.4984737260339, 1.34752684212464, 0.106461438449047,      -0.108424579472268, -2.85991432039569, -0.230115838261355, -1.54479536845993,      -1.23693649938367, -1.53704616612456, -1.04825100254239, 0.142768659484482,      0.28358135516745, -0.236302896321009, 0.708743856986942, -0.507503006972081,      0.401550711842527, -0.16928449007327, 0.867816722958898, 0.487459858572122,      1.35172112260613, 0.14742652989871, 0.742155288774287, 0.348552056119878,      -0.82489952485408, 0.0366834636917457, -0.731010479377091, 0.979544093857171,      1.4161996712129, 0.661035838980077, 0.600235250313596, -1.10872641912335,      -0.212101744145196, -0.919575135240643, -0.813993077336991, -0.547068540188791,      -0.0260198210967738, -0.0962240349391501, -0.251025721625606,      0.894913664382802, -0.21993004239326, 0.0628839847717805, 1.77763503559622,      0.718459471596243, 0.984412886705251, 0.77603470471174, 0.486187732642953,      1.78012655684609, 1.31622243756713, 1.29635178661133, 0.427995111986702,      0.993748401511881, 0.387623239882247, 0.42006794384777, -0.815889182132972,      -0.897540332229183, -1.041943103505, 0.379425827374942, -1.00707718576756,      -0.889182530787803, 0.148432805676879, -0.287928359114935, -0.747152636892815,      -1.41003790431546, -0.611571256991109, -1.02569548477235, -1.02700056733181,      -1.45808867127733, -1.47973458605138, 2.23643966561508, 2.69397876103083,      0.81841473415516, 2.12167589051282, 0.267133799544379, -0.326215175076418,      -1.08788244901967, -1.18733017947214), X4 = c(0.050598970482242,      -0.0279694583060402, 0.999225143631274, 0.199872317584803, 0.779316284168575,      -0.3552692229881, -0.232161792808608, -0.333479851296274, -0.748169603107953,      -0.57785843363913, -0.480747933235349, -0.740466500603612, -0.618559437949564,      -0.591541699294345, -0.538855647639331, 0.431376763414175, -0.327931008191724,      -0.469416282917978, -0.659224551441466, -0.55285236403596, -0.637082867133913,      -0.780321541069982, -0.40539035027884, -0.54024676972473, -0.185562290173831,      0.054439450703482, 0.624097793456316, 0.24018937319873, -0.264194638773171,      -0.389590537012038, -0.42771343162755, -0.738790918078674, -0.122411831542788,      0.600119921164627, 0.0442597161778152, -0.0955011351192086, -0.521259643827527,      -0.550050365103255, -0.504566887441653, -0.506571005423286, 0.523650149759566,      0.341920916685254, -0.396343801985993, -0.366532239883921, -0.739276449002057,      -0.56054127343218, -0.587601788901296, -0.56798329186843, -0.454937006653748,      -0.672730639942183, -0.564864467446687, -0.678853515419629, 0.573072971483937,      0.596973680548765, 0.0403978228634349, 1.93617633381248, 2.54301964691615,      0.363075891004736, 0.0205658396444095, 0.560923287570261, 1.24212005971229,      2.32518793880728, 2.69979166871713, 0.626868716830008, 0.219463581391793,      0.236477261174534, -0.115429539698909, -0.49754151674106, -0.40827433350644,      -0.0433283703798658, -0.578451015506926, -0.714208713291922,      -0.802387726290423, -0.836794085697031, -0.471800405613954, -0.668030208971065,      -0.610945789491312, -0.780838257914176, -0.411360572155088, -0.494388869332376,      -0.63231547268951, -0.743853022088574, 4.90627675753856, 3.38455016460328,      0.859445571488139, 2.42212262705776, -0.324759764820016, -0.541581784452693,      -0.485324968098865, -0.770539730529603), X6 = c(0.0150287583709043,      -0.151984283645294, -0.347950002037732, 0.379891135882966, 0.129107019894704,      -0.314047917638528, -0.516381047940779, -0.751192211830495, -0.884460389494645,      -0.462363867892961, -0.397583161539858, -0.559528880497725, -0.842987555132397,      -0.922797893301111, -1.01175722882932, 0.32346425626624, -0.610909601293237,      -0.605155952259822, -1.29840867980623, 0.0793710626694382, -0.806959976634144,      -0.674523251142452, -0.960113466801064, -0.783836535852452, -0.0665321645536412,      0.482235339656537, -0.319499220427413, -0.115345965733089, -0.30806448545927,      0.251747727063608, -0.305013811851957, -0.931916656036151, 0.415032839884745,      0.337184728843034, 0.0584335852357015, -0.0712185313438638, 0.78632612201797,      0.490831043388539, 0.8902425262631, 0.160088439571744, 0.90343086944952,      0.928495373121098, -0.389259569427933, -0.304578433259833, -0.593364448723133,      -0.411333868741105, -0.882691663964141, -0.91208274239495, -0.708633954450382,      -0.339396626779965, -0.420927315080057, -0.421383857909298, 0.407474183771483,      0.629710767351175, -0.438726438495567, 2.40977730548689, 2.47250810430208,      0.783562677342961, 0.781304150896319, 0.563221804716475, 1.85514126067038,      1.30723846671955, 1.94869625911545, 0.876751836832149, 0.626629859119409,      0.067113945916172, 3.54280776301513, 0.0082773667305384, -0.311414481848668,      -0.732779325538588, -0.594477082903005, -1.0385239418576, -1.04141739541776,      -0.99472304141247, -0.599659297534257, -0.804801224448196, -1.13096932958525,      -0.641957537144073, -0.722959119516237, -0.671146043591047, -0.714432955420477,      -0.766750949574034, 1.20993739830475, 2.79011376379402, 2.64532317075082,      2.54251033029822, -0.539516582572252, -0.6419726544563, -0.663768795224503,      -0.644829826467113)), .Names = c("Y", "ID", "Sess", "X3", "X4",      "X6"), row.names = c(NA, -90L), class = "data.frame")


title: Interpreting average linkage clustering results in SPSS

text: I am trying to cluster 6 PC scores for a 1600 cases data (so its a 1600 x 6 matrix). I am using average linkage clustering technique that would be followed by a K mean clustering. The purpose of Average linkage clustering is to provide seed values for K mean clustering. Would somebody explain me how to extract the results of average linkage clustering from Agglomeration Schedule?


title: In what ways do data mining and machine learning competitions aid/detract from these academic fields and their commercial applications?

text: The Netflix competition kicked off a plethora of competitions to develop the best methods to learn information from given datasets. These include predicting traffic flow, predicting the success of grant applications, estimating the time patients spend in the hospital, and constructing good representations of data in an unsupervised manner. Some are run independently, others are hosted by platforms such as Kaggle and TunedIT. Clearly, these competitions can spur interest in the field and some specific applications. However, they may also cause a large amount of resources to be devoted to minutely optimizing a narrow application. To what extent are these competitions beneficial, and to what extent are they detrimental? How can competition organizers structure these competitions to maximize the potential academic and commercial benefits, and what common pitfalls should they avoid?


title: 3-class classification with different costs

text: Consider the following scenario: I have an auto dealership, and I want to decide whether to buy a used car. If it's a good buy, I would make a substantial amount of money of it. If it's a bad buy, I would lose a substantial amount of money. If it's neither ("nop") , I would lose a small amount of money. Naturally, I start with training and testing sets, and then use a simple 2-class classifier, that aims to predict whether it's a good buy or a bad buy ( I used gbm ,whose latest version, by the way, now includes also multinomial regression!). Once I have my classifier, I can look at my testing set, and inspect the cars that the classifier told me to buy. with these cars, I would calculate my expected future earnings:                  number of good buys X profit - number of bad buys x loss - number of 'nop' buys x smaller loss.      So, if you think about it, I actually shouldn't be trying to classify. I should try and find an area in my features space in which the given criterion is maximized. However, this sound like a combinatorial problem. So (finally), my question is: what algorithm / method is there that I can use to maximize my cost function? Thanks a lot!


title: Significant change in estimate with new observation?

text: Suppose I have some yearly observations (e.g. rates) $x_1$,$x_2$...,$x_{t-1}$. I estimate some parameter of interest as $\bar{x}_{t-1}=\frac{1}{t-1}\sum_i{x_i}$. Now at year $t$, I get a new observation $x_t$, so that I can compute $\bar{x}_{t}$.   1. How can I "decide" if I need to update my estimate from $\bar{x}_{t-1}$ to $\bar{x}_{t}$?   2. I get new observations for subsequent years, how can I "decide" to update my estimate or not?


title: How to assess effect of intervention in one state versus another using annual death rate data?

text: I am a beginner in statistics. I have these unpublished data (cases and deaths) of a disease for 7 years (2004-2010) from 2 neighbouring states. The study was started in 2004. State 1 and State 2 received different treatments. The death rate is high in one State. I want to prove that treatment given in State 1 is superior compared to that given in State 2. Time series was not accepted because of many possible confounding variables. I have SPSS and Comprehensive Meta-analysis software. Please guide me.   1. What is the best design? Prospective cohort study or Case-control study?   2. What is the best effect measure (Odds ratio, risk ratio, any other)?   3. What is the best statistical test (?Chi square etc)   4. Can I use meta-analysis in this case?   5. Any other information you think that would be of use.              State 1         State 2      Cases   Deaths  Cases   Deaths     2004    1125    5   2024    254     2005    1213    5   1978    209     2006    1003    4   2294    217     2007    1425    6   2312    249     2008    1172    4   1528    197     2009    1092    3   1683    204     2010    1316    4   2024    218


title: Predictive Variable Selection

text: You may find this a basic question. I am trying to model a time series model problem. My understanding is that the most parsimonious model (least variables) with the best result is best because most models don't deal well with a large number of irrelevant variables and that these variables will only introduce noise into the model, or cause it to over-fit. My question is when shall I decide to remove or not a variable? I have been using trial and error of the model and add a variable if the error gets smaller, however this is highly inefficient when dealing with a large amount of variables. I guess beyond this basic procedure some other methods must exist and be more efficient. I have tried from package leaps the function regsubsets,however seems not very effective as per the results. Would you happen to know of other packages or functions that may help do variable selection.


title: How to choose between Pearson and Spearman correlation?

text: How do I know when to choose between Spearman's $\rho$ and Pearson's $r$? My variable includes satisfaction and the scores were interpreted using the sum of the scores. However, these scores could also be ranked.


title: Post hoc comparison in two way ANOVA with covariate using R

text: Does any body know how to run a post hoc comparison in a 2X2 ANOVA with covariate in R. multcomp package seems very nice but I could not find clear answer or example to my question with this package. Thanks so much in advance


title: Estimating parameters of a normal distribution: median instead of mean?

text: The common approach for estimating the parameters of a normal distribution is to use the mean and the sample standard deviation / variance. However, if there are some outliers, the median and the median deviation from the median should be much more robust, right? On some data sets I tried, the normal distribution estimated by $\mathcal{N}(\text{median}(x), \text{median}|x - \text{median}(x)|)$ seems to produce a much better fit than the classic $\mathcal{N}(\hat\mu, \hat\sigma)$ using mean and RMS deviation. Is there any reason to _not_ use the median if you assume there are some outliers in the data set? Do you know some reference for this approach? A quick search on Google didn't find me useful results that discuss the benefits of using medians here (but obviously, "normal distribution parameter estimation median" is not a very specific set of search terms). The median deviation, is it biased? Should I multiply it with $\frac{n-1}{n}$ to reduce bias? Do you know similar robust parameter estimation approaches for other distributions such as Gamma distribution or the Exponentially modified Gaussian distribution (which needs Skewness in parameter estimation, and outliers really mess up this value)?


title: Meta-analysis in R, using metafor package

text: How must I syntax the `rma` function from metafor package in order to get results in the following real-life example of a small meta-analysis? (random- effect, summary statistic SMD)               study,    mean1,  sd1,    n1, mean2,   sd2,   n2          Foo2000,    0.78,   0.05,   20, 0.82,    0.07,  25          Sun2003,    0.74,   0.08,   30, 0.72,    0.05,  19          Pric2005,   0.75,   0.12,   20, 0.74,   0.09,  29          Rota2008,   0.62,  0.05,   24, 0.66,    0.03,  24          Pete2008,   0.68,  0.03,   10, 0.68,    0.02,  10      I tried metafor's tutorials and R function's examples with a poor outcome for the time being :(


title: Residualized Coefficient Using Solely Variances and Covariances

text: I just want to confirm that I'm on the right track. I need to derive an expression for b^Y1.2 that refers solely to variances and covariances (that is, only to bi-variate and mono-variate terms). I have come up with the following: b-hatY1.2 = (CovY1)/(Var1) - (Cov1,2/V1)*(b-hatY2.1) That should be correct, but I just wanted some feedback.


title: Curious Sample Correlation Property

text: I'm examining correlations in a data set with a large number of variables but small sample sizes. To get a feel for how these quantities behave, I generated some random data and looked at the distribution of correlations:               n = 4     y = matrix(rnorm(1000 * n), 1000, n)     x = matrix(rnorm(1000 * n), 1000, n)     p = as.numeric(cor(t(x),t(y)))     hist(p)      To my surprise, the distribution is almost perfectly uniform: ![histogram](http://i.stack.imgur.com/XbDeb.png) Does anyone have an explanation for this phenomenon? It makes some sense in that for n=2 we have either p=1 or p=-1, and as n->infinity the distribution becomes normal, so this distribution falls somewhere inbetween. But why uniform? I'm stumped.


title: How to handle non existent (not missing) data?

text: I've never really found any good text or examples on how to handle 'non- existent' data for inputs to any sort of classifier. I've read a lot on missing data but what can be done about data that cannot or doesn't exist in relation to multivariate inputs. I understand this is a very complex question and will vary depending on training methods used... Eg if trying to predict laptime for several runners with good accurate data. Amongst many inputs, possible variables amongst many are:   1. Input Variable - First time runner (Y/N)   2. Input Variable - Previous laptime ( 0 - 500 seconds)   3. Input Variable - Age   4. Input Variable - Height . . . many more Input variables etc & Output Predictor - Predicted Laptime (0 - 500 seconds) A 'missing variable' for '2.Previous laptime' could be computed several ways but '1. First time runner' would always equal N . But for 'NON EXISTENT DATA' for a first time runner (where '1. First time runner' = Y) what value/treatment should I give for '2. Previous laptime'? For example assigning '2. Previous laptime' as -99 or 0 can skew the distribution dramatically and make it look like a new runner has performed well. My current training methods have been using Logistic regression, SVM, NN & Decision trees


title: Difference in difference using spss

text: I have 2 groups of patients (control and treated) and 1 variable measured at two time moments (before and after treatment). Using paired T test I am able to determine the mean (with 95%CI) difference in each group (at baseline and after treatment). Using split-plot ANOVA I can say that the effect of the treatment is not different between the two groups (from the Test of Between Subjects Effects). But I am not able to determine the difference in difference (mean with 95%CI). I forgot to say that I am using SPSS. Thank you


title: Best estimator for small samples

text: I have a sample of an unknown distribution for which I would like an estimate of the true value and the expected error on that estimator. Since the underlying distribution is not necessarily Gaussian and the sample size is small I hesitate to use mean and standard deviation. What is a reasonable estimator and an upper bound of the error on the estimator? By upper bound I mean something like 95% CL but would accept anything justifiable. The error on the estimator will go into a fitting method so it is much more important to have an quality estimate of the error than an unbiased estimator of the value itself. Please include the range of sample sizes that this should be applied over before a more estimator like mean, median or mode is acceptable.


title: How to check if a group is similar to another, in two different classifications

text: Suppose I have two different classification, e.g.                    C    R     s1  c1   r2     s2  c1   r2     s3  c3   r2     ...      Where sn are the classified samples and the columns of the matrix represent the different classifications (one called `R` \- reference, the other `C` \- clustering). I want to check if, for example, c1 samples are significantly classified (read: the most part) as r2. If so, I say that r2 samples are _enriched_ in c1 (I do not know if it's the correct term). To do so I thought about checking if a r-class is enriched in an c-class. I thought about checking if r-class is more present in c-class samples with respect to all the other r-classes combined. I checked significance using fisher exact test on samples count. My contigency table is:               a = samples in ci AND in rj     b = samples in ci AND NOT in rj     c = samples not in ci AND in rj     d = samples not in ci AND NOT in rj                 rj   ~rj      ci     a     b     ~ci     c     d      Doing so, I expect to find if the distribution of samples is significantly different. If so, and if a and d are larger than b and c (read: samples count of rj is maximized in ci and samples count of the other rj-classes is maximized in the other ci-classes), I think I can assume that I can say rj is enriched in ci. Does this make sense? Thanks


title: Streaming k-means

text: I want to perform something like streaming/online/out-of-core kmeans clustering on large data. Here is simple idea:   1. Break all data into N chunks.   2. Read from disk 1st chunk and calculate centroids using original kmeans algorithm.   3. Use centroid of previous chunk to initialize next chunk centroids.   4. Do steps 2-3 N times. But how it will differ from original kmeans clustering by result? Maybe I can improve it somehow? Also I found mini batch kmeans but it seems it can't work with large data.


title: Regression with multiple dependent variables and 2 sets of multiple independent variables

text: In my data I want to be able to determine which set of IVs (there are 8 variables total, 4 in each set) does a better job of predicting a set a DVs (there are 4). I know I could use multivariate regression to use all DVs in the same analyses, but I don't know how to group the IVs into 2 sets. On the other hand I could look at incremental validity using linear regression and putting the different IV sets into different blocks, but in SPSS (maybe that is my problem...) you can only have one DV.


title: How should I choose WLS weights with categorical IVs?

text: I have a 2X2 design (2 categorical independent variables (IVs) with 2 levels each), and a single dependent variable (DV). I have 91 samples in total with these cell counts: 30, 19. 20, 22. My data for the DV is heteroscedastic across one of my IVs (the variance is 2.5 with N=50 at one level, and 1.2 with N=41 at the other level, Levene's Test sig. is 0.012). I'd like to run Weighted Least Squares regression to try to manage this difference, but I'm not sure how I should be estimating the weights. Intuitively, it seems like I might be able to compute the variance of each of the 4 cells in my design, and then for each cell, set the weight for every element in that cell to the inverse of the variance. Is that a reasonable thing to do, or should I try something else instead?


title: What's the difference between regression coefficients and partial regression coefficients?

text: I've read here that > When the independent variables are pairwise orthogonal, the effect of each > of them in the regression is assessed by computing the slope of the > regression between this independent variable and the dependent variable. In > this case, (i.e., orthogonality of the IV’s), the partial regression > coefficients are equal to the regression coefficients. In all other cases, > the regression coefficient will differ from the partial regression > coefficients. However, the document did not previously explain what the difference between these two types of regression coefficients is.


title: Noncentrality Parameter - what is it, what does it do, what would be a suggested value?

text: I have been trying to brush up on my stats knowledge, especially in relation to Sample size determination and Statistical Power Analysis. But it seems that the more I read the more I need to read. Anyway I found a tool called G*Power which seems to do everything I need but I am having an issue understanding Noncentrality Parameter, what is it, what does it do, what would be a suggested value etc? The information on wikipedia etc is either incomplete or I'm not doing a very good job understanding it. I am conducting a series of two tailed z-tests if that is any help. p.s. Can anyone add better tags to this question?


title: Does the mean and standard deviation need to be reported when writing up the simple effects analysis of a ANOVA interaction?

text: I'm currently carrying out a four-way mixed factorial ANOVA and in the process of writing up the results. I've done simple effects analysis on all the significant interactions between the factors, however I have noticed that the output doesn't include standard deviations. I've had a look at past research and it seems that sometimes the means and SD's are reported when writing up the post hoc tests but sometimes they aren't. This is how I have reported the findings so far in which I've just left gaps for where I thought the SD's should go (if they are needed) ... Simple effects analysis revealed that Corrugator activity was only significantly higher when facial stimuli was exposed for 56ms during the empathic narrative condition (M = 2.88, SD = …) compared to the control narrative condition (M = -1.4, SD = …), F(32) = 4.77, p = .04, η2 = 1.3 (see Figure 1) If I do need to include SD's how do I go about including them in the analysis as there is no option (that I am aware of) that allows SD's to be in the simple effects outputs? Will I have to do subsequent analysis to find out? Thank you for any help that you might be able to provide! Emma


title: On the stopping criterion of coordinate descent method for linear SVM with $\ell_1$-regularization

text: I am trying to implement the coordinate descent method to solve the dual of linear SVM problem, but blocked at the stopping criterion. The **dual** of linear SVM problem is: $$\min f(\mathbf{x}) = \frac{1}{2}\mathbf{x}^\top K \mathbf{x} - \mathbf{1}^\top\mathbf{x},\quad \mbox{s.t. } 0\le x_i\le C \quad i=1,\ldots,m$$ where $C$ is a positive constant and the matrix $K$ has the elements $K_{ij}=y_iy_jx_i^\top x_j$ where $y_i\in\\{-1,1\\}$ are the labels and $x_i\in\mathbb{R}^n$ are the data. At iteration $k$ we perform $m$ inner iterations where the $i$-th inner iteration updates $x_i$ by solving: $$x_i^{k+1} = \arg\min_{y} f(x_1^{k+1},x_2^{k+1},\ldots,x_{i-1}^{k+1},y,x_{i+1}^{k},\ldots,x_{m}^{k})$$ under the constraint $0\le y\le C$. My question is, what is the stopping criterion of this algorithm if we want to   1. Obtain an $\epsilon$-accurate solution $\mathbf{x}$, i.e. $f(\mathbf{x}) - f(\mathbf{x}^*) \le \epsilon$ where $\mathbf{x}^*$ is the true optimal solution? For the moment I take $f(\mathbf{x}^k) - f(\mathbf{x}^{k+1})<\epsilon$ but $f(\mathbf{x}^k) - f(\mathbf{x}^{k+1})$ is only the lower bound of $f(\mathbf{x}^k) - f(\mathbf{x}^*)$.   2. Obtain an $\epsilon$-accurate duality gap? **Update: for this problem, strong duality holds, thus an $\epsilon$-accurate solution yields an $\epsilon$-accurate duality gap.** Thank you in advance.


title: Determining significance and variance in A B C D testing

text: **Scenario:** An email campaign where four different email designs (treatments) are sent to four different populations of equal size in an attempt to find which performs better. The results returned are:               Treatment     Population        Clicks      A             120 (Group 1)     45      B             120 (Group 2)     35     C             120 (Group 3)     68     D             120 (Group 4)     52      **Goals:**   1. To figure out if any treatments perform (statistically) significantly better then others.   2. To find the amount of random variance expected in the test.  **Problem:** While I can easily see treatment C appears to perform better, I'm aware that there will be a number of clicks simply due to 'chance' (or factors not related to the treatment) in each test. Given the small population, my gut tells me if that if I could run the test again it's quite possible that the results could indicate another winner. Essentially I'm trying to determine how much variability I should expect and factor into the results to determine if the winner is actually a winner or just ahead due to chance. **Approaches:** I have only basic stats understand and have tried the following ideas:   * looked at approaching this through standard deviation, confidence limits and ANOVA, however given that there is only one data point per treatment these don't seem like the correct tools.   * thought about conducting an A|A|A|A test next time, where the same treatment is sent to different groups and see what the variance is, then use that as fuzz factor Any suggestions on the correct tool to apply to achieve what I'm after would be hugely appreciated!


title: Estimation of Density Function from a Transformed CDF

text: Suppose I can observe $x_1,...,x_n$ as the realization of the random variables $X_1,..,X_n$. Using $x_1,...,x_n$, I can estimate the empirical cumulative distribution function (CDF), $F_n(x)=\sum_{i=1}^n\frac{I(x_i\leq x)}{n}$. Now, with a given $\lambda$, I can transform this CDF by using the Wang transform which is $F^*(x)=\Phi\big[\Phi^{-1}(F_n(x))-\lambda\big]$, where $\Phi(.)$ is the cdf of standard normal distribution. **Question** : How can I estimate the $f^*(x)$ (i.e. the probability density function under the new transformation) using $F^*(x)$? Is there any package in R to do that?


title: Associating a probability to disease propagation in regions of a map

text: I am trying to model disease propagation. I am considering a map separated into regions. One region is infected. There is data over time, with the number of people travelling into and out of the infected region and which region they travel to (a table of transport of individuals over the whole map). This is an O-D matrix (wiki description). There is also some data about the prevalence of this disease in the region. The questions I would like to ask are:   1. What are some classical options for modelling disease propagation to different regions with this spatial data which will produce easy to interpret probabilities? (eg. SIS/SIR models which require simulations that I can just average over, or a time series model such as ARM with some measured parameters?)   2. What specifics of the disease must I know before hand? (eg. parameters of infection in different circumstances?)   3. What is the common name for these types of applications in research? (eg. temporal-spatial disease modelling? etc.)


title: Propagation of errors using simulation

text: If I construct confidence intervals for a calculation based on simulating inputs rather than using error propagation formulas, would this be considered as belonging to 'Monte Carlo' methods? I would have several inputs into the calculation and they would each have standard errors that I would use to generate random numbers from a distribution using R. I would use many inputs and then get the standard error of the outputs. Using the error propagation formulas would be possible, but it would be tedious and prone to mistakes in my opinion. Is this method usually acceptable and what name would I use to best describe the method?


title: Correcting for outliers in a running average

text: We have a daemon that reads in data from some sensors, and among the things it calculates (besides simply just reporting the state) is the average time it takes for the sensors to change from one value to another. It keeps a running average of 64 datapoints, and assumes that runtime is fairly constant. Unfortunately, as demonstrated by the below graph, the input data isn't the most pristine: ![](http://random-hall.mit.edu/blog/wp-content/uploads/2009/04/graphpng.png) (Each line represents a different set of data; the x-axis doesn't really mean anything besides a vague historical time axis). My obvious solution for dealing with this would be to create a histogram of the data and then pick the mode. However, I was wondering if there were other methods that would yield better performance or would be more suited for operation with a running average. Some quick Wikipedia searches suggest algorithms for detecting outliers may be also suitable. Simplicity is a plus, since the daemon is written in C. **Edit** : I scoped out Wikipedia and came up with these various techniques:   * Chauvenet's criterion: using the mean and standard deviation, calculate the probability a particular datapoint would happen, and then exclude it if the probability is actually that bad is less than 50%. While this seems to be well suited for correcting a running average on the fly, I'm not quite convinced of its efficacy: it seems with large data-sets it would not want to discard datapoints.   * Grubbs' test: Another method that uses difference from the mean to standard deviation, and has some expression for when the hypothesis of "no outliers" is rejected   * Cook's distance: Measures the influence a datapoint has on a least squares regression; our application would probably reject it if it exceeded 1   * Truncated mean: Discard the low end and the high end, and then take the mean as normal Anyone have any specific experience and can comment on these statistical techniques? Also, some comments about the physical situation: we're measuring the average time until completion of a mechanical washing machine, so its runtime should be fairly constant. I'm not sure if it actually has a normal distribution. **Edit 2** : Another interesting question: when the daemon is bootstrapping, as in, doesn't have any previous data to analyze, how should it deal with incoming data? Simply not do any outlier pruning? **Edit 3** : One more thing... if the hardware does change such that the runtimes do become different, is it worth it to make the algorithm sufficiently robust such that it won't discard these new runtimes, I should I just remember to flush the cache when that happens?


title: Can someone give a clear-cut idea of $E(X|X<Y)$?

text: If $X$ and $Y$ be two i.i.d. random variables, then what should $E(X|X<Y)$ essentially look like? P.S.-Is the denominator equal to $\frac12$? (the two r.v.'s being iid serving as the motivation behind this.)


title: Discrepancy between calculated and simulated error propagation

text: Just curious if anyone knows why there is a tiny but definite difference between the standard deviation as calculated by the error propagation formula and the simulation. It is only a 0.2% difference and doesn't cause me any problems personally but I am curious which of the two would be considered more accurate? Adding another zero to the rnorm function causes my computer to crash. But I repeated this code 10 times or so (very small variance) and I am confident that the mean is in fact 0.2% higher than it should be.               # This takes about 6 seconds to run     x <- rnorm(10000000,mean = 100, sd = 10)     y <- rnorm(10000000,mean = 100, sd = 10)     z <- x*y     # SD of x*y using simulation     sd(z)     # SD  of x*y using formula     (((10/100)^2+(10/100)^2)^.5)*(100*100)


title: Statistical test for a class versus multiple classes

text: We are given n classes of real values. The number of elements in each class is not too high and can be different between classes. We want to test whether there exist a class with significantly higher (average) values from all other classes or not. The case is somehow similar to ANOVA, but slightly different. What test do you propose?


title: Calculating mutual Information

text: I am making a spam filter for an assignment. One of the steps outlined in a paper that I have read indicates finding the mutual information for all the words in the text and then choosing the 500 or so with the highest MI. The equation given in the paper (and elsewhere) is: $MI(X_i, C) = \sum{P(X_i, C) \log{\frac{P(X_i, C)}{P(X_i)\cdot P(C)}}}$ where $X_i$ is a feature and $C$ is a class. It is not clear in the paper but $X_i$ will be a word and $C$ will be either spam or not-spam. Can I just be clear of the calculation that I should be doing? I think it is as follows for the word "viagra": $P(X_i, C)$ in this case will be either the frequency that the word "viagra" appears in all of the spam emails divided by the number of words in the spam emails, or the same for the ham (non-spam) emails. $P(X_i)$ is the probability of the word occurring at all in both spam and ham emails. $P(C)$ is the prior for a class. So basically we are doing two calculations (one for the spam and one for the ham) and adding them together. Is this right do you think? My understanding is that mutual information is like finding the information gain for a feature. Is this right? Hope this makes sense. Many thanks in advance.


title: Efficient solution of fmincg without providing gradient?

text: I'm working on multiclass logistic regression model with a large number of features (`numFeatures`>100). Using a maximum likelihood estimation-based cost function and gradient, the `fmincg` algorithm solves the problem quickly... However, I am also experimenting with a different cost function and do not have a gradient... Is there a good way to speed up the calculation process? E.g. a different algorithm or `fmincg` setting?


title: Omitted dummy variable coefficient in OLS

text: I am running an OLS regression using dummy variables built from categorical variables. Say, race became race1, race2 and race3. I omit race1 in order to escape the dummy variable trap and run OLS and gain some coefficients for race2 and race3. What would be the coefficient for race1 in the regression equation? Or I just should not include the omitted dummy variable into the final equation to estimate the predicted value?


title: Are bias weights essential in the output layer, if one wants a universal function approximator (or non-linearly separable problem solver)?

text: I am learning about ELM (Extreme Learning Machines) and **it appears** to have no bias weights at the output layer. Besides that, just to clarify, the kind of ELM I am refering to are topologically no different than 2-layer MLPs (or 3-layer if you count the input as a layer). So, I am wondering whether I understood wrong about the absence of output layer bias weights in ELMs or an MLP can "survive" without it.


title: Creating alert system based on standard deviation and mean: How to deal with small sample size, trends and seasonality?

text: I'm trying to build a daily alert system to let me know when something unusual has happened in my analytics data that might require further investigation. Right now I'm taking 7 weeks worth of data, creating an average for each day of the week (ie Average Conversion Rate for Monday). I then compare how much this Monday has changed from the regular Monday average in terms of standard deviations. If it's +/-2 stdevs, then it warrants investigation. A few concerns I have :   * How many weeks should I take into account? My concern is that if I take too many then I might have a trend emerging that skews my alerts. If I take too few, I might have high variance and the system might fire alerts needlessly.    * Would it be better to take a rolling average for say last 7 Mondays and compare this Monday to that rather than the mean?   * What would the best way be to deal with seasonality - ie bank holiday Mondays which happen occasionally but would skew averages and standard deviations.   * Some websites might have little traffic per day. How should I deal with them? Would I have to take a longer time frame to calculate averages? How would I decide how many to take?


title: How to draw samples from a Bayesian nonparamatric density estimation? [DPpackage]

text: I am trying to compute a Kernel Density from high dimensional data ($n > 2$). The underlying (generative) model is assumed unknown. The goal is to draw samples from this estimate, in a sense creating synthetic data akin to the underlying distribution. I was hoping to use Bayesian nonparametrics for this, since the complexity of the model can then grow with the amount of observed data. The intuition was to use the Dirichlet Process with a Gaussian kernel. DPpackage supports this functionality, but I have little clue on how to interpret the results or sample from them. Does anyone have a short explanation on how to interpret the returned values from the `DPdensity` function, and possibly sample from them?


title: How similar are my 2 data sets?

text: I am kind of stuck with an easy question: I have two data sets with experimental data. The data sets do not have the same size. I would like to show that these data sets are possibly coming from the same experiment. I tried a two-sample $t$-test; it shows that the data are significantly different. Is there a way to generate something like a $p$-value for similarity instead of difference? Update:   Here an example:   Date set 1 (Vector): 1 1 2 3 1 2 1 3 4 1 Mean: 1.9   Data set 2 (Vector): 2 2 1 2 2 1 1 2 2 3 2 2 Mean: 1.83 How would you now show that this data sets are from one experiment?


title: binary repeated measures with 3 groups and 3 condititions

text: I'm using SPSS and I'm having touble with deciding what statistical method should be used: there are 3 groups with different participants and treatments and they were measured in three different experimental conditions, but one of the variables is binary (another one is continuous, derived data). Should I use glmm or generalized estimating equations (or something entirely different)? Would it be better to analize the binary variable together with the other one or separately? Thanks in advance and sorry for my english!


title: Normalize binary features for logistic regression

text: I am applying logistic regression and I am using a mixture of continuous and binary features. My question would be how do you go about binary features when trying to normalize. I would very much appreciate any help.


title: Combining observed Gaussian error with common fractional model error

text: I'm currently trying to incorporate model errors into a likelihood function to fit the model to some results. Let's just call them $x_{m,i}$ for the model and $x_{o,i}$ for the $N$ observations, indexed by $i$. Each observation has a given Guassian/normally-distributed error $\sigma_{o,i}$. For the model, I know that, in any given model, all the $x_{m,i}$ have the _same_ fractional error. That is, a single fractional error $\sigma_f$ is (for now...) Gaussian/normally-distributed, and each observation has error $\sigma_f\times x_{i,m}$. How do I incorporate this into a likelihood function for fitting the model? I'm not sure how much detail is needed for this to make sense, but basically I'm trying to do a $\chi^2$-minimization with $$\chi^2=\sum_{i=1}^N\frac{(x_{o,i}-x_{m,i})^2}{\sigma_{o,i}^2+\sigma_{m,i}^2}$$ but I'm not sure what to do for $\sigma_{m,i}$. My instinct is that just putting in $\sigma_{m,i}=\sigma_fx_{m,i}$ is wrong, because that still allows the model errors to be, in some sense, drawn independently, even though there's only really one random process for the fractional model error, rather than an independent process for each observation. I tried to formulate something for $\sigma_{m,i}$ with propagation of error by writing $$\sigma_{\text{total},i}^2=\left(\frac{\partial(x_{o,i}-x_{m,i})}{\partial x_{o,i}}\right)^2\sigma_{o,i}+\left(\frac{\partial(x_{o,i}-x_{o,i}\frac{x_{m,i}}{x_{o,i}})}{\partial \frac{x_{m,i}}{x_{o,i}}}\right)^2\sigma_{f}^2 \\\=\sigma_{o,i}^2+x_{o,i}^2\sigma_f^2$$ but, as you can see, that just got me back to this answer I don't believe. I suspect a better way of thinking about it is in terms of a perfectly correlated error component. Well, this is what I'll think about next, anyway! Thanks in advance for any advice. I'm not a statistician, nor particularly well-versed in statistics, so I apologize if my description is incomprehensible to an expert. I'll try to clarify if necessary.


title: ANOVA type for dependent samples through time

text: The problem is how best to analyse data which involves multiple dependent measures, taken on multiple/repeated occasions and determine whether there's a difference between two groups. I'm trying to analyse data from a study involving cochlear implants. Now a cochlear implant has an array with multiple electrodes, and for each electrode it's possible to measure various different parameters. These parameters can be, and are, typically measured over time for the purpose of understanding how the parameter changes through time. Per good study design, in the problem I'm considering there's the control group and the experimental group. To note is that there's no cross-over element in the design. To be specific:   * Two groups of N test subjects: control and experimental   * 12 different values for any "measurement set", from a single person at a single point in time   * Multiple measurement points (in time), the repeated measures part   * Does the control and experimental group differ, generally?   * Is there a significant difference at the different time points?   * Are there any differences between the 12 values through time? Thus my questions:   1. My thought is that this is a repeated measures MANOVA analysis type. Have I understood the situation correctly?   2. How should the measurements at multiple time points be taken into account? Or maybe this is catered for with the (multiple) repeated nature of the analysis. Thanks in advance!


title: Generalizing results - what is the population?

text: Take as an example the well-known Boston housing dataset. In the original paper, the regression model was used to model the amount people were willing to pay for improved air quality. But the 'people' whose behaviour was being modelled was never specified. Presumably, it was not limited to 'residents of Greater Boston'. Only implausibly would it be 'humans anywhere on the planet'. Nor, 'all US residents'. So what inference are we invited by the authors to draw about the scope they intend for their theory? Perhaps 'all residents of large US metro areas'. But - what about other First World metro areas? The first two sections of the Boston paper scrupulously avoid any reference which might suggest that they had a particular 'population' to which their 'sample' was to apply. There is also the dimension of time - the data came from the 70s, whereas the theory is has no time limitations spelt out. Plus - going back to first principles - how could the Boston data be treated as a sample of a population including, say, the Chicago metro area when no data from the Chicago area was considered? Besides, lots of the covariates in the Boston regression weren't samples at all, but totals - so, for example, the proportion of blacks was supplied from census stats. My question: how can we conceive the Boston data as supporting the geographically unlimited propositions in the paper? Is there a mathematical way of showing equivalence or approximation of the Boston data analysis to other cities? Or is it just assumed as a matter of common sense - such that the paper's authors would be amazed to think anyone could have any doubts on the matter?


title: What is the best strategy to train and validate classification using PLS-[classifier] in caret package?

text: I have 4 clusters (see plot below) extracted from data of medical samples `N=218` measured for 11 genes/predictors `P=11` by this method: first PCA analysis validated to have 2 important PCs that explained 75% of the data, then different clustering algorithms, distances, linkages (in hierarchical approach only) were compared: the majority support the presence of 4 distinct clusters. Taking the scientific hypothesis into consideration, clusters out of $K$-means algorithm were found the most plausible and were the most balanced clusters too: class #1 `n=12`, class #2 `n=21`, class #3 `n=79`, and class #4 `n=106`. Projecting the observations on plane 1-2 component scores, revealed the below scatter plot with each cluster color coded. ![enter image description here](http://i.stack.imgur.com/09sTF.png) **The aim is to find a global optimum classifier using R after doing PLS to the data.** Knowing that these 4 clusters were actually the product of latent PCA components, it was natural to think of PCR as a next step to predict classes, but that approach turned out to be sub-optimal for two reasons: first, results do not related to probibilities (0-1), second, it does not relate well with the classes as the outcome variable. As many know, this would be better solved with PLS-DA method + softmax to find probabilities of class (0-1). However, many reports confirm the superiority of using LDA as a second step using the _scores_ of PLS, given that same standardization parameters (mean, sd) be used of the training set on the holdout-test set, even using the PLS projections out of the training set on the test set in order to get the _scores_ which would be the _actual_ holdout-test set to validate the classifier in question. From the methodology point of view, this path is potentially encompassed with many dangers and subtle errors when one is un/misinformed about the tools used in context. The `caret` package which is unique of its kind given the consistent infrastructure it provides to train and validate an array of different models making use of _de facto_ standard respective R-packages, and hence `caret` promotes itself as a road map to a validated modeling leveraging off R rich libraries. As heart to blood vessels, so `caret` to other packages in my opinion. That being said, unwatchful playing with the heart could cost you dearly, and might lead also to a stand-still or a _model-arrest_ of your data. R is free, many free books out there, but buying `caret` only book paid off, i.e. `Applied Predictive Modeling`. The help files, companion website (very appealing btw), are great resources but they won't substitute the text inside the book IMHO. However, in the book, I couldn't find a direct answer to the PLS-Classifier two step method amid others. The potential with `caret` is immense, thanks to Max Kuhn and his colleagues, that primarily encouraged me to post this question. Back to the example above and the methodology of wish:   **Data splitting:** Training set (77% `n=168`) for 10K-cross-validation: tuning (model-specific parameters, feature selection `P=11`, and cost to deal with imbalanced clusters). For CV this would be roughly `n=150` for fitting the model using differnt parameters of wish and 'n=17`for evaluation of parameters (I would call the`n=17' the CV-test to avoid confusion later on). Repetition = 5, so this will make 10 folds x 5 times = 50 training folds (`n=150` each) and 50 CV-test folds (`n=17` each). Holdout-test set (23% `n=50`). **Q1** I know that one can do parameters' tuning along with feature selection at the same time (i.e., parallel), but how to evaluate the cost/weights if one would like to evaluate cost-sensitive models (SVM, CART, C5.0) using the PLS scores to counteract class imbalance? **Q2** What is the alternative approach when reserving a separate data set for cost evaluation (i.e., `evaluation set`), as recommended, is not possible given the small sample size in this case? can one do tuning of model parameter, feature selection, and cost for imbalance all three at the same time? if not what is the best practice in this case? **Q3** Given the small sample size, is bootstrapping preferred to CV? if yes how would it be implemented to do exhaustive tuning like above for the PLS scores? **Q4** Given the imbalance above, is there a way to ensure that each CV training fold would include the minimum number of `hard` class(es) in order to have good estimation on the CV-test fold? is there any argument to pass to ensure presence of the small classes each time fold would be generated? **PLS special notes** This is the approach in my mind (please correct me if I am missing something somewhere during the course): In each CV iteration on the many CV-traning folds, there should be a unique PLS projection matrix for each iteration that would be used in the next second step of getting PLS scores for the the respective CV-test set inheriting the same standardization parameters (`mean`, `sd`), this means that two things would be inherited; the PLS projection and the standardization parameters (mean and sd) in order to apply them to the CV-test folds, this way, given the example here, 50 values would be returned hoping to reach the best parameter in question. One complication though, there should be an argument to specify the desired number of PLS components to retain and to be used in calculation of scores out of each CV-test fold (better to be pre-defined in a previous tuning step may be). My expectation, is that after deciding on the best model, there should be a way to get the PLS projection matrix for the whole training set (i.e.`n=168`) along with (`mean`, `sd`) to apply them on the holdout-test to validate the best model. So in total, there would be 50 different PLS projection matrices, means, sds from CV step and 1 extra frothe whole training set, am I right?   Feature selection in this method would entail two things: predictors space and the PLS components space. **Q5** How to perform these two selections (predictor and PLS component) in `caret`? this is because feature selection here is different than otherwise since here we deal with scores rather than the observations themselves to determine best predictors that to construct the PLS components. **Note:** When one is happy with the best final model, it would be recommended to fit the model on the whole data set `n=218` to get the correct estimates withe the least uncertainty. A similar procedure is implemented in `caret::train()` function that can be fed with `preProc` argument to specify the type of desired pre-processing of data (most are mentioned in the help system but I couldn't find **PLS** among them, better if with an argument to specify the desired components similar to PCA pre-processing). I am aware of the fact, that inheriting pre-processing parameters to holdout test set and to CV-test, can only be performed using the `predict.train()` function, as opposed to calling the generic `predict()` function to the `$finalModel` that won't inherit pre-processing parameters. **Q6** How to implement this strategy (if correctly described) to train and validate the two-step PLS-[classifier] methodology using PLS scores subspace instead of observations making use of `caret` infrastructure? Thanks in advance.


title: How to calculate variance over a distribution of time duration proportions

text: I have events that take a proportion of an observed time slice. The observed variable is the number of events happening at any one moment of time. Therefore the mean is merely the sum of proportions, the minimum is the count at a time where there are the least overlapping. ![time series diagram](http://i.stack.imgur.com/DamQL.png) In the above diagram we have data that is in terms of start and end times for each 'instance'. We also have the proportion of overlap with the blue box which are the numbers in the table below the diagram. However what we are really trying to calculate are descriptive statistics on the random variable which is the number of instances present at any given time. My question is, how would I calculate the variation of this kind of dataset? * * * Thinking about this question a bit more, I realize that what I'm actually asking for is help in identifying the distribution that my data falls in. A better visualisation of the distribution above is the following diagram: ![diagram 2](http://i.stack.imgur.com/c6MoR.png) The data is generated by counting for each moment in time the events (above called instances) that are in progress. These start and end randomly. Note that despite appearances the above is not a histogram as the time variable is continuous. * * * Doing some more reading around, it seems that my data could be easily characterised as a Birth-Death Process as the event count goes up or down depending on events beginning and ending. I still don't see how to calculate the variance however.


title: Analysis of likert scale data for a radiology study on dose reduction techniques

text: To all, I have performed a study where 50 patients underwent a computerized tomography (CT) scan. Each scan was reconstructed using 3 different algorithms (no noise reduction, medium noise reduction, high noise reduction). These reconstruction algorithms create noise reduction which is an issue when trying to reduce radiation dose. Since all 50 patients had 3 reconstruction, there were 150 studies total. These studies were then anonymized and randomized. 8 separate radiologists then read each study and graded them on a Likert scale of 1-5 (1-awful, 5-ideal) for numerous qualitative categories including diagnostic confidence, perceived image noise, perceived contrast enhancement, and perceived image artifacts. I want to compare the readers scores between the different reconstructions. For instance, does the high noise reduction algorithm lead to improved diagnostic confidence compared to the other reconstructions. Also, in general, most studies were subjectively graded as either a 4 or 5, leading to negatively skewed data. So long story short, what test would you recommend for analyzing this data. I am thinking a Friedman test. Is this correct? If so, would you run a post-test of any kind. Any help would be great appreciated. Thanks, SK


title: Should I be using a Welch or standard T-test?

text: I'm really confused here. I am trying to compare two sets of data which differ in sample size. One has 257 entries and the other has 598. They are measures of scores on two different assessment types where the only possible scores are 1 through 4, with half points allowed. The means are a half point different, the standard deviations are .2 different, and the skews are different. I'm positive that the difference is significant. What test can I use to prove that the difference is significant? I think it's a 2 tailed, independent t-test, but apparently a Welch's test is used when the sample sizes are different? The answer to this is probably very simple, but I've spent enough time trying to figure this out that I'm more confused than when I started.


title: Comparing nested, non-linear models

text: I would like to compare the fit of two non-linear regression models: 1) $$ Y = (\Pi^{10}_{i=1}\beta_i^{x_i})^{1/\Sigma \beta} $$ 2) $$ Y = \begin{cases} (\Pi^{10}_{i=1}\alpha_i \beta_i^{x_i})^{1/\Sigma \beta} & \text{for a predefined subset of datapoints} \\\ (\Pi^{10}_{i=1}(1-\alpha_i)\beta_i^{x_i})^{1/\Sigma \beta} & \text{otherwise} \end{cases} $$ Can I still use the F-statistic to compare the fit of these two models?


title: SVM for unbalanced data

text: I want to attempt to use Support Vector Machines (SVMs) on my dataset. Before I attempt the problem though, I was warned that SVMs dont perform well on extremely unbalanced data. In my case, I can have as much as 95-98% -1's and 2-5% 1's. I tried to find resources which talked about using SVMs on sparse/unbalanced data, but all I could find was 'sparseSVMs' (which use a small amount of support vectors). I was hoping someone could briefly explain: 1) how well SVM would be expected to do with such a dataset   2) which, if any, modifications must be done to the SVM algorithm 3) what resources/papers discuss this Thanks!


title: Why we shouldn't be obsessed with unbiasedness

text: In my Bayesian statistics class, my professor makes the remark that we should not be obsessed with unbiased estimator. First: I understand this statement in the sense of trading biasedness for smaller variance, i.e. "good" estimator = estimator with smallest Mean Squared Error (MSE). This idea is seen in Frequentist also. But second, and the confusing part: My professor then said something about an unbiased estimator is picking up noise and not signal. Since most of the observed data are just noise, we want to shrink our estimate to 0. I do not quite understand what is noise and signal: 1) Don't all data points come from the same distribution, with the difference due only to sampling variation? 2) Why can't unbiased estimator (say, sample mean of a normal model) differentiate noise and signal? If it can't, what sort of estimator can? An example would be great, and a social science example would be best!


title: How to subset alternatives in nested multinomial logistic regression?

text: I am trying to predict whether or not captains in a particular groundfish fishery choose to fish on any given day and what variables may influence that decision. Originally I had planned on using vanilla logistic regression, where each captain decided to fish or not fish on any given day of the season in a particular fishery. However, this is not a realistic model as many vessels participate in multiple fisheries. Thus, I want to run a nested multinational logistic regression model with the structure at the bottom. Note that Group 1 and 2 exist to avoid violating the IIA assumption. Because each species has a different season and within each species captains have variable amounts of fishing quota lbs (when the quota is gone you can no longer fish even if the season is open), not all nodes are available at the same time. I would like to include a dummy variable for whether or not an individual node is an option, but I am not sure if that makes sense or even if this method (multinomial logistic) is still viable. I could try something like this Coding of semi-numerical variables i.e whuber's response, but I have a lot more subsets and nests within those subsets. Edit* I have since looked into, but at the moment have not implemented, a varying choice logit model. The description of the model can be found http://web.mit.edu/teppei/www/research/dchoice.pdf. I am using R and was planning on using lmer4, but I am not sure if I will be able to. ![enter image description here](http://i.stack.imgur.com/4kC9H.jpg)


title: two component non standard mixture (normal + unknown)

text: I have some univariate data which might be well modeled as a two component mixture where the first component is normal with unknown mean and variance and the second is some unspecified continuous distribution. The mean of the second distribution is greater than the mean of the normal. The goal is to classify the observations into each of these two distributions. Does this suggest a method to anyone? All suggestions much appreciated.


title: Convergence of empirical distribution parameters of a sequence of generated normal variables

text: I am generating a sequence of normal random variables (using the routines from boost C++ library). How fast would you expect the mean and the variance of the sequence converge to the actual variance? In my case the true values are:               mean=-0.01208333333333333     variance=0.0075      For `10000000` generated variables I get:               Mean=-0.0120757,      var=0.00750125      which intuitively seem quite far (especially variance) for this ridiculously large number of samples. Should I expect faster convergence (i.e. there is some issue with my code), or is it what you would expect?


title: Comparing fuzzy RD estimates to an "OLS" analogue

text: In IV, I've sometimes seen the estimates compared to OLS to give a sense of how the LATE may compare to ATE. What would the analogue be in fuzzy RD? What is the benchmark estimate? I'd imagine the analogue would be the structural equation without instrumenting for the endogenous variable? Would you include the RD polynomial in that? Would you limit it to a narrow bandwidth? Does the following make sense? Suppose the RD is for a test-score cutoff for a acceptance to an academic program, and the fuzzy RD looks at the impact of being apart of that program on future earnings. The OLS baseline would be regress wages on a dummy for participation in the academic program and the level of the test score? The instrument is not The IV baseline would be regress wages on participation in the program and the level of the test score, and instrument participation with a first stage that regresses participation on a dummy for whether the test score was above the threshold plus the level of the test score? Then, the RD nails that by doing 2 things: First, it replaces the level of the test score with a range of flexible polynomials in the level of test score: piecewise linear, quadratic etc ... Second, it limits the sample to narrower and narrower bandwidths around the discontinuity. Presenting estimates of numerous permutations between bandwidth and polynomial order? Is that a sensible framework for comparing ATE, LATE, and Weighted LATE (ie fuzzy RD)?


title: Error term/Innovation process in ARCH/GARCH processes?

text: I am wondering about the distribution of the error term/innovation process in a ARCH/GARCH process and its implementation, I am not sure about some points. The basic assumption is $r_t=\sigma_t*\epsilon_t$ where the $\sigma_t$ is the volatility, modeled by ARCH/GARCH and the $\epsilon_t$ are mostly assumed to be N(0,1). Now my questions are:   1. More sophisticated models drop this assumption. So I can say, e.g. $\epsilon_t$ follows a generalized hyperbolic distribution. So the mean does not need to be zero and the variance does not need to be equal to 1. This is correct, right?   2. If I use the rugarch package: It supports different distributional assumptions. But I am not getting the following: So they also drop the assumption of mean zero and variance one? Or are they using something like a "standardized" version?   3. Suppose I want to fit a GARCH(1,1) assuming, that the $\epsilon_t$ follow a generalized hyperbolic distribution, but the mean does not have to be zero and the variance does not need to be one. Is rugarch doing a jointly parameter estimation? So in my final output, do I get the parameters of the GARCH process and the parameters of my generalized hyperbolic distribution? My last question is, how can I implement this? I guess I have to use the following command:               ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1),      submodel = NULL, external.regressors = NULL, variance.targeting = FALSE),      mean.model = list(armaOrder = c(1, 1), include.mean = TRUE, archm = FALSE,      archpow = 1, arfima = FALSE, external.regressors = NULL, archex = FALSE),      distribution.model = "norm", start.pars = list(), fixed.pars = list(), ...)      the distribution.model has to be set to `ghyp`. Is this assuming a mean of zero and a variance of one? I think no, right? How can I use the hyperbolic distribution for distribution.model?


title: $X \sim $ Geometric$(p)$, $Y \sim$ NegativeBinomial$(2,p)$: find the MVUE for $p(1-p)$

text: Let $X$ and $Y$ be independent random variables with pdf \begin{align} f(x)&=p(1-p)^x,x=0,1,2,\ldots \\\ f(y)&=(1+y)p^2(1-p)^y, y=0,1,2,\ldots \end{align} Find the MVUE for $$ p(1-p). $$ I know $X \sim $ Geometric$(p)$ and $Y \sim$ NegativeBinomial$(2,p)$.   So the complete sufficient statistic is $$ T=(X,Y). $$ Let $$ g=\frac12 X\cdot 1_{\\{Y=1\\}}, $$ then $$ E(g) = p(1-p), $$ and by Lehmann-Scheffe $$ E(g|T)=g, $$ is the MVUE. Is this correct?   I find strange that $g$ doesn't change in $E(g|T)$.


title: Quantifying the overlap in sort task results with cells < 5

text: We have data from two sort tasks where the same items (N=87) were assigned to different types of group ('g' & 'd'). We want to compare the overlap between the assignments to the 'g' groups (N=13) to those of the 'd' groups (N=16). To give you an idea of the data I have created a truncated frequency table _(apologies for the table format - not sure how best to show it on here)_ : In this example 9 items placed into group 'g4' were also placed into 'd1'.                 | g1 | g2 | g3 | g4 | g5 |     ----------------------------     d1|  1 |  0 |  0 |  9 |  0 |     ----------------------------     d2|  0 |  0 |  0 |  0 |  4 |     ----------------------------     d3|  0 |  7 |  4 |  0 |  2 |     ----------------------------     d4|  0 |  0 |  8 |  0 |  4 |     ----------------------------      The size of the 'g' groups differ and so the potential number of matching items assigned to the 'd' groups depends on the size of the 'g' group. For instance, g4 actually has 18 items assigned to it but only 9 match with d1, g3 has 3 items assigned and none match d1. So the frequencies alone are not particularly representative. We wanted to have 1) a measure of the overlap between all of the 'g' groups with all of the 'd' groups taking into account the possible overlap, and 2). a global measure of the association between 'g' groupings and 'd' groupings. I am assuming (please correct me if I am wrong) we can't use Chi-Squared as most of the cells are less than 5 and we can't use Fisher's exact test as it isn't a 2x2 table. **My questions:** 1) Are there any other test we can perform to quantify the relationships between these assignments? 2) When I was playing with the data and ran a chi-squared test in R it gave me **Pearson's Residual values** for each GxD combination - are these residual values OK to use independently of the Chi-Square test? -- _They seem to represent what we need for our first requirement, i.e. providing an index of the overlap between groups taking into account group size - but I don't want to use these incorrectly._ 3). I saw the option of computing p-values by Monte Carlo Simulation in the chisq.test function in R. Could this simulation be used to overcome the problem of cells < 5 with the data?


title: Removing economic effects from time series

text: I am analyzing the number of donors cancelling commitments as a monthly time series which varies significantly with economic indicators, some according to internal data, and maybe according to season. I want to see the variance after adjusting for the economy and season (assuming the remainder of variance is due to the non-profit's behavior). Is a good way to do this to fit a linear regression model in R pseudo-code:               fit <- lm(cancels ~ economy + season)      And then use the residuals?


title: Regression with t distributed autocorrelated errors

text: I am new, so please bear with me, until its correctly formulated:   **In short:**   I am doing simultaneous non-linear regression (parameter estimation) of two different parametric models to three different, but related sets of measurements in an experiment. Some of the parameters in my parameter vector, exist in both models. _Current method (Coneptually):_   * I assume normally distributed errors.   * I estimated the weight (the standard deviation of error destribution).   * With the objective function the weighted sum of squared errors is calculated.   * My objective function calculates the log-likelihood from the weighted sum of squared errors.   * For all individual regressions, the log-likelihood is calculaan add the  _Current method (Mathematically)_ Each statistical model in the regression looks like this: $ L(p,\phi,\sigma^2) = − /2* (^2)+ 0.5∗(1−\phi^2 )− 1/(2^2 ) * _1 (,\phi)$, with $ S_1(p,ϕ) = (1−\phi^2 ) − (_1−(_;))^2+S_2(p,ϕ) $ $ S_2(p,ϕ) = ∑_{=2}^ (_ - \phi _{−1} −(_;)+\phi (_{−1} ;))^2 $ $\phi$ - parameter for autocorrelation   $\sigma$ - parameter for standard deviation   $p$ - parameter vector   $y$ - observations   $(_;)$ - model response **Questions**   1. I assume autocorrelated errors, but I need to take cross correlation into account, too. How does the statistical model look like to account for cross correlation?    2. How does my statistical model look like, if I use t-distribution or cauchy distribution instead. These have the possibility to describe fat tailed errors. A reference to a publication or textbook would suffice. Because sadly, I can't derive it myself.


title: How to interpret output from least trimmed squares estimate and compare it to OLS?

text: I have to compute and compare the least squares method on a model to the least trimmed method. The LS model results were:               Coefficients:                    Estimate Std. Error t value Pr(>|t|)         (Intercept)    -39.9197    11.8960  -3.356  0.00375 **      Air.Flow         0.7156     0.1349   5.307  5.8e-05 ***     Water.Temp       1.2953     0.3680   3.520  0.00263 **      Acid.Conc.      -0.1521     0.1563  -0.973  0.34405          The LTS results were:                 (Intercept)     Air.Flow      Water.Temp      Acid.Conc.      -3.580556e+01  7.500000e-01    3.333333e-01    3.489094e-17       How do I interpret the results for LTS? I know `Air.Flow` and `Water.Temp` are significant and `ACid.conc` is not. But don't know about LTS coefficients since there are no standard errors.


title: Statistical comparison of two means with a range not starting at 0

text: I don't know the exact term for this, so googling didn't work. I will explain exactly what I need below. I want to compare two sets of values such as blood sugar or blood pressure, where the values never start at 0. Please consider the example below: **Update** I compared the duration of motor block after spinal anesthesia. The results are not distributed normally: Min. 1st Qu. Median Mean 3rd Qu. Max. 75.0 140.0 160.0 157.2 175.8 280.0 90.0 166.2 190.0 193.3 210.0 295.0 And Wilcoxon rank sum test with continuity correction resulted in W = 622.5, p-value = 1.475e-05. So, is there any special concerns for comparing these groups aside from applying t-test or Mann-Whitney-U, regarding the fact that legal range does not start at 0? **@Peter Flom** Well, after reading your last comment, I tried it on R and saw that the two p values are the same: > tension1<-c(160,180,170,150,145,176,198,200) > > deviation1<-c(20,40,30,10,5,36,58,60) > > tension2<-tension1+12;deviation2<-deviation1+12 > > ks.test(tension1,tension2) Two-sample Kolmogorov-Smirnov test data: dev1 and dev2 D = 0.375, p-value = 0.6601 alternative hypothesis: two- sided > t.test(tension1,tension2) Welch Two Sample t-test data: tension1 and tension2 t = -1.1792, df = 14, p-value = 0.258 > t.test(deviation1,deviation2) Welch Two Sample t-test data: deviation1 and deviation2 t = -1.1792, df = 14, p-value = 0.258


title: How to simulate a binary response variable based on two non-interacting continuous variables

text: I want to simulate a binary response variable which depends on two normally distributed continuous variables, and I want to have more 1s than 0s in the response variable. I wonder how this can be done such that a logistic regression will not identify a significant interaction term. My current approach in R looks like this:               n = 1e5     x1 = rnorm(n)     x2 = rnorm(n)     y = x1+x2+rnorm(n)     y = ifelse(y > 2, 1, 0)     df=data.frame(x1=x1, x2=x2, y=y)     summary(glm(y ~ x1*x2, df, family=binomial(logit)))$coefficients      This usually results in a highly significant interaction term, even though the y is just the sum of x1 and x2. So how can I simulate a y which depends on both x1 and x2, but not on their interaction?


title: Partitioning sum of squares

text: I am working through the proof for partitioning sum of squares.   Could someone please explicitly explain how one goes from line 14 to 15 in the following PDF: http://capone.mtsu.edu/dwalsh/4380/438PARTN.pdf That is, how one moves between these lines - \begin{align} \rm{SSTO} &= \sum_{i=1}^r \sum_{j=1}^{n_i} (Y_{ij}-\overline{Y}_{..})^2 \\\ &= \sum_{i=1}^r \sum_{j=1}^{n_i} [(Y_{ij}-\overline{Y}_{i.})+(\overline{Y}_{i.}-\overline{Y}_{..})]^2 \end{align} I'm actually reading through Casella and Berger's "Statistical Inference" (pg 536), but the same jump occurs in both. Sorry, I was never the best at algebra... would someone please explain this.


title: Find the normal probability $P(|X|<1)$ using z values?

text: If I have a random variable that follows a normal distribution i.e.      $X \sim N(-3,4)$      and I want to calculate $P(|X| <1)$ how would I go about doing so, using the z values? Seeing as it's the modulus, would I only need to consider one half of the normal curve and thus calculate $P(X<1)$?


title: How to verify extremely low error rates

text: I am faced with trying to demonstrate through testing an extremely low error rate for a sensor (no more than 1 error in 1,000,000 attempts). We have limited time to conduct the experiment so we anticipate not being able to obtain more than about 4,000 attempts. I see no problem showing the sensor does not meet the requirement, as even one error in 4,000 attempts will yield a 95% confidence interval for the error rate with a lower limit greater than 0.000001. Showing that it does meet the requirement, however is the problem, as even 0 errors in 4,000 attempts still results in a lower bound greater than 0.000001. Any suggestions would be greatly appreciated.


title: Are there examples of R packages changing dramatically between version such that the results of a statistical function were significantly different?

text: Im trying to understand how people use R packages and was wondering if there are documented cases where R packages have produced different answers. Clarification: The motivation behind this question comes from an effort I'm involved in where the goal to understand the importance of provenance in the analytical methods and how it facilitates reproducible research. While R is big in the science community at present, and R packages are versioned in CRAN, without detailed information [especially version numbers], someone trying to reproduce a body of work in the future might come to a different conclusion that the original work (even with the original data). Example: Paper by John Doe says "we used R 2.3.1 and package glmulti to fit our models". 10 years from now, someone might use a new version of glmulti (no one knows what version was used in the original) which might produce a much different conclusion. My question: Are there examples of such a thing happening already? Version 2 or R package produces a much different result that version 1.


title: Interaction in Survival analysis

text: I am a physician and clearly not a statistician (I try to understand why and how to perform the right analysis, but I don't understand the formulas). I am using SPSS v.19 in analyzing my data. I try to analyze if a biomarker is associated with an increased mortality risk in a subpopulation from my study and I don't know if my analysis is correctly performed. This the result (part of it) for the entire population: ![enter image description here](http://i.stack.imgur.com/THJwk.jpg) This is for the one of the subgroups: ![enter image description here](http://i.stack.imgur.com/zT7SR.png) This is for the other subgroup: ![enter image description here](http://i.stack.imgur.com/a0Q8B.jpg) In this situation, I wanted to see if the association with the mortality risk in the second group was truly significant and I performed a Cox survival analysis with an interaction term: Biomarker value*LOT_1: ![enter image description here](http://i.stack.imgur.com/EFMwA.jpg) Is this analysis correct? And if it is, from this output I should understand that there is no statistically significant association with mortality for the studied biomarker in the second subgroup (no difference for this association in between the two subgroups)? Thank you. LE: This is the output including also the main effect of LOT_1. ![enter image description here](http://i.stack.imgur.com/mVFkP.jpg)


title: Interpretation of p-value in comparing proportions between two small groups in R

text: I wish to compare proportions between two small groups. In group 1, 2 of 18 students received honours. In group 2, 9 of 21 received honours. This is a small sample, so if I use the t-test:               x1 <- c(rep(1,2),rep(0,16))     x2 <- c(rep(1,9),rep(0,12))     t.test(x1,x2, alternative="two.sided", conf.level = 0.95)      I obtain a p-value of .0239 and CI of (-0.590, -0.0445): If I had used the z-test (a.k.a. chi-squared test with 1 degree of freedom):               prop.test(x=c(2,9), n= c(18,21), alternative="two.sided", conf.level=0.95)      I get a p-value of .0659 and a CI of (-0.626, -0.00921) I had always thought that for small samples one must use a t-test since the normal approximation is not valid for small samples, causing exaggerated significance under z-test. So why is the difference significant under the t-test, but not under the z-test? And in the t-test, why does the CI not include zero when the p-value is less than alpha?


title: Quadratic programming when the matrix is not positive definite

text: http://cran.r-project.org/web/packages/quadprog/quadprog.pdf R package `quadprog` seems to be able to solve the quadratic programming problem only when the matrix $D$ is positive definite. However, there is a case when the matrix $D$ is not positive definite. such as \begin{eqnarray} \min(x^2 + y^2 - 6xy) \\\ \text{subject to}\quad\quad x + y &\leq& 1,\\\ 3x + y &\leq& 1.5,\\\ x,y &\geq& 0\. \end{eqnarray} How can I solve this kind of problem?


title: Support Vector Machines and Recommender Algorithms

text: The recommendation problem is this: Suppose that we have a matrix whose columns are items (e.g., movies) and whose rows are users. A small part of the matrix is filled with rating values. That is values (e.g. between 1 and 5) which represent how much a particular user likes a particular item. The problem is to predict the missing ratings in the matrix. After prediction the highly predicted items will be recommended to the users. So far the best methods are based on matrix factorization techniques. SVMs are one of the best off-the-shelf prediction algorithms. I tried to think of ways of using SVMs for the rating prediction problem. However, I could not make much progress. Also, I could not explain why SVMs are not good for this problem. For example, one can build a dataset as a list of user, item, rating triples.               123, 3214, 4     214, 1282, 1     ...      Here the first column contains user ids, the second column contains item ids and the last column contains the ratings (the class label to be predicted). You can think of it as a regression problem or a multi-class classification problem. Build an SVM classifier then try to predict a new rating given a user and an item. But I don't think that this is a good thing to do. Any ideas? Can SVMs be used for this problem? If yes, how? If no, why? Thanks


title: What is the difference between multiple regression & mutivariate regression?

text: I have data on GDP growth as a dependent variable and growth in main production sectors of Pakistan such as mining, electricity, communication, manufacturing and electricity. I am supposed to run a regression on it. Is there any difference between multiple regression and multivariate regression? If so, than what is it? Which is suitable for my data?


title: In OLS on a census, should statistical significance be ignored?

text: Let's say I'm trying to tell whether my company's sales have had a downward trend over the last 10 years. Let's also say I have available to me the last 10 years of quarterly sales. If I run a bivariate regression with time as the independent variable and sales as the dependent variable, and if the regression coefficient for time is negative, do I need to examine the regression coefficient's t-statistic? I mean, aren't these data based on a census of my population of interest? If they're a census, why would I need statistical testing? Basically I'm wondering how often my knee-jerk reaction to apply statistical testing has been wrong-headed.


title: Is the value of $\alpha$ the same for all support vectors (SV) in the dual and what is the reason for it if they do or don't?

text: Consider the dual with no offset and not slack. In the dual we have that for data points that are support vectors: $$\alpha_t > 0$$ and that the constraint is satisfied with equality (since a support vector!): $$y^{(t)}\sum^n_{t'=1} \alpha_{t'}y^{(t')}x^{(t')} \cdot x^{(t)} = 1$$ However, just because $\alpha_t > 0$ its not obvious to me if all the support vector have the same value of $\alpha_t$. Do all of the support vectors have the same value of $\alpha_t$? I was kind of looking for a mathematical justification or intuition for my question. My guess is that they would have the same $\alpha_t$ value in the absence of slack. My other guess was that, $\alpha_i$ intuitively, says the important of that point, right? (maybe thats wrong) If thats true then, not all support vector have the same value of $\alpha_i$ because not all support vectors are essential support vectors. What I mean is, that some support vectors will lie on the margin boundary but removing them will not necessarily change the margin (for example, consider the case where you have 3 SVs in a line, you can remove the middle one and keep the same boundary). So do essential support vectors get larger value of $\alpha_t$? Also, I was wondering, what would happen in the case of slack and offset to the value of $\alpha_t$. I was guessing that closer points to the decision boundary, get larger value of $\alpha$? However, are points inside the margin even consider support vectors?


title: Questions about confidence intervals

text: I'm working on some homework and cannot figure out this question so I thought I'd come here for some help. > Part 1: A recent study of food portion sizes reported that over a 17-year > period, the average size of a soft drink consumed by Americans aged 2 years > and older increased from 13.1 ounces (oz) to 19.9 oz. The authors state that > the difference is statistically significant with P < 0.01. Explain what > additional information you would need to compute a confidence interval for > the increase, and outline the procedure that you would use for the > computations. (Select all that apply.) > > a.) t and degrees of freedom could be used to find the confidence interval. > In this case we could compute SED and use degrees of freedom to find t*. > > b.) Standard deviations and degrees of freedom could be used to find the > confidence interval. In this case we could now find the P-value, which could > be used to find SED. > > c.) Sample sizes and a more accurate P-value could be used to find the > confidence interval. In this case we could determine standard deviations and > the confidence interval in the usual way. > > d.) Degrees of freedom and a more accurate P-value could be used to find the > confidence interval. In this case we could determine t, then calculate SED > and t*. > > e.) Sample sizes and standard deviations could be used to find the > confidence interval. In this case we could find the interval in the usual > way. * * * > Part 2: Do you think that a confidence interval would provide useful > additional information? Explain why or why not. > > a.) Yes, the confidence interval could give us useful information about the > average size of soft drinks. > > b.) Yes, the confidence interval could give us useful information about the > variability between sample participants in the study. > > c.) Yes, the confidence interval could give us useful information about the > magnitude of the difference. > > d.) No, the confidence interval could give us no more useful information > because it cannot tell us the sample size in the study. > > e.) No, the confidence interval could give us no more useful information > because the P-value already tells us that the interval does not contain 0. For part 1, I think a is among the correct answers, but I can't tell if any others would work as well. And for part 2, I've eliminated b as a possible solution.


title: Mode estimation in high dimensions

text: Suppose we have a sample $\boldsymbol{x}_i$ for $i$ in $1,\dots, n$, from a $d$-dimensional unimodal density $f(\boldsymbol{x})$. I would like to estimate the mode of $f(\boldsymbol{x})$. The mean-shift algorithm discussed in this related post can be used to estimate the mode, but given that it is based of a non-parametric estimator of the density, I think it will requires a huge sample size $n$ if $d > 5$. Are there other options? Thanks


title: Different formulations for SVM with slack variables (primal)

text: I have seen two different ways to formulate the SVM optimization but I was not sure what the difference was between them or if there was any difference. First formulation: $$min \frac{1}{2}||\theta||^2 + C \sum^n_{i=1}\xi_i $$ $$\text{s.t. }\ \ y^{(i)}(\theta \cdot x^{(i)} +\theta_0) \geq 1 - \xi_i$$ Second formulation: $$min \frac{\lambda}{2}||\theta||^2 + \frac{1}{n} \sum^n_{i=1}\xi_i $$ $$\text{s.t. }\ \ y^{(i)}(\theta \cdot x^{(i)} +\theta_0) \geq 1 - \xi_i$$ Are they the same formulation? The thing that is confusing me is that I think what they mean that: $$min \frac{1}{2}||\theta||^2 + \frac{1}{\lambda} \sum^n_{i=1}\xi_i $$ where: $$C = \frac{1}{\lambda}$$ Is the same as: $$min \frac{\lambda}{2}||\theta||^2 + \sum^n_{i=1}\xi_i $$ Am I right that there should not be and 1/n? Or does it make a difference?


title: Dependent variables in regression

text: I have two variables in a regression problem, where I predict the end interest rate for a loan application.   * x1: risk band (A+,A,B,C)   * x2: initial rate (6%, 7%, 8%, 9%) Based on these two variables, I'm trying to predict what the actual end rate, y, will be (i.e 6.3%). Risk band is dependent on initial rate and vice versa. i.e A+ -> 6%. I would need to convert the risk band to 0s and 1s for each degrees of freedom in the category. I have read the linear regression requires that the input features be independent. In my case, they are not. Would it make sense to include both in the model? So a possible input feature vector might look like this: [0,0,1,9] -> C risk band, 9% initial rate [0,0,0,6] -> A+ risk band, 6% initial rate I can z-score the initial rate to make the ranges smaller. * * * Given the answer below, If I choose to represent initial rate as a categorical feature as well, then I would have a problem because they would be scalar multiples of the risk band. For example [0,0,1,0,0,1] -> C risk band, 9% initial rate [0,0,0,0,0,0] -> A risk band, 6% initial rate However, If I use an ordinal feature for the intial rate, then there shouldn't be a problem, because they are not scalar multiples of risk band. For example, [0,0,1,9] -> C risk band, 9% initial rate [0,0,0,6] -> A+ risk band, 6% initial rate Is that correct?


title: How to compare 2 classifiers after testing trained model in Weka?

text: I'm new to Weka (and machine learning) so this question would be a bit silly. So I have 2 models built using J48 and RandomForest (both run with 10-fold cross-validation mode) on a 40,000-tuple training set. I also have 6 different smaller test sets, each has around 8000 tuples. Now I want to compare how these 2 classification algorithms perform after testing their models (I have the predictions reported). The problem is, I don't have a good understanding of the reported predictions so I don't know how to compare the prediction of each to see which one is better in terms of accuracy.


title: Similarity between objects based on tags (binary features)

text: I have five millions of objects each of them having one or more tags. How do I compute statistically sound similarity score between each pair of the objects taking into account that:   1. There are 100 millions of tags most of them used only once.    2. Some tags are used very often, say one tag may be used on 1% of whole dataset.   3. Some objects are heavily tagged with hundreds of thousands tags on them while others may be tagged only a couple of times.  Second question: how do I cluster objects taking into account 1-3? My guess is that k-means and other popular clustering techniques won't do much here. I've tried k-means already with simple distance defined as number of similar tags on the objects and clusters are so vague to the point being almost meaningless.


title: Is the t-test really robust to shape of distribution?

text: Can I use the t-test to compare means between groups (N=300) if data is not normally distributed?


title: Probability Interpretations

text: How comprehensive is the following book - What interpretations are missing? Interpretations of Probability, Andrei Khrennikov, 2009, de Gruyter, ISBN 978-3-11-020748-4 http://www.degruyter.com/cont/fb/ma/detailEn.cfm?isbn=9783110207484&sel=pi Contents:http://www.degruyter.com/files/pdf/9783110207484Contents.pdf


title: How can I assess GEE/logistic model fit when covariates have some missing data?

text: I have fit two generalized estimating equation (GEE) models to my data: 1) Model 1: Outcome is longitudinal Yes/No variable (A) (year 1,2,3,4,5) with longitudinal continuous predictor (B) for years 1,2,3,4,5. 2) Model 2: Outcome is the same longitudinal Yes/No variable (A), but now with my predictor fixed at its year 1 value i.e. forced to be time invariant (B). Due to missing measurements in my longitudinal predictor at a few time points for different cases, the number of data points in model 2 is higher than in model 1. I would like to know about what comparisons I can validly make between the odds ratios, p-values and fit of the two models e.g.:   * If the OR for predictor B is bigger in model 1, can I validly say that the association between A and B is stronger in model1?   * How can I assess which is the better model for my data. am I correct in thinking that QIC/AIC pseudo R squareds should not be compared across models if the number of observations is not the same? Any help would be greatly appreciated.


title: Comparing Cox Proportional Hazards Models (variable selection)

text: I am using a cox proportional hazards model to run a survival analysis in r on a number of non-nested, distinct covariates such as Age, Blood Type, Cancer, etc:                A, B, C, D, E          When I run the model on the omnibus null hypothesis:               surv ~ A + B + C + D          The effects of all of the covariates are insignificant because the number of subjects that have measurements for every covariate is relatively small. However, when I isolate single or other combinations of covariates in different cox models:               surv ~ A         surv ~ A + C     surv ~ B + D      I'm showing significant effects because the sample set is larger (i.e. the number of observations discarded by the model shrinks). What I'm having difficulty understanding is how to do the following:   * Comparing the different cox models for the best fit, i.e. is `surv ~ A + B + D` a better model than `surv ~ A + C` ? Should I be comparing the likelihood, wald or logrank scores?   * Is it possible to run every possible combination of covariates to determine the best model? I have about 15 covariates.   * More broadly, is this tactic the best approach to optimizing for both significant covariates and overall model "cost"? I will be attaching a cost to each distinct cox model i.e. using covariates `A + B + C` in the model costs \$100 while using covariates `A + B` costs \$75 and using only covariate `A` costs \$10. I'd like to look at the cost for each combination of covariates vs. the accuracy for each cox model.  Thanks very much for your help!