title: How to calculate standard errors of a non-linear model prediction?

text: I'm trying to understand how to show the prediction error of a model fit in R using the non-linear least squares function `nls`. Although there is an argument `se.fit=TRUE in the`predict.nls` function, the help explains "At present se.fit and interval are ignored." I'm able to calculate the confidence interval of the fitted parameters, but am unsure if this information can be used to calculate prediction std. error. Below is an example: ### Example data:               df <- structure(list(L = c(13L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,      15L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,      16L, 16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L,      17L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L,      19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,      20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L,      21L, 21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 22L, 22L, 23L, 23L,      23L, 23L, 23L, 23L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,      24L, 24L, 24L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 26L,      26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 27L,      27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L,      28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L,      28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L,      28L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L,      29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L,      30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L,      30L, 30L, 30L, 30L, 30L, 30L, 31L, 31L, 31L, 31L, 31L, 31L, 31L,      31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 32L,      32L, 32L, 32L, 32L, 32L, 33L, 33L, 34L, 34L, 34L, 34L, 34L, 34L,      34L, 35L, 35L, 35L, 35L, 37L, 38L), mat = c(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, 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, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L,      0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L,      1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L)), .Names = c("L", "mat"), class = "data.frame", row.names = c("1",      "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",      "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",      "25", "26", "114", "115", "27", "28", "29", "30", "31", "32",      "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43",      "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54",      "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65",      "116", "117", "66", "67", "68", "69", "70", "71", "72", "118",      "119", "73", "74", "75", "76", "120", "121", "77", "78", "79",      "122", "123", "124", "80", "81", "82", "83", "84", "85", "86",      "125", "126", "127", "128", "129", "87", "88", "89", "90", "130",      "131", "132", "133", "134", "91", "92", "93", "94", "95", "96",      "135", "136", "137", "138", "139", "140", "141", "97", "98",      "99", "100", "142", "143", "144", "145", "146", "147", "148",      "149", "150", "151", "101", "102", "103", "104", "105", "106",      "152", "153", "154", "155", "156", "157", "158", "159", "160",      "161", "162", "163", "164", "165", "166", "167", "168", "169",      "170", "171", "172", "107", "108", "109", "110", "111", "173",      "174", "175", "176", "177", "178", "179", "180", "181", "182",      "183", "184", "185", "186", "187", "188", "189", "190", "191",      "192", "112", "193", "194", "195", "196", "197", "198", "199",      "200", "201", "202", "203", "204", "205", "206", "207", "208",      "209", "210", "113", "211", "212", "213", "214", "215", "216",      "217", "218", "219", "220", "221", "222", "223", "224", "225",      "226", "227", "228", "229", "230", "231", "232", "233", "234",      "235", "236", "237", "238", "239", "240", "241", "242", "243",      "244", "245", "246", "247", "248", "249"))          df2 <- aggregate(df, by=list(as.factor(df$L)), FUN="mean")[,-1]     names(df2) <- c("L", "pmat")      ### Fit the NLS model               fmla <- "mat ~ 1 / (1 + exp(-a*(L-L50)))"     fit <- nls(       fmla,       data=df,       start=c(a=0.3, L50=25)     )     summary(fit)     fit.ci <- confint(fit)     fit.ci          newdat <- data.frame(L=seq(0,50,.1))     pred <- predict(fit, newdat)      Below, I compare this fit to a GLM model using a binomial distribution. `predict.glm` allows me to easily calculate the std. error of the prediction. You can see that the predictions between NLS and GLM are quite similar on this data. ### Comparison to binomial GLM               fit2 <- glm(mat ~ L, data = df, family=binomial)     pred2 <- predict(fit2, newdat, type="response", se.fit=TRUE)     plot(df, xlim=c(0,50))     points(pmat ~ L, df2, pch=20, col=8)     lines(newdat$L, pred)         lines(newdat$L, pred2$fit, col=2)         polygon(x=c(newdat$L, rev(newdat$L)), y=c(pred2$fit+pred2$se.fit, rev(pred2$fit-pred2$se.fit)), col=rgb(1,0,0,0.2), border=NA)     legend("right", legend=c("NLS", "binom. GLM"), col=1:2, lty=1, bty="n")     legend("left", legend=c("raw", "binned"), col=c(1,8), pch=c(1,20), bty="n")      ![enter image description here](http://i.stack.imgur.com/UZXec.png)


title: Pybrain trainer runtime error

text: I am trying to use pybrain for the first time on 20x20 handwritten digit MNIST data. I keep getting a runtime error when I setup the trainer and I can;t figure out what the problem is? code:               import sys     import os     from pprint import pprint          import numpy          X = numpy.loadtxt('/Users/enpasant/Octave/mlclass-ex4-005/mlclass-ex4/Xs.csv', delimiter=',')     y = numpy.loadtxt('/Users/enpasant/Octave/mlclass-ex4-005/mlclass-ex4/Ys.csv', delimiter=',')     y = numpy.reshape(y,(y.shape[0],1))     print X.shape,y.shape          y = numpy.reshape(y,(y.shape[0],1))          from pybrain.tools.shortcuts import buildNetwork     from pybrain.datasets import SupervisedDataSet     from pybrain.datasets import ClassificationDataSet     DS =  ClassificationDataSet(X.shape[1] , nb_classes=10)     DS.setField('input',X)     DS.setField('target',y)     TrainDS, TestDS = DS.splitWithProportion(0.75)     print len(TrainDS)               from pybrain.structure import TanhLayer     from pybrain.supervised.trainers import BackpropTrainer     net = buildNetwork(400, 25, 10)     trainer = BackpropTrainer(net, dataset=TrainDS)     trainer.trainUntilConvergence()      output:               (5000, 400) (5000, 1)     3750     Error in sys.excepthook:     Traceback (most recent call last):     File "/Applications/TextMate.app/Contents/SharedSupport/Bundles/Python.tmbundle/Support/sitecustomize.py", line 44, in tm_excepthook     io = fdopen(error_fd, 'wb', 0)     OSError: [Errno 9] Bad file descriptor          Original exception was:     Traceback (most recent call last):     File "untitled", line 36, in <module>     trainer = BackpropTrainer(net, dataset=TrainDS)     File "/Library/Python/2.7/site-packages/PyBrain-0.3-py2.7.egg/pybrain/supervised/trainers/backprop.py", line 35, in __init__     self.setData(dataset)     File "/Library/Python/2.7/site-packages/PyBrain-0.3-py2.7.egg/pybrain/supervised/trainers/trainer.py", line 23, in setData     assert dataset.outdim == self.module.outdim     AssertionError


title: Forming training set for Multinomial Naive Bayes

text: Is it true that Multinomial Naive Bayes requires equally by count training data for each class to get best performance? For example, we forming classifier for three classes - Japan, China, Korea. For Japan available 500 training data, China - 300, Korea - 100. So to train Multinomial Naive Bayes we need to get only 100 data for each class inspite of we have 500 for Japan, 300 for China. I know that we can train Naive Bayes with all available data, but in this case we got worse performance, is it true?


title: Have population, use inferential statistics? Also, non-normal dependent variable, what to do?

text: Background: We are looking at parental leave in Iceland. We are particularly interested in whether the economic crisis and the resulting changes in parental leave legislation affected the time taken for parental leave. We have reason to believe that the effect of the crisis / new laws will be different for mothers and fathers (who have equal rights to a leave), will depend on income and education, and that there might be an interaction between factors (e.g., that the length of leave for fathers would be independent of income before the new laws, but would start to depend on income after the laws were passed).   1. We actually have not just a sample, but the entire population (about 50000 kids). Do we even need inferential statistics at all? Can we just describe the results either numerically or graphically, because whatever difference there is, that is the actual difference in the population?   2. If we do need to do inferential statistics, then we have some potential problems:  First, the dependent variable (length of leave) is not even close to normally distributed. Instead it is multimodal, e.g., people are likely to take 0 days, 30 days, 60 days, 90 days etc. but not, e.g., 3 days or 34 days. I cannot transform this distribution to look anything like it is normally distributed. I initially considered using some kind of non-parametric test that looks at differences in medians, but the problem is that the medians might actually always be close to same (e.g., 90 days) but the distribution is changes nonetheless. I then considered binarizing the dependent variable (e.g., takes less than standard leave vs. takes standard leave or more). This would allow me to use logistic regression and the weirdness of the distribution would be gone. I am fine with this. However, I am interested not just in main effects (e.g. main effect of time and main effect of income) but also in interactions (e.g., the interaction between time and income). I am not sure how to deal with interactions in logistic regression, especially since I might have to treat the factors as categorical (e.g., I am not expecting the length of leave to, say, linearly increase or decrease with time -- I am expecting a curvilinear relationship between length of leave and time). What to do? I mainly use SPSS for analysis, in case that is relevant.


title: Trend of Cramer Rao Lower bound with signal-to-noise ratio?

text: Intuitively, CRLB should decrease with increase of signal-to-noise ratio (SNR). A plot of CRLB (Y axis) and SNR on X axis should be somewhat a decreasing exponential curve. Thus, when the signal has more of information and less noise, CRLB decreases. I have not found any reference for such a plot, but shall be extremely grateful for an explanation as to what this trend shows and why is that happening. My other question is, about a confusion : IF CRLB exists then the estimator is an efficient, unbiased estimator? Or is it just an unbiased estimator? What does efficient mean? Thank you


title: Bayesian approach to infer whether a patient has a particular disease

text: Suppose for each of $p$ symptoms we have a measurement of the degree of severity from $1$ to $5$. We have a dataset $\{(x_i,y_i)\}_{i=1}^n$, where $x_i=(x_{ij},j=1,\cdots,p)$ are the measured severities of the symptoms of patient $i$, and $y_i=0$ or $1$ indicates whether or note patient $i$ has the disease. Let $\pi$ be the prior of class $Y_i=1$. Then the joint probability over the class variable $Y_i$ and observation $X_i$ is $$ P(X_i,Y_i)=\pi^{Y_i}(1-\pi)^{1-Y_i}\prod_{k=0}^1\prod_{j=1}^p\prod_{v=1}^5\theta_{kjv}^{I(Y_i=k)I(X_{ij}=v)} $$ where $\theta_{kjv}$ is the probability that $X_{ij}=v$ given $Y_i=k$. Now I am asked to take a Bayesian approach to derive the posterior distribution of parameters $\theta$ and $\pi$, given the dataset $D=\{(x_i,y_i)\}_{i=1}^n$, under the assumption that the prior of $\pi$ is $\text{Beta}(1,1)$ (i.e. uniform) and the prior of $\{\theta_{kjv}\}_{v=1}^5$ is $\text{Dirichlet(1,1,1,1,1)}$ (i.e. uniform on the simplex $\sum_v\theta_{kjv}=1$). Also, assume the parameters are independent. Could anyone show me how to do this? My try is as follow: $$ \begin{align*} f(\pi|D)&=\frac{P(D|\pi)f(\pi)}{P(D)}\\ &\propto P(D|\pi)\\ &=\int_{\theta\in\text{simplex}}P(D,\theta|\pi)d\theta\\ &=\int_{\theta\in\text{simplex}}P(D|\pi,\theta)f(\theta)d\theta\\ &=\mathbb{E}_\theta[P(D|\pi,\theta)] \end{align*} $$ It is the most I can do, how can I move on to have a closed form for $f(\pi|D)$? It is requires to write it in $n_k$ and $m_{kjv}$ where $n_k$ is the number of data in class $k$, and $m_{kjv}$ is the number of patients in class $k$ with severity level $v$ for symptom $j$.


title: Regression - volatility vs return

text: I am attempting to estimate a linear model as: $$ y = a +bX +e $$ I have a series of annual returns and I would like to estimate the effect of volatility on losses. My Null hypothesis is that volatility in period $t$ should have no linear relationship with losses in period $t+1$. So in actual fact, my model is: $$ \min[0,y]t = a+bX_{t-1} +e $$ where $x$ is the standard deviation of $y$ in period $t-1$. My questions are as follows:   1. What is the best way to test if my model is correctly specified?   2. Given that $X$ is rolled over in each period, should I be worried on the impact of serial correlations?    3. In the case of serial correlations, can this be better reduced by using higher lags of $X$?


title: How to define the maximum k of the kNN classifier?

text: I am trying to use kNN classifier to perform some supervised learning. In order to find the best number of 'k' of kNN, I used cross validation. For example, the following codes load some Matlab standard data and run the cross validation to plot various k values with respect to the cross validation error               load ionosphere;     [N,D] = size(X)     resp = unique(Y)          rng(8000,'twister') % for reproducibility     K = round(logspace(0,log10(N),10)); % number of neighbors     cvloss = zeros(numel(K),1);     for k=1:numel(K)         knn = ClassificationKNN.fit(X,Y,...             'NumNeighbors',K(k),'CrossVal','On');         cvloss(k) = kfoldLoss(knn);     end     figure; % Plot the accuracy versus k     plot(K,cvloss);     xlabel('Number of nearest neighbors');     ylabel('10 fold classification error');     title('k-NN classification');      The result looks like ![Result](http://i.stack.imgur.com/JtlQH.png) The best k in this case is k=2 (it is not an exhaustive search). From the figure, we can see that the cross validation error goes up dramatically after k>50. It gets to a large error and become stable after k>100. My question is what is the maximum k we should test in this kind of cross validation framework? For example, there are two classes in the 'ionosphere' data. One class labeled as 'g' and one labeled as 'b'. There are 351 instances in total. For 'g' there are 225 cases and for 'b' there are 126 cases. In the codes above, it chooses the largest k=351 to be tested. But should we only test from 1 to 126 or up to 225? Is there a relation between the test cases and the maximum number of k? Thanks. A.


title: In general, would you always prefer feasible GLS to OLS?

text: I know that GLS estimators only have exact distributions asymptotically, so the efficiency gains in finite samples are not all that clear. But otherwise, I'm struggling on how to attack this discussion.


title: Calculate conditional probabilities on a logistic regression

text: I am analysing data of an digital advertising tool. I have access to all the exposures for each of the 6 channels (x1, x2, .., x6) of an online campaign. The outcome is the conversion (y=1 or 0) on the website. My goal is to calculate an attribution model and credit each of the channels with the number of conversions that he deserves. I started by creating a logistic regression between the numbers of exposures by channel (x1,..,x6) to explain the conversion.   Coef S.E. Wald Z Pr(>|Z|)   Intercept -6.7910 0.0938 -72.36 <0.0001   x1 0.0756 0.0193 3.92 <0.0001   x2 0.0698 0.0091 7.64 <0.0001   x3 0.0020 0.0002 11.05 <0.0001   x4 0.0172 0.0057 3.03 0.0024   x5 0.0304 0.0045 6.82 <0.0001   x6 -0.0132 0.0042 -3.17 0.0015 In order to give credit to each of the channels, I’d like to use the Shapley Value, applied to the different conditional probabilities. So I am calculating now these conditional probabilities but I am not sure how that my formula is the right one. For example, if I want to know the probability of users to convert if exposed to channel 1, then the formula will be: P(Y|X1) = mean(exp(-6.7910 + 0.0756*x1)/(1+exp(-6.7910 + 0.0756*x1))   For all the x1 values of my data. To calculate P(Y|X1 and X2) then I will to the same but with X1 and X2 values. I am not entirely sure of the method, could someone confirm to me that it is the right way to calculate conditional probabilities on a logistic regression ? Thanks,


title: How to approach a real machine learning problem?

text: As a ML beginner, I don't how to begin with this problem, or more generally, are there any typical steps to take when approaching a real problem? If I have some domain relevant knowledge, then I could make use of it to analyze the data and do the modelling. But if I don't, then all I have is just the data, a matrix of input and a vector of response (if supervised). My question is if I don't have domain knowledge, what to do before I randomly try out some learning algorithms and tune the parameters? Am I supposed to analyze the data per se to find out any clue and how?


title: Choosing between MANCOVA and RDA to factor out treatment effects

text: I am working on a data set of intertidal organisms and am looking for a multivariate technique that will help me determine why some groups have different species assemblages than others. I was thinking multivariate analysis of covariance (MANCOVA) might be a good method because it takes care of the experimental design (which is 12 experimental blocks, 3 treatments and all 36 sites being sampled once a year for 15 years), the raw data, and covariates (which will hopefully include data about water temperatures, salinity, etc). However, it seems like perhaps people use MANCOVA more to factor OUT the covariates to focus on the experimental design, whereas I really am more interested in doing the opposite (factoring out the effect of the treatments). This is because the interesting stuff we're noticing happened in spite of the treatments - we think because of unexpected climate stuff. Is this an appropriate method to use for my purposes or not? I already have a principal component analysis (PCA) that gives a visual representation of some of the dominant assemblages, and am thinking if MANCOVA won't work my next best bet might be doing a redundancy analysis (RDA) instead.


title: What is the difference between the studentized and non-studentized Breusch-Pagan test?

text: I have a multiple regression model with 4 independent variable. I want to test heteroscedasticity with the Breusch-Pagan test. If I do a studentized BP test the p-value is 0,173, with a non-studentized test, it is 0,89.


title: Power of Holm's multiple comparison testing compared to others

text: Following one-way ANOVA, there are many possible follow-up multiple comparison tests. Holm's test (or better, the Holm-Sidak) test has lots of power, but because it works in a stepwise manner, it cannot compute confidence intervals. Its advantage over the tests than can compute confidence intervals (Tukey, Dunnett) is that is has more power. But is it fair to say that the Holm method _always_ has more power than the methods of Tukey and Dunnet? Or does it depend...?


title: What is your favorite data visualization blog?

text: What is the best blog on data visualization? I'm making this question a community wiki since it is highly subjective. Please limit each answer to one link.


title: What's the purpose of window function in spectral analysis?

text: I'd like to see the answer with qualitative view on the problem, not just definition. Examples and analogous from other areas of applied math also would be good. I understand, my question is silly, but I can't find good and intuitive introduction textbook on signal processing — if someone would suggest one, I will be happy. Thanks.


title: Multiprocess models

text: I would like to use multi-process models. in particular, I would like to apply this model: Kulu, Hill. "Migration and fertility: Competing hypotheses re- examined." European Journal of Population/Revue européenne de Démographie. In the fertility equations I don't understand how the estimation procedure differs between the two models (conceptually I get the difference but I am rather confused when it comes to the actual methodological application). In particular, to let the errors vary across the three migration equation is the estimation procedure extended in a discrete-time setting? Additionally, i don´t understand whether the author distinguishes between a first birth after first migration and first birth after second migration or whether he is just considering first birth after whatever migration. Can someone help me to better understand? Textbooks about the topic are welcome thank you!


Let $X_{i}$, $i=0,1,\cdots$ be a sequence of random variables generated by $X_{i+1}=\rho X_{i}+W_{i}$, where $\rho$ is constant and $W_{i}$ are i.i.d random variables. Suppose $X_{0}$ is independent of $W_{i}$. Is $X_{i}$ a martingale, submartingale, or supermartingale?


title: How to set the labels of One-Class SVM on LibSVM?

text: I try to perform one-class SVM Schölkopf on LibSVM. But all the predicted labels are "-1 " which is very strange. Below is my code. I use MATLAB and libSVM.               % ex5  with iris dataset          load iris.mat ;          % choose the first class of data as target data      % and other two classes of data as negative data     % ex5  with iris dataset          load iris.mat ;          % choose the first class of data as target data      % and other two classes of data as negative data          target_data = data(label==1,:)          % select part of target data as training data           train_data =  target_data (1:30,:) ;          train_label =  ones(30,1) ;          % use the remaining  data as test data           test_label = [  ones(20,1); -ones(100,1) ]          test_data = data (31:150 ,:) ;          % train           model = svmtrain ( train_label,train_data,'-s 2   ' ) ;          % test           [prdicted_label,accuracy ] = svmpredict( test_label,test_data,model)             And the predicted_labels all are -1. What is the reason?


title: Continuous-by-continuous interaction using simple slopes when I have nuisance variables in my model

text: I've created a multiple regression predicting 1 DV from 4 IVs (all mean- centered) + 1 interaction term between two of the IVs. The interaction term is significant. I'd like to get a nifty graph comparing the effect of one of the variables included in the interaction on the DV at particular values (say, $\pm1 SD$) of the other variable (i.e., the moderator). I understand how to do this when I only have 3 IVs (predictor, moderator, and interaction), but how can I do this when I have 2 other IVs not involved in the interaction? Let me know if additional info or exact numbers will help.


title: Is this a correct use of "optimal"?

text: During fine-tuning some classifier parameters, is it correct to talk about "optimal" configuration? For example, if perfect classification is impossible, but the error converges to some relatively good small value, is it correct to call it an "optimal" solution, or is it "suboptimal" because there's no single optimal one? If there are several configurations that perform best and equally good, would it be correct to call either of them "optimal"?


title: R vs. SPSS Kolmogorov-Smirnov test

text: I was recently given the following data.               7 1 6 6 2 8 3 7 4 6 3 3 1 8 4 5 9 3 2 4 1 4 4 3 6 3 4 2 3 2 7 2 7 4 2 4 4 4 3 2 2 4 5 8 6 4 3 1 5 5      I then ran a KS test as follows in R:               ks.test(mydata,"pnorm",mean(mydata),sd(mydata))      and got the following:               One-sample Kolmogorov-Smirnov test          data:  mydata      D = 0.1832, p-value = 0.06985     alternative hypothesis: two-sided       Trouble is when I do the same test in SPSS, as shown here: https://statistics.laerd.com/spss-tutorials/testing-for-normality-using-spss- statistics.php I get a p-value of 0.000 but exactly the same test statistic. Does anyone have an idea why this is? I'm not an expert in SPSS at all (don't even have it installed) but the fact that the D-statistic matches makes me think I ran it correctly and am somehow misreading it. **Addendum**   RioRaider's hint was the right way to go. After manually applying Lilliefors Significance Correction as described here: http://home.ubalt.edu/ntsbarsh/stat-data/SPSSSAS.htm Formula:               Exp[ -7.01256D2(n + 2.78019) + 2.99587D(n+2.78019)1/2 - 0.122119 + 0.97498/n1/2 + 1.67997/n ]       I get exactly the p-values as described. Is there a way I can flag his comment as a +1?


title: coding survey data for cosine similarity and euclidean distance?

text: I want to know how to code survey data such that a similarity function can be applied on it. Say I want to use cosine similarity. All the search results and QA I've found while in my search deal only with the similarity between _documents_ , with vectors consisting of word frequencies or tf/idf. What about survey data? What is the sensible/common/useful way of coding survey data such that similarity can be compared? (Is it even sensible to use functions like cosine similarity for this?) My data is record data, purely categorical, neither binary nor numerical. Should I code it into numerical data? My data looks like this (3 sample records):               Do you like Technology?  | Current GPA       | Institute name     Y                        | Band 1 (3.75-4)   | UUIC     N                        | Band 3 (3.0-3.5)  | ADU     N                        | Band 2 (3.5-3.75) | UUIC      etc. These are just 3 questions, my survey had a lot more questions, but I hope you get the idea. Is it sensible to code the data into numerical vectors, where for eg. I represent yes/no values as binary variables, and assign numbered categories to other values? In which case the above 3 records would become:               (1, 1, 1)     (0, 3, 2)     (0, 2, 1)      Where UUIC = 1, ADU = 2, and the GPA bands are represented simply by 1, 2, 3, 4, etc.. And then apply cosine similarity or euclidean distance? Would this make sense? I've been searching for similar examples for a while now but everything that comes up seems to be about _document similarity_. There doesn't seem to be much beginner's help on how to deal with _survey data_.


title: Question about homoscedasticity test

text: one of the reviewer of a paper of mine suggested to perform a homoscedasticity test between the results of two experiments, testing the same thing in two conditions. The experiments consisted in ratings along a 7 points likert scale. One of the experiment results were distributed on a large range of values, while the other showed a tendency towards the center of the likert scale range (i.e. 3.5). The reviewer argued that the latter behavour could be due to the fact that participants answered more randomly, and he suggested that to verify this possibility, an homoscedasticity test should be performed comparing the two conditions. Now, I would like to understand this comment. In relation to my case, what does it mean that the variances of the two conditions are significantly different? What instead if they are equal? Secondly and more importantly, which test for homoscedasticity do you suggest I perform? I use R. Can you suggest the function more suitable for my case? I saw that there are many.


title: Geometric intuition for why an outer product of two vectors makes a correlation matrix?

text: I understand that the outer product of two vectors, say representing two detrended time series, can represent a cross-correlation (well covariance) matrix. I also know that the inverse of a correlation matrix represents the partial correlations between two variables. Geometrically, I know that the partial correlation between two variables is the angle formed by the projection of their residuals when regressed against all other variables onto the surface perpendicular to all other variables. I'm wondering how these two relate. That is, I know the interpretation of the inverse of matrix (partial correlation) but not the matrix or its construction.


title: Can we still use OLS on truncated a $Y$ if its conditional distribution is normal?

text: I was recently reading about Heckman selection models, and got sidetracked by how little I knew about truncated data. I was reading these slides, and on page 78 Baum mentions that if part of a sample has a truncated $Y$ variable, we can’t even use regression analysis on that $Y$ to make inferences about the sub-sample (whose $Y$ variable is truncated). I know that having a truncated $Y$ means the distribution will be skewed, but if the conditional distribution of $Y$ is normal, can we still use basic OLS? Also, does a tobit model work on the assumption that the conditional distribution is truncated?


title: Calculating the space/distance between scale points (IRT)?

text: I am working with a seven item Likert Scale (1, 2......6, 7) for a classroom management. How can I calculate the distance between these seven scale points? This is for instrument reliability purposes, I want to see if people find, for example, a larger difference between choice 4 and 5 and a shorter distance/difference between choice 5 and 6. There is an IRT method to figure this out, but I am coming up blank on it. Any/all help is much appreciated.


title: How to determine the correlation between 2 time series while controlling for a 3rd?

text: I would like to determine the relationship between two variables after controlling for a third. Specifically, I want to know if the prices of mercury and gold over time are correlated with each other more than or less than they are correlated with a generic metals price index. Is there some special signal in the relationship between these commodities or are they just moving in line with the overall index? I have a data set with 32 annual values (1980-2011) for mercury price, gold price, and metals index. I am using R. Below I have posted time series plots of each of the three variables. ![enter image description here](http://i.stack.imgur.com/9gF4O.png)


title: 3 independent var with 2 levels

text: 1. There are 3 independent variables (carb, protein, fat) in my experiment, each with two levels (low and High). I am trying to figure out if 3-Way ANOVA would be the best method. Example: Carb (low and high); Protein (low and high) Fat (low and high)   2. How do I figure out the sample size of the participants in this in-between subjects experiment?    3. What is the essence of conducting a pre-test and a post-test in such an experiment? Is it really required? Is it the same as post-hoc?? Thanks


title: How to understand moments for a random variable?

text: Wikipedia says that the name of concept comes from physics, but I cannot find any similarity between these two concepts.


title: How to complete matrix multiplication with an index that runs across all nodes?

text: Not sure if this would be the right place to ask...but I am attempting to calculate a topological overlap measure (TOM) matrix using R. I'm running into a problem with a likely simple solution. Basically, I am not sure how to interpret the matrix notation for the TOM matrix. An example would be: $l_{ij}=\sum_{u}a_{iu}a_{uj}$ where the index u runs across all nodes. Would this mean I am multiplying entire rows by entire columns, and then finding the sum of these multiplications?


title: Simple question about the asymptotics of estimators

text: Consider any arbitrary estimator called $\hat{M}$ (e.g., regression coefficient estimator or specific type of correlation estimator, etc) that satisfies the following asymptotic property: $$\boxed{\sqrt{N}(\hat{M}-M) \overset{d}{\to}\mathcal{N}(0,\sigma^2)}\,\,\,\,\,\,\,\,\,\,\,\,(1)$$ which implies that our $\hat{M}$ is consistent. We also have a consistent estimator $\hat{\sigma}$, which gives rise to the asymptotic property: $$\displaystyle \ \ \boxed{\frac{\sqrt{N}(\hat{M}-M)}{\hat{\sigma}} \overset{d}{\to}\mathcal{N}(0,1)}\,\,\,\,\,\,\,\,\,\,\,\,(2)$$ I'm wondering if I can use the $z$- or $t$-test **just like normal** for _any_ such $\hat{M}$ that satisfies the above? Let $Q$ be defined as the test statistic: $$\displaystyle \ \ \boxed{Q_\hat{M} = \frac{\hat{M}-M_{H_0}}{\sqrt{\frac{1}{N}\hat{\sigma}^2}}}\,\,\,\,\,\,\,\,\,\,\,\,(3)$$ My goal is to do the following hypothesis test: $H_0: M = 0$ $H_a: M \not= 0$ yet the only information I have access to is $(1)$ and $(2)$, whence my question. * * * $$\underline{\text{Update}}$$ The current answers suggest that I can't always robustly $z$- or $t$-test for any such $\hat{M}$. I am reading the relevant sections of _All of Statistics_ (Wasserman), as well as _Statistical Inference_ (Casella & Berger). Both state that, if: $$\displaystyle \ \ \frac{\sqrt{N}(\hat{M}-M)}{\hat{\sigma}} \overset{d}{\to} \mathcal{N}(0,1)$$   * then " _an approximate test can be based on the wald statistic $Q$ and would reject $H_0$ if.f. $Q < -z_{\alpha/2}$ or $Q > z_{\alpha/2}$_" (in Casella & Berger, page 492, "10.3.2 Other Large-Sample Tests")   * or, in (Wasserman, page 158, Theorem 10.13) " _Let $Q = (\hat{M}-M_{H_0})/\hat{se}$ denote the observed value of the Wald statistic $Q$ $\big($where $\hat{se}$ is obviously equal to my $\sqrt{\frac{1}{N}\hat{\sigma}^2}$$\big)$. The p-value is given by:_ $$p = 2\Phi(-|Q|)$$ This contradicts the existing advice since they do not state any other necessary assumptions to be able to do this legitimately (to the best of my ability to comprehend). Either;   * I have failed to understand existing answers.   * I have failed to express my original question clearly.   * I have failed to read these chapters properly.   * They are excluding thoroughness for pedagogical purposes. I would appreciate some assistance on which option is correct. Thanks. $\big($Please go easy I am new to stats :)$\big)$. * * * Another dimension is that my intended application is $n = 3000$, so perhaps the finite sample problems are less relevant?


title: How to use libSVM for one-class SVM problems?

text: I plan to use libSVM for a one-class svm problem, but I'm not sure about the meaning of `nu` in `svm_parameter`. Does it mean the probability that a test point lies _outside_ of a set S (estimated from the training data) equals `nu`?


title: How do you detect if a given dataset has multivariate normal distribution?

text: I'm looking at Fisher's LDA on various datasets on UCI ML repository and trying to see where LDA might perform badly. One reason I can think of is if the data distribution is not a multi-variate normal distribution. This is from the fact I read in a book where you apply LDA on multivariate normal distribution. Is that thought process correct? When might LDA give bad results?


title: What's the optimal way of measuring the error of quantitative estimates, across question sets?

text: My research relates to situations in which individuals make a quantitative estimate and then make a second quantitative estimate, perhaps after some manipulation. Minimally I am interested in comparing: A) The error of a first estimate   B) The error of a second estimate   C) The error of the average of the first and second estimates. I want to compare the efficacy of different manipulations, and also the extent to which changing question sets affects the efficacy of a particular manipulation. The question sets vary widely; in some question sets you can get responses that are huge outliers (e.g. answers to a question like "How many kilograms of fuel were on the Space Shuttle orbiter at take-off?") while with other question sets big outliers can't really exist (e.g. answers to a question like "What percentage of the world's airports are in the United States?"). I've seen the Mean Squared Error frequently used in this context (e.g. this paper), and also the Root Mean Square Error. Others have favoured a method based on the median (e.g. this paper), and I've also seen the Mean Absolute Deviation used (e.g. this paper). I'm keen to know which method should be preferred, or if no method is preferable over the others how I could make some principled argument that that is the case.


title: Chi Square Goodness of Fit vs. Igloo

text: I am very new to Chi square goodness of fit tests but have done a fair bit of research. Basically, I have the following 12 data points: 7392, 7656, 7241, 6164, 4984, 4664, 15262, 17053, 21814, 5094, 4581, 10500. And I would like to test the assumption that a Lognormal($\mu$=9.11,$\sigma$=0.51) is a good fit to this data. I have some output from the program "Igloo", which states that the Chi square statistic is 3.10. I am struggling to derive this number. I know there are a number of ways to define e.g. the size and location of the "bins", the expected frequencies, etc. The other information from Igloo is that the number of bins chosen is 3. That's all! I am struggling to justify using 3 bins, and struggling to find the size that the 3 bins should be! If I can't replicate the 3.1, then a method which gets close to this answer would suffice. Any help on this would be greatly appreciated. I appreciate it is a very vague question.


title: Effects of replacing missing data with mean or median?

text: When preparing data for use with Covariance-based Structural Equation Modeling (CB-SEM), what are the different effects of replacing missing data with mean or median? When is one of them more appropriate than the other? Is any of them generally recommended?


title: When the response variable is a proportion (r/n) in proc logistic, how to express the model in terms of $\pi, \alpha, \beta$?

text: A (SAS documentation file) (page 1906) on "The LOGISTIC Procedure" gives the following procedure-               proc logistic;     model r/n=x1 x2;     run;      Here, n represents the number of trials and r represents the number of events. Question 1: The response variable - r/n will produce proportions, therefore, it's not either 0 or 1. In that case, what kind of logistics model it is? Ordinal? Question 2: If there are 4 groups formed according to very high, high, medium and low levels of exposure in predicting a disease outcome, and if I want to use r/n as response variable using the proc logistic (as mentioned above) with a single exposure variable, what will be the model expression in terms of pie, alpha and beta? In other words, I have a 4x2 table with 4 levels of exposure, and 2 outcomes- disease and non-disease. Hope at least one expert will come forward to answer this simple question. Many thanks in advance!


title: Multiple dependent binomial tests

text: This may well be a textbook question - I am not a statistician and just need to solve a practical problem - I will be grateful for any pointers to literature. There is a $M\times N$ matrix $P$ defined by (matrix) multiplication of vectors $X$ and $Y$ (with lengths $M$ and $N$); that is, $P=XY'$. Here $X$, $Y$, and $P$ are hidden. There is another $M\times N$ matrix $T$, with values i.i.d. (normally or uniformly, not very important). $T$ is observed. $M\times N$ binomial experiments are conducted, using $P$ as probability of success, and $T$ for number of Bernoulli trials. The resulting matrix of numbers of success $S$ is observed. The goal is estimating $P$ (alternatively, estimating $X$ and $Y$ up to multiplicative constants). A sloppy method of just dividing marginal sums of $S$ by marginal sums of $T$ to arrive at estimates of $X$ and $Y$ (up to a constant) works well for tight distributions of $T$ (or if rank of $T$ is otherwise close to $1$). I thought that solving a system of $M\times N-1$ "almost-linear" equations (with elements of $T$ serving as coefficients, elements of $X$ and $Y$ as variables, and marginal sums of $S$ as constants - I can elaborate if their form is not obvious) will give better estimates for $X$ and $Y$. Turned out, the estimates are often worse than the "sloppy" ones. How should I exploit the dependency structure among elements of $P$? Is hierarchical Bayesian the best approach here? $P$ are usually small, so despite $T$ being of the order of thousands, typical $S$ elements are rather small ($0$ is not uncommon, though not that usual). I understand this affects accuracy of estimates dramatically, but I still believe I am not extracting the fullest accuracy possible. Oh, and both $M$ and $N$ are in the range $10-100$.


title: Permutation paired test with signed statistic

text: I have one paired data from 7 patients: $$(X_i ,Y_i)=\{(5,1),(-2,-4),(4,2),(-6,-3),(-3,-3),(-5,2),(3,1)\},\quad (i=1,\ldots, 7)$$ where $X$ means placebo $Y$ means the test drug. My classmate asked me can we have permutation test for the sign statistics to identify $X < Y$. I think the coin function can help me shuffle all the observation into $2^7$ permutations, and then we are not sure if we can just implement Wilcoxon rank sum statistics for each permutation to see the $p$-value? So the question is how to shuffle all these data and implement the test?


title: Comparing Gini coefficients: Variance estimation etc. needed?

text: In a project on software measurement, we plan to use aggregating statistics (e.g., Gini) to describe the concentration of certain observed program attributes (size) among program units (e.g., modules, functions). This is a purely exploratory, descriptive effort (e.g., no hypothesis testing involved). We are aware of the limitations of e.g. Gini (sample bias, symmetry, different distribution/same Gini etc.) and the various sample/population estimation approaches to inequality measures such as Gini (jackknifing estimates & bootstrap estimates such as here). For a given program, we compute several Gini coefficients, each representing a different "perspective" on the program structure:                                  P1        P2               ...        P_n     Program            0.5578098 0.7981829        ...        0.7981829        1. As a purely descriptive statistic, is it acceptable to compare the observed Gini coefficients directly, to state Gini(P1) <|>|... Gini(P2) etc.? Should we take corrective actions (e.g., small-sample bias correction?!)   2. Is it necessary to report estimate of variance and/or standard error in this case? In our view, we compute the statistic on the actual population data rather than some sample. Besides, we do not run hypothesis tests on some confidence interval(s). If so, how can this be achieved best (var, se, quantile error bounds, ...)? Thank you for your critical comments/suggestions!


title: SPSS/Stata Repeated Measures - Contrast question

text: I have data from 8 subjects collected at 5 time points. I would like to see if there is a difference between time points (1 + 2) and (3 + 4 + 5) using either SPSS or Stata. When in SPSS I use 'Analyze > GLM > Repeated Measures' I am able to specify the 5 time points, but I can not create contrasts. I can only select the predefined contrasts. I guess I could just average time point (1 + 2) and (3 + 4 + 5) and compare the means in a paired-sample t-test. **Q1)** Are there any disadvantages with averaging the time points and using a paired t-test (e.g. loss of power?) **Q2)** How could I model these contrasts in SPSS or Stata?


title: Find the values that let me split x in equally numbered groups

text: I have some data: x = c(0.0636931012704989, 0.0421033586843369, 0.0214721845269228, 0.038624779052474, 0.033235347720021, 0.0307709640052774, 0.0606316442337883, 0.0289958997354256, 0.0573689255248752, 0.044448618126441, 0.0164384406368751, 0.0357802024114138, 0.0133831663454831, 0.0329126088838874, 0.051468481535772, 0.0302927520748979, 0.169023097551101, 0.0311182034472608, 0.117954563949638, 0.0593735124639249, 0.0110273691186214, 0.0286868155187407, 0.0522496344459376, 0.0298854332250777, 0.0626024185545931, 0.106914814257457, 0.0396392209572767, 0.0111988702332617, 0.0276848766950444, 0.0142151606082234) I would like to find the values, in R preferably, in order to split x equally numbered groups. ![enter image description here](http://i.stack.imgur.com/EcHhh.png)


title: If I have T total records, how big should my sample size be for a valid analysis?

text: I have a bunch of records, T in total. I want to know how many of these I can get away with analyzing in order to extrapolate the analysis to the entire population T. I know this is a basic question and largely depends on how much error I can accept, but can anyone tell me the math?


title: Feature-based malware detection using Weka

text: My aim is to detect malware by testing it with the extracted features. I have extracted the binary features of 150 malware files and 150 benign files by using the hex editor tool for unpacking the files and Netbeans for generating a dataset in an Excel sheet. Now I have Features of malware and benign datasets in separate sheets. Now I want to test a new file with the existing features to say whether it is a malware file or a benign file. For this, I'm working with the Weka tool. Can I achieve my goal with this tool? If not, which tool is better for my project? I just want to compare the new file automatically with those two malware and benign datasets when given as input and it should say if it is a malware file or a benign file.


title: How to interpret categorical variables in an OLS when only one category is statistically significant?

text: I am running a simple OLS.   * **Dependent Variable:** Population Change In A Congressional District After An Election   * **Independent Variable:** Who won the election: Democrat, Republican, or, Independent (Categorical) I've created dummies variables for the different categories and included two of them in two regression models:   * **PopulationChange = DemocratWin + RepublicanWin + e**   * **PopulationChange = DemocratWin + IndependentWin + e** However, when I run the model, only DemocratWin is statistically significant. Substantively, how do I interpret that?


title: I am use caret package to do model comparion and met error message

text: Thanks to Zach's code. I run R code as below and get error that cold not find "trainControl and train function"               #Use the same CV-folds for each model     require(caret)     myControl <- trainControl(method='cv', number=10, index=createFolds(y, k=10))          #Fit models     model_rpart <- train(X, y, method='rpart', trControl=myControl)     model_ctree <- train(X, y, method='ctree', trControl=myControl)     model_rf <- train(X, y, method='rf', trControl=myControl)      Would you give me some suggestion since I am a R beginner. Thank you.


title: Generate correlated AR process for given correlation between demand series

text: How can I generate two correlated $AR(1)$ data series with given correlation between $d_{1,t}$ and $d_{2,t}$, $r_{12}$, where $\rho_{12}$ is correlation between the two error series $$d_{1,t}=\mu+\phi_1d_{1,(t−1)}+e_{1}(t)\quad\quad d_{2,t}=\mu+\phi_2d_{2,(t-1)}+e_{2}(t)$$ $e_1$ and $e_2$ are not the same and $r_{12}$ is desired correlation between $d_{1,t}$ and $d_{2,t}$. The relations between $r_{12}$ and $\rho_{12}$ is: $$r_{12}=\rho_{12}*\sqrt{(1-\phi_1^2)*(1-\phi_2^2)}/ (1-\phi_1\phi_2)$$ one example for: $\phi_1=.3$, $\phi_2=-.8$ AND $r_{12}=-.8$; the problem is, for this example $\rho_{12}$ is $-1.7$, now how can I generate the error terms where the correlation between the two error series $\rho_{12}$ is not between $-1$ and $+1$. please help me if i'm wrong somewhere. THANKS YOU


title: Extracting fixed effects table output from lme object

text: I am trying to fit a model which is similar to the               fm2orth.lm<-lme(distance~age,data=OrthoFem, random=~1|Subject)      and did               summary(fm2Orth.lm)      my output had these parts   1. Aic/BIC    2. Random Effects    3. Fixed Effects table   4. Correlation. My question is I am interested in only the third output above and tried               fm2Orth.lm$coefficients       and had coefficients of the individual subjects, which is not what I want, can someone show me how I could get only the fixed effects table with the p-values.


title: Level shifts in an ARIMA model: fewer are better?

text: I've created an ARIMA (0 1 1)(0 1 2)12 model for our electricity costs. I have 44 months of monthly data, and have used X-12-ARIMA to do the model fitting. (I get similar results in R, it's just easier to use X-12-ARIMA and it gives a LOT of diagnostic output.) X-12-ARIMA allows you to enter a level shift easily, for example LS2008.Jul enters a level shift at July, 2008. Unfortunately, between the utility and the state commission, rates get jerked around more often than I'd like. I know there were rate changes in July of 2008, September of 2009, May of 2010, and July of 2011. Is there any reason to not enter all four of these level shifts? When I do the first and last (LS2008.Jul and LS2011.Jul, both close to the data boundaries), the results look good. When I insert the middle two (LS2009.Sep was a provisional jump up and LS2010.May a reversal), the ACF output shows issues. If it makes any difference, at the same time as LS2010.May, there were refunds in May-August that I've entered as USER (exogenous) variables. Should this be? My analogous thinking is that a LS releases tension in a spline fit, as it were, and too many just supply too much slack. Even if I know that they are there. Is there any reason to not include an LS when I know there should be one?


title: What is the "root mse" in stata?

text: I have a question that has been confusing me ever since I took econometrics last year. What does the "root MSE" mean in stata output when you regress a OLS model? I know that it translates into "root mean squared error", but which variable's mean squared error is it after all, and how is it calculated? Can anybody provide a precise definition and formula, and explain why it is helpful to have that value?


title: Linear Discriminant (decision region) singly connected and convex?

text: I am reading about linear discriminants, and have encountered a phrase that I have no idea about. The phrase says that a decision region constructed in a certain way, is "singly connected and convex", I want to know in really simple terms what this means. I have an idea of what convex means in that a convex function has a definite optimum. But I don't know what singly connected means and I'm not sure if convex means the same thing here. This question might be very simple but I am very new to the field of statistical learning and am struggling with the concepts related to optimization theory so any help would be much appreciated. Thanks in advance!


title: Given a mean and standard deviation, what can you determine of a non-normal distribution?

text: I have the following question: > Given a mean of 11 and a standard deviation of 2, with a non-normal > distribution, can you determine the % of numbers that are between 8 and 12? My guess is no, because the distribution could be weighted on the left and right... but I'm getting thrown off by the sd of 2. Not sure what that implies for this question.


title: How to perform text mining on an online news article?

text: I am new to data mining and currently working on an online news article from TOI. My aim is to get some useful information out of this text which is not clear when you read the article and the most important things mentioned in it. I am done with document processing on the article "TOI manifesto" which I have made it in a simple text format. I have applied tokenization, filtering and ngrams and got the result as shown in the images. ![doc processing](http://i.stack.imgur.com/Ta8bJ.png) ![result](http://i.stack.imgur.com/VXI9I.png) I don't know how to proceed further, which algorithm / model would be best to apply for mining this article and what should be my next step?


title: Can one force randomness in a sample?

text: I recently spoke to a large qualitative analysis company who were working on in-depth analysis of the potential customers of train company. I asked them how they chose the people to include in their sample, they said, that they picked them at random from list of people in the region, but they also said that had predefined criteria, that the sample needed have equal distribution of both genders as in the region, and similarly for ethnicity, income, education and so on. I asked them whether they believed that the list they picked from was biased, but this was not the reason, it was rather to make the sample more random. Does this make sense, or do they risk interfering with the randomness, and hence make it less random? Their sample was around 100 people, but I am interested in an answer also for very few people.


title: Continuous-states hidden Markov chain

text: How to deal with HMC that has continuous states? Any papers, links, materials that explain the solution?


title: Math behind multivariate testing for website optimization

text: I am looking for theoretical resources (books, tutorials, etc.) to learn about making sound statistical inferences given (plenty of) multivariate website conversion data. I'm after the math involved, and cannot find any good non-marketing stuff on the web. The sort of questions I want to answer: how much impact does a single variable (e.g. color of text) have? what is the correlation between variables? what type of distribution is used for modelling (Gaussian, Binomial, etc.)? When using statistics to analyze results - what should be considered as a random variable - the web-page element that gets different variations or the binary conversion-or-no-conversion outcome of an impression? There's plenty of information about different website optimization testing methods and their benefits\pitfalls, plenty of information about multivariate statistics in general, do you guys know of resources that discuss **technical statistics in this specific context of website optimization** ? Thanks for any info!


title: Basic Instrumental Variables Dummy Variable question

text: I have been reading Mostly Harmless Econometrics and have been won over by the virtues of splitting my instrument into groups. In particular, I have a continuous variable as an instrument representing percent and would like to split it into 3 groups (0-33%, 34-66%, 67-100%) and use dummies to indicate each group, i.e. if the continuous variable is between 0-33% $\text{dummy}_1=1$ and the others 0, if between 34-66% $\text{dummy}_2=1$ etc. If I run 2SLS would I need to drop the first group (0-33%) in the first stage? And if so, what would be the interpretation of the constant in the first stage and reduced form? Also would I be able to over identification test? What would be the number of over identification restrictions? Yes split into 3 groups and use a dummy for each one as a separate instrument. In particular a higher percentage represents a higher randomised treatment intensity, so this would be 3 randomised treatment intensities. I was hoping I could include all 3 in the first stage without facing linear dependence.


title: Calculating entropy of a binary matrix

text: Given a matrix whose entries consist of only 1's and 0's, I would like to come up with a measure of how "ordered" the matrix is, in some sense. This exact question was posed here: Measuring entropy/ information/ patterns of a 2d binary matrix in which the top-rated answer posted by whuber provided what I'm looking for, except that I didn't understand one key detail. Referring to his answer, he writes: 'Let's measure this randomness with their base-2 entropy. For $a_1$, the sequence of these entropies is $(1.92, 1.5, 1.58, 1.0)$. Let's call this the "profile" of $a_1$.' Unfortunately, I've tried but failed to come up with the same entropy values, so I'd appreciate it if I could get some help on how to arrive at those numbers.


title: How should I combine the data with different 0/1 ratio to build the logistic regression model?

text: I have a question: Suppose I have a logistic regression model, with data from 2008 to 2011. Because of some reasons like I cannot collect all data, only a sample(or subset) of the data was used to build the model, for example, bads and goods ratio is 1:4. And now I have all 2012 and 2013 data. The bads goods for 2012-2013 data is like 1:10. The question is: shall I combine these two data together directly to renew the model? Or shall I sample from the 2012-2013 data to keep the 1:4 ratio and then combine with old data to renew the model? As I think, this is like an oversampling question. By theory, it should not effect too much. Does anyone have real experience and how did you deal with this? Thank you!


title: Given $n$ different univariate non-normal sample sets calculate for a new sample, $x$, which it most likely belongs to

text: Say you have $n$ different, non-normal, potentially overlapping data sets of samples. Maybe their densities look something like: ![](http://i.stack.imgur.com/OjZ6q.png) and you are given a new sample $x$, how would you decide to which of these sample sets it would most likely belong? I would assume it is possible to do a (kernel) density estimation of each of the sample sets, interpolate them to functions $p_n$, and then calculate $p_n(x)$ for each, selecting the one with the highest probability. However both the decision rule and the method seem, well, off. Can anybody help me in getting a better intuition of this problem? **EDIT** Say the data comes from a continuous scoring system, like a depression scale, and I have annotated data for many different subjects. So I can get the density plots for "severe", "mildly" and "non-depressed" subjects. Now I have a new sample and wish to know (based on the score) which the person most likely belongs to.


title: Where to start with nonlinear regression?

text: What are some good introductory materials for learning non-linear regression modeling? Bonus points for anything freely available. Thanks!


title: How to bootstrap non-negative matrix factorization results?

text: I have RNA-seq data from 9 samples and around 15,000 genes. I know that these 9 samples consist of varying proportions of two cell types, each with their own expression profile. I am using non-negative matrix factorization with a rank of 2 to determine the cell-type ratios in each sample. In other words, I use NMF to decompose my 9 x 15,000 matrix into a 2 x 15,000 matrix and a 2 x 9 matrix. I want to obtain confidence intervals for my estimates of the 2 x 9 matrix. From what I understand about bootstrapping, I could do this by resampling from my 9 observations a bunch of times and keep track of each estimate for each observation. My question is, since NMF doesn't have any notion of observations vs. features, could I also estimate the variance by resampling genes (features) rather than observations?


title: range normalized with unclear highest value

text: I tried to normalize many set, which each set consists of 12 values. I want to normalize it into range between 1 and -1. However, the highest value is unclear because the value can be trillion, or even more. I code it using java. can somebody help me what to do? I tried to normalize it using standard deviation but the result is not satisfying. thanks in advance


title: Sample size when comparing several AUCs

text: How many individuals do I need to compare 3 AUC which are computed from the same set of patients? Is there a R program or some other program available? There are many programs for computing the sample size when comparing 2 AUC, but I have found none for more than two.


title: Am I specifying my lmer model correctly?

text: I've scoured Google and this site and I am still confused about the lmer function in the lme4 library. I have some data collected from different psychiatric wards, which have a multilevel structure. To simplify, I'll pick two level 2 and two level 1 variables, although I actually have a few more. Level 2- WardSize [this is the number of people on the ward] & WAS [this is a measure of how "nice" the ward is] The grouping variable that tells R who's in which ward is called "Ward" Level one- Gender [this is gender, obviously] & BSITotal [this is a measure of symptom severity] Outcome is Selfreject, which again is what it sounds like. I have this formula: help=lmer(formula=Selfreject~WardSize+WAS+Gender+BSITotal+(1|Ward)) I'm hoping this means "each individual has a score related to their own Gender and symptom severity, and also a ward-level effect relating to the size of the ward and how "nice" it is" Is this correct? The thing that's confusing me is that I can't see how R can tell which are level 1 and which level 2 variables, except for the ward level intercept given at the end. If anyone could explain the notation so an idiot like me can understand that would be even better. Many thanks!


title: Should I be concerned if the cells of values obtained from bootstrapping are correlated?

text: I have data from 40 people, each measured once in 2 conditions. For reasons tangential to my question I'm not comfortable using parametric analyses to compare the two conditions, so I'm employing bootstrapping. However, when I run a bootstrap analysis (resampling people with replacement, computing condition means, then repeating 1000 times), the condition means I obtain appear to be correlated across iterations. I noticed this because when I plot the 95% confidence interval in each condition the intervals capture each others' means, but when I compute a difference score between the conditions means within each iteration, the 95% confidence interval of this difference score excludes zero. Going back to the results of the bootstrapping, I find that the correlation of the conditions across the bootstrap iterations is .7! Should this affect my confidence in the interval generated for the difference score? Alternatively, does this tell me something (possibly about individual differences) about my data otherwise?


title: Analyzing a 2x3 repeated measures design using a logit mixed model

text: An experiment I conducted recently used a 2 (between participants) x 3 (within participants) design. That is, participants were randomly allocated to one of two conditions, and then completed three similar tasks each (in a counterbalanced order). In each of these tasks, participants made binary choices (2AFC) in a number of trials, each of which had a normatively correct answer.   In every trial, participants were presented a distractor, which was assumed to bias responses towards one of the alternatives. The tasks differed only in presence and magnitude of this distractor (i.e. no distractor vs. distractor of small and large magnitude). I would like to examine the error rates (deviations from the normatively correct answer) across these conditions. I hypothesize that the error rate will increase when a distractor is present, but will not increase further when the magnitude of the distractor is increased. Also, I expect this increase to differ between the between-subjects conditions. The latter interaction is the central interest. From discussions here and from the literature (Dixon, 2008; Jaeger, 2008), I gather that logit mixed-models are the appropriate analysis method, and that, in R, the lme4 package is the tool of choice.   While I could compute some basic analyses (e.g. random intercept model, random effects ANCOVA) with lme4, I am stuck as to how to apply the models to the design in question -- I have the feeling that I am very much thinking in terms of HLMs, and have not yet quite understood the entirety of mixed effects models. Therefore, I would be very grateful for your help. I have two basic questions:   1. In a first analysis, I would like to consider the error rates in only those trials in which participants were biased towards the _wrong_ answer. The first model would therefore look only at trials in which the bias would point away from the correct answer. If my observations were independent, i would probably use a model like this: `correct ~ condition + distractor + condition:distractor` ... but obviously, they aren't: Observations are grouped within a task (with a constant presence of a distractor) and within participants. My question, then, is this: How do I add the random effects to reflect this?   2. (If I haven't lost you already :-) ) Would it be possible to include all trials (those where the bias would be into the direction of the wrong and of the correct answer), and include this difference (i.e. direction of the bias) as a trial-level predictor? In my imagination of HLMs, such a predictor (at the level of the trial) would depend on on the magnitude of the distractor present (at the level of the block), which again would depend on the condition of the participant (plus, possibly, a unique factor for each participant).   The interactions would then emerge ›automatically‹ as cross-level interactions. How would such a model be specified in the ›flat‹ lme4 syntax? (Would such a model make sense at all?) Ok, I hope all of this makes sense -- I will gladly elaborate otherwise. Again, I would be most grateful for any ideas and comments regarding this analysis, and would like to thank you for taking the time and trouble to respond. **References** Dixon, P. (2008). Models of accuracy in repeated-measures designs. _Journal of Memory and Language, 59_ (4), 447-456. doi: 10.1016/j.jml.2007.11.004 Jaeger T. F. (2008). Categorical data analysis: Away from ANOVAs (transformation or not) and towards logit mixed models. _Journal of Memory and Language, 59_ (4), 434-446. doi: 10.1016/j.jml.2007.11.007


title: What should I know about designing a good Hybrid/Hamiltonian Monte Carlo algorithm?

text: I am designing a Hybrid Monte Carlo sampling algorithm for PyMC, and I am trying to make it as fuss free and general as possible, so I am looking for good advice on designing an HMC algorithm. I have read Radford's survey chapter and Beskos et. al.'s recent paper on optimal (step size) tuning of HMC and I gathered the following tips:   * Momentum variables should be distributed with covariance $C^{-1}$, where $C$ is generally something like the covariance matrix of the distribution (for simple distributions), but could conceivably be different (for funny shaped distributions). By default I am using the hessian at the mode.   * Trajectories should be calculated with the leapfrog method (other integrators don't seem to be worth it)   * Optimal acceptance rate is .651 for really large problems and otherwise higher.   * Step size should be scaled like $L\times d^{(1/4)}$, where $L$ is a free variable and $d$ is the number of dimensions.   * Step size should be smaller when there are light tails or otherwise regions with odd stability characteristics. Step size randomization can help with this. Are there other ideas that I should adopt or at least consider? Other papers I should read? For example, are there adaptive step size algorithms that are worth it? Is there good advice on trajectory length? Are there in fact better integrators? Someone please make this a community-wiki.


title: Dealing with Categorical variables in Multiple Regression

text: I have a data having 2 continuous and 4 categorical variables. Each categorical variable has 3 levels. I want to know how to include the variables in the model. I am using SPSS Variables:   1. Sales - Dependent Variable   2. Retail Price   3. Location : TX, CA, PY   4. Display Location in Store : DL1, DL2, DL3   5. Price Point : 2.5, 4.5, 1.0 I have created 9 dummy variables. 3 for Location, 3 for Display Location in Store and 3 for Price Point. I am confused on how to use these independent variables. How to add the independent variables into the models. Should I add all the nine dummy variables or two dummy variables of each categorical variables along with Retail Price. I searched on net but , did not get any solution. I would appreciate any kind of help. Thanks in advance.


title: Levenberg-Marquardt adjustment

text: I understand that Levenberg-Marquardt adjustment is about adding $\lambda \times\; I$ to second derivative of the log likelihood function, where $I$ is the identity matrix $$-[l''(\theta | Y) - \lambda I] \quad \text{such that it becomes positive definite}$$ **My question is:** How should we choose $\lambda$ ? if $\lambda = 10$ is just enough to make the overall expression to become positive definite and I deliberately choose $\lambda = 100$ (which still makes the expression positive definite), what impact would it have on the convergence of parameters ? So, how is the lambda chosen (wisely)?


title: Minimum number of levels for a random effects factor?

text: I'm using a mixed model in `R` (`lme4`) to analyze some repeated measures data. I have a response variable (fiber content of feces) and 3 fixed effects (body mass, etc.). My study only has 6 participants, with 16 repeated measures for each one (though two only have 12 repeats). The subjects are lizards that were given different combinations of food in different 'treatments'. My question is: can I use subject ID as a random effect? I know this is the usual course of action in longitudinal mixed effects models, to take account of the randomly sampled nature of the subjects and the fact that observations within subjects will be more closely correlated than those between subjects. But, treating subject ID as a random effect involves estimating a mean and variance for this variable. Since I have only 6 subjects (6 levels of this factor), is this enough to get an accurate characterization of the mean and variance? Does the fact that I have quite a few repeated measurements for each subject help in this regard (I don't see how it matters)? Finally, If I can't use subject ID as a random effect, will including it as a fixed effect allow me to control for the fact that I have repeated measures? Edit 1: I'd just like to clarify that when I say "can I" use subject ID as a random effect, I mean "is it a good idea to". I know I can fit the model with a factor with just 2 levels, but surely this would be in-defensible? I'm asking at what point does it become sensible to think about treating subjects as random effects? It seems like the literature advises that 5-6 levels is a lower bound. It seems to me that the estimates of the mean and variance of the random effect would not be very precise until there were 15+ factor levels.


title: Dynamic Panel/GMM in R with group:time fixed effects?

text: Is there a solution coded in R to estimate models of the form $$ y_{igt} = \alpha_i + P_{gt} + \beta_1y_{igt-1}+ \beta_2y_{igt-2} + X_{igt}'\gamma + \epsilon_{igt} $$ ?? `plm` offers the `pgmm` package, which implements the Arellano-Bond estimator, but it doesn't seem able to handle FWL-based demeaning of factors other than the cross-sectional unit, or the simple addition of time dummies. `lfe` on the other hand doesn't seem to be able to handle dynamic panel GMM estimators. I've got N= 2000, n = 65K, G = 39, and t = 25, so including `time:group` effects as simple factor variables is not an option, particularly given that I'll need to fit multiple models to figure out what lag structure I need to wipe out the autocorrelation. If nothing is coded yet, can anyone recommend any clever workarounds?


title: Maximizing expectations vs Mode maximization

text: In many statistical problems, I see the following formulation for maximizing rewards: Assuming that my total reward $R$ is the sum of individual rewards $R$: $\underset{\textbf{x}}{\operatorname{argmax}}E[R_T] = \underset{\textbf{x}}{\operatorname{argmax}} \left(E[R_1] + E[R_2] + E[R_3]\right)$ where $R_i = r_i \cdot x_i$ is the reward on every bucket $i$, assuming a reward rate of $\alpha_i$ and a number of actions per bucket $x_i$ I assume that maximizing expectations is not the same as maximizing the mode of the distribution of $R_t$. Is that correct? For example, in graphical models, it is common to run `max-sum` algorithms to find the values that correspond to the **MAP** of the distribution (i.e. mode maximization). If so, **when does it make sense to maximize the mode (i.e the probability of $R_t$) vs the expected value (i.e. $E[R_t]$), and viceversa?**


title: Random Effects Tobit Model

text: I'm studying the effect of various criminal case and court district characteristics on sentence lengths. I was planning on running xttobit in Stata because I have individual defendants/cases within federal court districts and a censored sentence length outcome. I have both an upper and lower limit. A prison sentence of "0" months actually means a non-imprisonment sentence (e.g. probation, fine, etc). In addition, I capped my sentence length measure at "life" (or, 470 months - a common cutoff in this line of research). So now my distribution is bounded at both 0 and 470. Should I run a tobit model [with or without -vce(cluster district)-] before using xttobit to conduct diagnostic and specification tests? Also, I have a binary variable that indicates whether or not someone received imprisonment (0=no prison, 1=prison). Can I run this model with logistic and then run the sentence length model, or should I just use the marginal effects for the sentence length model bc this tells you the odds of imprisonment (given the "0" values)?


title: Multiparameter log-likelihood function in R and power analysis (by using simulations)

text: I would like to estimate power of the following problem. I am interested in comparing two groups that both follow e.g. Weibull distribution. So, group A has two parameters (`a1, b1`) and two parameters has group B (`a2, b2`). By simulating random variables from distribution of interest (assuming that e.g., `a1=1.5*a2`, and `b1=0.5*b2`), apply log-likelihood ratio test to test if a1=a2 and b1=b2, and estimate power of the test. The question would be what is the log-likelihood for the full model? That is, for reduced model (when a1=a2,b1=b2) log-likelihood is (in R):               LL.reduced <- function(par, data){sum(log(dweibull(data, shape=par[1], scale=par[2])))}       what is the likelihood for the full model, when `a1!=a2`, `b1!=b2`, i.e. when 4 parameters have to be estimated? Is it possible to estimate it buy building two log-likelihoods for separate groups and add it together (i.e., `LL.full <\- LL.group1 + LL.group2`?               To summarize:      # n (sample size) = 500     # sim (number of simulations) = 1000     # alpha  = .05     # Parameters of Weibull distributions:         # group 1: a1=1,     b1=20        # group 2: a2=a1*1.5, b2=b1*0.5     # LL.reduced <- function(par, data{sum(log(dweibull(data, shape=par[1], scale=par[2])))}     # LL.full    <- ?????


title: Question about trend and predictions in R

text: I have weekly data which I am considering as reference. So Mondays trend for example.. orders will kind of serve as valid data to make predictions for orders throughout the day for every Monday. I looked at the `lm()` and `predict()` functions in R, but it seems that they just predict based my previous data model. They don't take into account the level of orders on the current Monday. So I can't blindly use the prediction obtained from the historical model. What I want is that the algorithm refers to previous data for trend but also considers the current levels and then makes the prediction. Is there any way to do that? Or maybe am I looking at the wrong models?


title: Add classification filter in Weka

text: There is a filter in Weka called "Add Classification". One can choose a classification algorithm and then can get the classification labels provided by the model generated by the algorithm, the distribution of the classes, and errors it is making. Then, one can add these information (classification labels, class distributions, and errors) into the existing feature/attribute set. For instance, a dataset has 30 attributes (including class attribute) and 2 classes. Using this filter one can generate classification labels, 2 distributions, and classification errors. Adding them to the original dataset gives 34 attributes. My question is: Does this method of using one classifier's labels, distribution and errors as additional features for another classifier have any particular name? I would like to know more about the efficiency of this method. Can anyone help me please?


title: Does this way of calculating the probability that an event affected price movement make sense?

text: The (generalized version) of the problem I'm trying to solve: Say you're observing the surface of a puddle while it's raining. There are various sized ripples occurring wherever a raindrop hits the puddle, which happens with a regular frequency. Now imagine you throw a pebble at the puddle; how would you determine, by watching the waves, the probability of whether the pebble landed inside the puddle (creating ripples stronger than a raindrop would make) or outside (creating no ripples)? In other words, can you calculate the probability of whether an event (a pebble was thrown in) occurred? In more concrete terms, I have market time series data with each entry containing a high and low for that day. So, first I calculate the average true range (an exponential moving average) and divide by the median (high-low) price to get an average daily price movement in percent. I then calculate the average standard deviation of the movement using a simple moving average. On the day of the event, I measure the movement and calculate the probability using...Z-scores? This is the part I've gotten to and I'm not sure if it's in any way reasonable to look at it like this. Is this a valid way to solve the problem? Are there better ways? I'm looking for a big picture view, not necessarily something which would be _extremely_ accurate but at least which would be a valid interpretation.


title: Correlation between two variables over time

text: I have 5 individuals. I measure variable X and Y of each one of them over time (on 5 different time points). I want to find out if there are any correlations between X and Y over time and show the same trend for all 5 individuals. The data looks like this (just an example, not my real data):               time 1:     X              Y     10            10        #individual 1     20            30        #individual 2     30            50        #individual 3     40            55        #individual 4     10            20        #individual 5          time 2:     X             Y              2             6         #individual 1          3             8         #individual 2     7             11        #individual 3     10            14        #individual 4     11            14        #individual 5     ... e.t.c.      **Which statistical analysis should I use to test if X and Y (and that X and Y correlates) show the same trends for all 5 individuals over time?**


title: How to get information of which features are selected by Random Forest while constructing Number of trees in Forest

text: I am working on Text classification problem (Sentiment Analysis), and i am having a dataset of movie reviews which contains 1000 positive and 1000 negative reviews, and for classification of those reviews in positive or negative, i am using Random Forest classifier(Tool : WEKA). In WEKA i am using the default values for the Number of trees which is 10, and Features will be based on SQRT(total features). Now, As WEKA gives the output of classification process something like this "Random forest of _N_ trees, each constructed while considering _M_ random features." Now my question, is there any possible way which can help me to get some information about those features, i mean i want to know which specific features were used for construction of trees. it will be really great help if anyone can help me with related source. Thank You.


title: Significance Test for Comparison of Variables

text: I am not sure how to ask this question without giving an example. I am trying to measure the "cleanliness" of office buildings. I have two variables that try to measure this. Variable one is a person who inspects a random sample of the offices in a building and rates the cleanliness of each office. These scores are averaged together to get an overall score for the building. Variable two is the number of janitorial requests to clean up a messy area in a building. Someone actually calls the maintenance office and says "there is a mess in this office can someone with a broom/mop come clean it up." This is the exact count at the building level, not a sample. The idea is that both variables are trying to measure the same concept, "cleanliness." As such, I expect some "congruity" between the two variables. The way I am trying to measure the congruity is by using Spearman's Correlation (ranking the buildings by cleanliness for each variable and applying the correlation formula to the ranks). Here are my questions: 1) Can Spearman's correlation coefficient be used as an indicator of whether or not the inspectors are accurately scoring the offices (assuming janitorial requests are an appropriate measure of cleanliness). 2) Does the fact that the building's score from the inspector is a function of a sample of offices from the building while the number of janitorial requests is a population measure affect the calculation of correlation coefficient? 3) If Spearman's correlation is an OK measure to use, what range of values would be considered "Strong", "Medium", "Weak"? Is this up to the end user of the data? 4) If I am only concerned with 10 buildings and I have data for all 10 buildings, do I need to use significance tests? For example, is the coefficient different from zero? My understanding of significance tests is that they are needed when data is being sampled from a population but not when you have the population level of data. Thanks!


title: Tukey's HSD after Repeated Measures One-Way ANOVA

text: We compare the performance of >2 classifiers. Samples are drawn from n-fold cross-validations on the same data sets. We are aware of the fact that cross- validation makes the samples dependent, because training sets overlap. Nevertheless, this is not our main concern here. Because the same data set is used for the cross-validation of each classifier, we decided for a test capable to handle paired samples: the repeated measures one-way ANOVA. Now we wonder if Tukey's HSD can be used as a post-hoc test here. Up to now, we are using paired t-tests for multiple comparisons with p-value adjustments.


title: Adaboost feature weight calculation

text: I thought I understood Adaboost, until code analysis made me realize that sample_weight is not an array of the feature weights... and after further investigation I am left confused as to how multidimensional features are even utilized in Adaboost without stumps or tree. Is the implicit definition of a weak learner the inability to calculate the specific feature weights of a sample? (ie if a learner knew each feature weight, it could simply attempt the classification on its own?) Therefore, an average weight-sum over a set of weak learners gives a sign that can (hopefully) approximate away any mistakes... If so, after delving into the code and a little theory, I am left confused as to how multidimensional features are utilized by weak learners in pure Adaboost. I am interested in knowing the feature weight dynamics as they change over the iterative set of weak learners. I do(believe) I understand Gini and entropy measures method to split/feature selection for trees. I am also under the impression that each feature is utilized in the iterations of Adaboost, perhaps because they are distribution(sum=1)... or are features selected iteratively as well... if so I see no evidence or mention in theory/code code applications I have reviewed. Perhaps there are some helpful resources? So to recap, is my intuition about feature weight complexity vs. approximation of an averaged sample_weight threshold correct? How does sample_weight relate to feature weight? When exponential loss is described as a population minimizer, is this in reference to features, or samples with negative labels? Finally, I apologize if these questions are naive, or have been explained elsewhere. I am really interested in weight dynamics more so than classification, and do not wish to go down the tree path... how are feature weights obtained/tracked using adaboost with exponential loss? Thankyou


title: Model selection in Gamma GLM with identity link

text: I have used Gamma distribution in my GLM with identity link. I want to do backward model selection with drop1 command in R as described in Zuur, 2009. But I get this warning message:               Error: no valid set of coefficients has been found:       please supply starting values             In addition: Warning message:             In log(ifelse(y == 0, 1, y/mu)) : NaNs produced      When using log link it works fine. Why am I getting this message when using Identity link? What is the best option for my model selection? Thank you!


title: Mean structure and mean/variance relationship in regression

text: I read multiple times that the it is crucial in regression to (1) look at the "mean structure", and (2) at the "mean/variance relationship". Alas, I have absolutely no idea what that means, and how it is connected to the assumptions of various regression (e.g. homoscedasticity and normal distribution in OLS, proportional odds assumption in ordered regression, etc.) Could you hint me at some readings?


title: ARMA conditional density

text: Consider the ARMA(1,1) process \begin{align} y_t=a_1y_{t-1}+b_1\epsilon_{t-1}+\epsilon_t, \end{align} and assume $\epsilon_t$~$N(\mu_t,\sigma_t)$. And $\mu_t$ and $\sigma_t$ are all known. $f(y_t|y_{t-1})$ and $f(y_t>c|y_{t-1})$ are easy to calculate. For instance $f(y_t|y_{t-1})$ is the Normal density with mean $a_1y_{t-1}+b_1\epsilon_{t-1}$ and variance $\sigma^2_t$. But how do I calculate $f(y_t|y_{t-1}>c)$ and $f(y_t>c|y_{t-1}>c)$? Many thanks!


title: Hierarchical linear regression analyses, but with logarithms?

text: I was ask to do a linear regression analyses of my variables. I have Reaction times. Can you recommend me a linear regression analysis (hierarchical?) and tell me whether I should transform my raw data into logarithms before that analysis? My reaction times are in milliseconds.


title: VIsualizing the effect that only one predictor has on the outcome (R)

text: So I have performed a logistic regression on a data set with multiple predictors. I want to graphically represent the relationship between the outcome and only one of the predictors. What would be the best way of doing this? I am having a hard time figuring out how to graphically show the effect that predictor X has on Y independent of the other predictor terms. I read something about partial residual plots but I do not know how to implement this in R or how it would apply to data with a binary outcome. Any input is appreciated. Thanks,


title: Which analysis for a set of (0/1)binary variables alone?

text: I have a dataset I would like to analyze and plot It consists of 100 binary variables (0/1) for about 2,000,000 observations There is absolutely no quantitative variable, nor anything I could use as an explained variable for a regression analysis. Actually, the dataset represents the patronage of 2 billion customers for 100 stores. It equals 1 if the consumer go to the store, 0 if he doesn't. With no further information. Consumers can of course visit several stores. As the variable look like factors (0/1), I thought I could go for a Mutliple Correspondence Analysis (MCA). However, the resulting plot consists of 2 points for each variable (one for 1 and one for 0) which is not easily interpretable. (or is there a method for not plotting certain points in MCA? - in R) I also tried to consider my dataset as a bipartite network (consumer-store). However, the plot is not really insightful, as I am especially looking for links between stores. (kind of "if a consumer go to that store, he probably also goes to this one...") So, I have a simple question: which method you would choose for computing and plotting the links between a set of binary variable?


title: Defining Empirical Risk Minimization

text: I am reading Machine Learning - A probabilistic Perspective by Kevin Murphy and in chapter 6.5 the author discusses Empirical Risk Minimization, and provides the following definition: $p_*(\mathbf{x}, y) \approx p_{\text{emp}} (\mathbf{x}, y) = \frac{1}{N} \sum^N_i \delta_{\mathbf{x}_i}(\mathbf{x})\delta_{y_i}(y)$ and then defines the empirical risk as follows: $R_{\text{emp}}(D, D) = R(p_{\text{emp}}, \delta) = \frac{1}{N} \sum^N_i L(y_i) \delta(\mathbf{x}_i)$ I am having a hard time understanding the two formulas above, but mostly, what is the difference between $\delta_{\mathbf{x}_i}(\mathbf{x})$, $\delta_{y_i}(y)$, and $\delta ({\mathbf{x}_i}) $?


title: Analyse unbalanced repeated measures 2x2x2x2 type II anova interactions

text: I have a run an unbalanced 2x2x2x2 Type II ANOVA in R and am having trouble following up on the results. Here is the output:                      Effect DFn DFd             F           p p<.05             ges     2           cond   1 127  3.2359349424 0.074414031       0.0110769653158     3             sf   1 127  1.6981345415 0.194889782       0.0058436648717     5             ba   1 127  1.5404865586 0.216833055       0.0012264293759     9             tt   1 127  1.9253260611 0.167700584       0.0054448755666     4        cond:sf   1 127  0.1846599042 0.668127012       0.0006387833954     6        cond:ba   1 127  5.8799698820 0.016721105     * 0.0046651103251     7          sf:ba   1 127  1.4992638464 0.223051114       0.0011936498636     10       cond:tt   1 127  0.5266890712 0.469337439       0.0014954062256     11         sf:tt   1 127  0.0768302867 0.782090431       0.0002184199961     13         ba:tt   1 127 11.5004885802 0.000927851     * 0.0087996011237     8     cond:sf:ba   1 127  0.0042138896 0.948344162       0.0000033589171     12    cond:sf:tt   1 127  0.1197411309 0.729888009       0.0003403692520     14    cond:ba:tt   1 127  0.3878677814 0.534539033       0.0002993221940     15      sf:ba:tt   1 127  0.0001339682 0.990783282       0.0000001034158     16 cond:sf:ba:tt   1 127  0.4820119706 0.488780454       0.0003719473651      **cond** and **sf** are between-subjects factors. **ba** and **tt** are within (or repeated). The unbalanced nature of the experiment has meant that I have used a type II anova. You can see that we have a suggestively significant main effect of **cond** (2) and two significant interactions (6 & 13). I have graphed the interactions and they seem logical, but I am sure that **cond** is contributing somewhat to the interactions. I am at a loss at how to proceed. I suppose I wish to do some kind of post-hoc analysis concentrating on the interactions. I have investigated a number of different R packages (afex, phia, contrasts etc), but have yet to work out what I am actually doing with these interactions. My data looks like this:               str(xx)     'data.frame':   524 obs. of  6 variables:     $ ba  : Factor w/ 2 levels "before","after": 1 1 1 1 1 1 1 1 1 1 ...         $ tt  : Factor w/ 2 levels "targ","calm": 1 1 1 1 1 1 1 1 1 1 ...     $ p   : Factor w/ 131 levels "1","2","3","4",..: 4 8 9 10 13 18 19 22 25 29 ...         $ cond: Factor w/ 2 levels "Control","Spider": 1 1 1 1 1 1 1 1 1 1 ...     $ sf  : Factor w/ 2 levels "Fear","No-Fear": 1 1 1 1 1 1 1 1 1 1 ...         $ eda : num  1.478 -0.56 -0.27 -0.902 -0.483 ...      Moreover consider Thompson (2006) : > As noted by Rosnow and Rosenthal (1989a), the cell means “are the combined > effects of the interaction, the row effects [a main effect], the column > effects [a second main effect], and the grand mean” (p. 144). By the same > token, simple post hoc tests of the cell means also do not yield insight > about the origins of interaction effects, because the interaction effects > are not uniquely a function of the cell means (Boik, 1979). So I guess post-hoc t-tests (with adjusted p-values) are out? **Update:** Following advice I have _found_ from @henrik (here) I have been investigating the phia package and the **testInteractions** function. I have been getting some results (for a type III anova - so not the type II I am after) but, again, I am way out of my depth here: e.g.,               > testInteractions(m2[["lm"]], pairwise = "ba", "cond", idata = m2[["idata"]],     adjustment = "none")                   Multivariate Test: Pillai test statistic     P-value adjustment method: none                           Value Df test stat approx F num Df den Df  Pr(>F)       after-before : Control -0.31910  1  0.045305   6.0268      1    127 0.01544 *     after-before :  Spider  0.14243  1  0.008045   1.0300      1    127 0.31208       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       **Latest Update:** I seem to be getting further bogged down with this. Using lmer as referenced here gives this:               fit.1 <-lmer(eda ~ 1 + cond * sf * ba *tt +(ba *tt | p), xx)      which leads to:               > anova(fit.1)     Analysis of Variance Table               Df  Sum Sq Mean Sq F value     cond           1 0.00330 0.00330  0.0158     sf             1 0.04236 0.04236  0.2035     ba             1 0.07408 0.07408  0.3559     tt             1 0.34651 0.34651  1.6649     cond:sf        1 0.05633 0.05633  0.2707      cond:ba        1 2.38017 2.38017 11.4359     sf:ba          1 0.45349 0.45349  2.1789     cond:tt        1 0.03832 0.03832  0.1841     sf:tt          1 0.03197 0.03197  0.1536     ba:tt          1 2.39403 2.39403 11.5025     cond:sf:ba     1 0.02595 0.02595  0.1247      cond:sf:tt     1 0.00667 0.00667  0.0320     cond:ba:tt     1 0.08078 0.08078  0.3881     sf:ba:tt       1 0.00003 0.00003  0.0001     cond:sf:ba:tt  1 0.10034 0.10034  0.4821      Which seems to be confirming the earlier finding of a significant interaction between cond * ba and ba * tt (but is it type II?). Still not sure if any of this is correct. **To clarify then: I am looking for advice about what to do next, in terms of understanding the significant interactions.**


title: How to report hazard ratios from a Cox proportional hazards model in English?

text: My understanding is that a hazard ratio from a Cox proportional hazards model compares the effect on the hazard rate of a given factor to a reference group. How would you report that to an audience that doesn't know statistics? Let's try to phrase an example. Say we enroll people in a study of how long before they buy a couch. We right-censor at 3 years. For this example we have two factors: age < 30 or >= 30, whether they own a cat. It turns out the hazard ratio of "owns cat" to the reference group (age < 30, "doesn't own cat") is 1.2, and significant (say p<0.05). Am I correct to say that means all of these: cat owners have more events (couch purchase) within 3 years, OR that time-to-event (couch purchase) is faster for cat owners, OR some combination of those two things? **Edit** : Suppose the event is their first purchase of the couch within the period (if one occurs). This model does not help us analyze multiple purchases within the time period.


title: Dividing controls between two groups of cases

text: I have a set of data that contains cases and controls. I would like to divide the cases by a particular phenotype before running my tests, to see how these groups differ (but not necessarily to directly compare them). My feeling is that I must also randomly divide the controls between these two groups of cases before I go any further. Is this true? I've searched a bit, but I think this may be such an elementary thing that I can't find it spelled out anywhere. Would appreciate your thoughts and expertise.


title: Measure to use in experiment

text: Suppose we have $100$ people and we obtain the average blood pressure $x$. Now suppose we have another group of $40$ people and they have an average blood pressure $\approx x$. How do we quantify the fact that the second experiment is more "efficient" than the first experiment? Would we do a power analysis? What would be the null hypothesis?


title: Random walks in multinomial case

text: **Model** : A vector $X=(X_1, X_2, X_3)$ that follows a trinomial distribution with parameters $p=1/3$ and $n$. (I have a coin with three sides $S1$, $S2$, $S3$). I flip the coin $n$ times. The coin has a probability $p=1/3$ to be flipped to the side $S1$, and similarly with $S2$ and $S3$. $X_1$ counts the number of times the coin is flipped to $S1$ ($X_2$ and $X_3$ are defined similarly). **Questions** :   1. Let $0 \leq \alpha \leq n$. I Want to find $T(\alpha)= P (((X_1 - X_2) \geq \alpha) \cap ((X_1 - X_3) \geq \alpha))$, or a lower bound on $T(\alpha)$   2. For what values of $\alpha$, The lower bound on $T(\alpha)$ does not depend on $n$ (is a constant) ? If $X$ was following a binomial distribution, the problem would have been easy to solve using random walks, but I do not know how to solve it in the multidimensional case. Any idea? Thank you. **Two-dimensional case** $X= (X_1, X_2)$ follows a binomial distribution with parameters $p=1/2$ and $n$. (A coin is flipped $n$ times. the coin is up with probability $p=1/2$ and down with probability $p$. $X_1$ counts the number of ups and $X_2$ counts the number of downs) Let $Y= X_1 - X_2$ (We have $n$ coin flips, when the coin is flipped up, we add (+1) to $Y$, when the the coin is flipped down, we add (-1) to $Y$). We can see this as a random walk, when the coin is up we go to the right and when it is down we go to the left. $\mathrm{Var}(Y)=n = \sqrt n ^2 $. $E(Y)=0$ By the Central Limit Theorem, $P(Y \geq \alpha) \approx 1 - \Phi(\alpha / \sqrt n)$ (normal distribution).


title: What are latest statistical challenges regarding missing data?

text: Could you suggest more recent reviews or book on the methods for handling missing data (i.e. multiple imputation)? I want to know the challenges in this field.