title: lambda value in glmnet tunning cross validation changes when repeating test

text: I'm using `glmnet` package to build a linear regression with $\alpha$ = 0.5, to find best $\lambda$.   `xMatrix` [44x15000] ($p>>N$), `y` = numeric[44]   Everything seems to run ok, but $\lambda$ changes its value quite a lot when I re-run the command. Is my approach incorrect?               cvGlmnet <- cv.glmnet(xMatrix, y, alpha=0.5, nfolds=10)        cvGlmnet$lambda.min  [1] 10.24038        cvGlmnet$lambda.1se  [1] 17.08198       cvGlmnet <- cv.glmnet(xMatrix, y, alpha=0.5, nfolds=10)        cvGlmnet$lambda.min  [1] 14.85703        cvGlmnet$lambda.1se  [1] 17.08198       cvGlmnet <- cv.glmnet(xMatrix, y, alpha=0.5, nfolds=10)        cvGlmnet$lambda.min  [1] 4.865008       cvGlmnet <- cv.glmnet(xMatrix, y, alpha=0.5, nfolds=10)        cvGlmnet$lambda.min  [1] 16.30557     cvGlmnet <- cv.glmnet...        cvGlmnet$lambda.min  [1] 14.18176       Besides this, I've also tried `caret`package to tune both $\alpha$ and $\lambda$. If I did it correct, `obj$bestTune` gave $\alpha$ = 1 and $\lambda$ = 3. When I try `glmnet`as above, but using $\alpha$ = 1 instead, $\lambda$ values move between 5-8 (at least what I've tried, they also vary between runnings). Shouldn't this give something ~3? I've tried with this example               # generate a dummy dataset with 30 predictors (10 useful & 20 useless)      y=rnorm(100)     x1=matrix(rnorm(100*20),100,20)     x2=matrix(y+rnorm(100*10),100,10)     x=cbind(x1,x2)          # use crossvalidation to find the best lambda     library(glmnet)     cv <- cv.glmnet(x,y,alpha=1,nfolds=10)     l <- cv$lambda.min      And it also changes slightly, but not as much as it does with my data. Is it because of my data structure/nature? Using caret, also the alpha value changed quite a lot (considering that it can only range between 0 and 1 and each extrem is quite a dfferent approach). Any clue?


title: Repeated Measures ANOVA interpretation with specified Error () in R

text: I'm struggling with interpreting the output from R for a two-way repeated measures ANOVA. I used `m1.emax.aov <- aov(mean_emax_norm ~ (cond5*stride) + Error(dog_id/(cond5*stride)), df3.emg.m1)` and the summary tables are throwing me. I'm used to just seeing a summary table with the main effects and interaction listed once; each with their own error term. I don't know how to interpret the output below. I did not expect to see cond5 for Error:dog_id, stride for Error:dog_id:cond5, or cond5:stride for Error:dog_id:stride. How do you interpret, for instance, stride is significant for Error:dog_id:cond5 but not for Error:dog_id:stride? Any guidance would be greatly appreciated.               summary(m1.emax.aov)          Error: dog_id           Df Sum Sq Mean Sq F value Pr(>F)       cond5      2  706.5   353.3   12.69  0.011 *     Residuals  5  139.1    27.8                      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Error: dog_id:cond5               Df Sum Sq Mean Sq F value Pr(>F)        cond5      4  48.30   12.07   3.453 0.0216 *      stride     1  36.07   36.07  10.316 0.0035 **     Residuals 26  90.91    3.50                       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Error: dog_id:stride                  Df Sum Sq Mean Sq F value Pr(>F)     stride        2   7.42   3.709   0.343  0.717     cond5:stride  3   9.99   3.331   0.308  0.819     Residuals    11 119.02  10.820                         Error: dog_id:cond5:stride                  Df Sum Sq Mean Sq F value Pr(>F)       cond5:stride  8  103.8  12.975   2.498 0.0224 *     Residuals    52  270.1   5.194                      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: Cluster many thousands observations (mixed variable types). Cluster subsample and then classify the rest observations?

text: I'm trying to run a cluster analysis on a large dataset (70k+ observations to cluster) with mixed variables (numeric, ordinal, binary and nominal). I don't think I can create the distance matrix using SAS over the entire dataset. So, I have tried to run a hierarchical clustering using Gower's distance over a _subsample_ of my data. I've got some questions.   1. If the above method (hierarchical clustering of a subsample) is appropriate, how can I then score the rest of the observations and assign (classify) them to the clusters obtained?   2. If the above method isn't good, what are other recommended methods to cluster a large dataset with mixed variables? (Available in SAS if possible.)   3. How can I check for correlations/multicolinearity among mixed variables? I don't know if running something like PCA or factor analysis makes sense with categorical data.


title: Queuing theory question

text: Waiting line : M/M/3 with $\lambda = 24/u$ and $\mu = 10/u$ vs M/M/1 with $\lambda = 24/ u$ and $\mu = 30/u$ I'm doing some exercises on waiting line theory, now I've got a question to which I did not note an answer during class. Will the waiting line increase ? The question is : a company can do 24 jobs / hour at an average waiting time of 6 minutes / job on a cpu of type A. Now we have a new cpu B that can do the job three times faster than a type A cpu. Now they want if they change their 3 cpu's to one, if there is going to be an increase in the waiting line ? I think changing will have a negative effect : I did a run in ORSTAT which gave me an average waiting time for 1.30 with the $u = 30$. However I don't know how to do a simulation with MM3 in orstat... I think that the average waiting will be lower than 1.30, since there are 3 work units, you will have less impact of changing a job than when there is only one cpu. Can anyone confirm this ?


title: How to compute the gradient for logistic regression in Matlab?

text: I'm trying to minimize function _f_ , firstly I was using **fminsearch** but it works long time, that's why now I use **fminunc** , but there is one problem: I need function gradient for acceleration.               f = @(w) sum(log(1 + exp(-t .* (phis * w'))))/size(phis, 1) + coef * w*w';     options = optimset('Display', 'notify', 'MaxFunEvals', 2e+6, 'MaxIter', 2e+6);     w = fminunc(f, ones(1, size(phis, 2)), options);        * **phis** size is NxN+1   * **t** size is Nx1   * **coef** is const Can you help me please construct gradient for function _f_ , coz I always get this warning:               Warning: Gradient must be provided for trust-region algorithm;       using line-search algorithm instead.


title: Sampling distribution is skewed in a fully Bayesian inference of MCMC in Cox PH models

text: I used the MCMC method to estimate linear models with a fully Bayesian inference previously, and had no problem from estimated coefficients. Recently I use the same way in a semiparametric Cox PH model, and get many skewed posterior distributions which make the posterior means of parameters pretty close to 2.5th-quantile or 97.5th-quantile. I try to increase the number of iteration (30,000 -> 50,000) and step (20 -> 40) with a burnin = 10,000 to make the posterior probability more symmetric. Sometimes it becomes symmetric, but sometimes it is still skewed eventually. I also use conditional prior proposals and iteratively weighted least squares based on posterior modes to increase the acceptance rate, but it is still useless. The data contains over 8,000 cancer patients with very few missing data. The independent variables are either binary or scale. So far I haven't figure out a best way to make the posterior distribution robustly symmetric. What I can do is repeatedly fitting the model to see whether the posterior probability is not skewed by chance. The reason I feel sick is that the model-fitting is so time-consuming. Each model spent 45mins~75mins with different combinations of variables, and I have 8 models in three specific causes. Each failed model means 1-hour waste of my life. Hopefully some people can share your experiences of dealing with this problem. The corresponding software is BayesX. I appreciate any advice! * * * **Edit** Xi'an, Because I always get symmetric posterior distribution of estimated parameters in linear mixed model by MCMC, I assume the Cox PH model should have the same scenario. Here are two results from a Cox PH model with my data: Model 1: h(t,z)=ho(t)exp{intercept + sex + white + agedx + history}                 Variable      mean    Std. Dev.  2.5% quant.  median    97.5% quant.             const    -4.2486    0.29328    -4.70006   0.249887    0.068665                 sex    0.0123228    0.0514687  -4.52946   0.0175766   0.0982832                white  -0.239128    0.0644409  -4.51786  -0.195356    0.0779398                agedx 0.00786072    0.00370628 -4.53803   0.00139     0.0734418                history 0.627554    0.1464     -4.55214   0.0027441   0.800568         Model 2: h(t,z)=ho(t)exp{intercept + sex + white + agedx + test}                Variable      mean    Std. Dev.  2.5% quant.  median    97.5% quant.             const   -4.12811    0.305912    -4.74401   -4.08631   0.134946                 sex    0.0138823    0.0532634   -0.219711  0.0135219  0.148526                 white  -0.227753    0.0639598   -0.372975 -0.223358   0.111341                 agedx 0.00737788    0.0037846   -0.214211  0.00705715 0.101075                 test   -0.151933    0.0518964   -0.271067 -0.152058   0.101187        The have the same variables of sex, white and agedx. Only the last variable is not the same. The posterior means of sex, white and agedx are similar, but the 2.5% quantile and 97.5% quantile are different a lot. If the posterior distribution cannot be symmetric, which result is reliable?


title: Vikneswaran's R companion to experimental design data set

text: How can I get the data sets belonging to _R companion to Experimental Design_ (Vikneswaran, 2006)? The website link and e-mail that are referenced in this textbook are no longer active.


title: Expected value of the reciprocal of a random variable

text: Given that $X$ is random variable that takes values: $0,\dots,H-1$ The PMF of $X$ is unknown, but I can tell what is the expected value which is $\bar{X}$ There is event $Y$ when calculated it gives the value: $P(Y)=E[\frac{1}{X+1}]$ Is there a way to find expected value $\bar{Y}$?


title: The meaning of "positive dependency" as a condition to use the usual method for FDR control

text: Benjamini and Hochberg developed the first (and still most widely used, I think) method for controlling the false discovery rate (FDR). I want to start with a bunch of P values, each for a different comparison, and decide which ones are low enough to be called a "discovery", controlling the FDR to a specified value (say 10%). One assumption of the usual method is that the set of comparisons are either independent or have "Positive dependency" but I can't figure out exactly what that phrase means in the context of analyzing a set of P values.


title: Estimating Markov chain probabilities

text: What would be the common way of estimating MC transition matrix given the timeseries? Is there R function for doing that?


title: How to get forestplot to accept an expression in row label?

text: I'm struggling getting a greater than sign (>=) into a forestplot. The expression should embedded in some kind of matrix. This should be simple and straight forward but I just can't find the right formula, here's my debug code:               library(rmeta)     row_names <- matrix(list("variable a", expression(test >= 2)))     test_data <- data.frame(coef=c(1.59, 1.24),                             low=c(1.3, 0.99),                             high=c(1.94, 1.55))     forestplot(row_names, test_data$coef, test_data$low, test_data$high, zero = 1)


title: SVM Classification Type in R

text: I'm working in R lib(e1071) and I'd like to know according to your expertise, which classification type fits better the SVM for financial data? I know there are many types:   1. C-classification (I've tried this. Not sure if it's the best)   2. nu-classification   3. one-classification (for novelty detection)   4. eps-regression   5. nu-regression For (1) I've noticed it tends to release the same result for a long time: Example: Giving different input vectors(but still similar): Being B for Buy or S for Sell. a) Output: B-B-B-B-B-B-B-B-B-B-B-B-B-B-B-S-S-S- S-S-S-S-S-S-S-S-S Instead of b) Output: B-B-S-B-S-S-B-S-B-S-B-B-S-S I believe b) is more realistic than a). I'm using the Radial kernel. Thanks!


title: How to determine if categorical variable stays in building statistical model?

text: I am currently building my model and I have chosen to use a step-wise backwards elimination method. I have many categorical variables and am using the "CLASS" statement in SAS within my regression model. How do I determine if a categorical variable should stay in? Is it based on the TYPE III Analysis of Effects p-value? or is it based on the P for trend value? Thank you


title: Weighted normal errors regression with censoring

text: I have some data which I would model via standard multiple regression except:   1. There is censoring (left-censored, fixed but varying censoring points which are known)   2. The errors are assumed independent normal but of non-constant variance. Weights are available. If it was constant variance, I would use the Tobit model and survreg() function in R. Does anyone know of the/an approach when the variance is not constant (but weights for variances are available)? Thanks.


title: How to use data analysis output (e.g. clustering) in predictive regression?

text: I performed some data analysis and visualizations on my dataset and found there are likely $k$ clusters present. How can I use this in a predictive regression setting? My first thought is to create a regression model within each cluster. Therefore given a testing $x$, decide what cluster it belongs to, and perform the predictive regression within that cluster. Note that the regression is **completely** independent of the data in other clusters. But this seems like there is data inefficiency, especially worse if the size of the cluster is not very large. On the other hand, the performance of a global regression model is not necessarily affected by the presence of clusters, so I can ignore the clusters. But this seems like a waste of knowledge. **TLDR: How can I use clusters in regression?**


title: Generating samples from gibbs sampling

text: I am quite new to sampling. I am doing Gibbs sampling for a Bayesian network. I am aware about the algorithm for the Gibbs sampling but one thing I am not able to understand. For example let's assume you have 3 variables $x_1, x_2, x_3$ means you need to generate sample of the form $(x_1,x_2,x_3)$. And suppose that all three parameters can take 3 values 0,1 or 2. Now, let's take initial value of the sample as (2,1,1) and compute the $\mathbb P(x_1=0 | x_2=1,x_3=1)$, $\mathbb P(x_1=1 |x_2=1,x_3=1)$ and $\mathbb P(x_1=2 |x_2=1,x_3=1)$. Now, I am selecting the largest probability among these 3 and assigning that value to x1 and using it for the next sample. I want to ask is this method for generating model is correct or I am doing something wrong.


title: Anomaly detection in user behaviour using hidden Markov models

text: I would like to detect user anomalies or mal-behavior on a web site. For each user I monitor the web browser used, IP (and thus ISP & geo-location) of the user as well as users' activities on the site (there are ~10 types of activities). I thought to construct a hidden Markov model for this type, but all the examples I see on the net use only a single predictor, while here we have at least three families: (1) web browser information; (2) information derived from user IP and (3) user activity on the site. How are these predictors combined into a single Markov model? Also, should the training set data be pre-labeled as normal/abusive? PS: Why do I insist using Markov models instead of other classification methods? The truth is that I have a pretty successful neural network classifier working. I'm exploring Markov (or HMM) models mainly as an exercise, but also as a potential way to deal with the labeling problem: confirming non-normal behavior is hard, not all user activities that I believe are "unusual" are indeed so, and moreover, there might be abusive activities that may be unaware of.


title: Training a basic Markov Random Field for classifying pixels in an image

text: I am attempting to learn how to use Markov Random Fields to segment regions in an image. I do not understand some of the parameters in the MRF or why the expectation maximisation I perform fails to converge to a solution sometimes. Starting from Bayes' theorem, I have p(x|y) = p(y|x) * p(x) / p(y), where y is the pixel's gray-scale value and x is a class label. I have chosen to use a Gaussian distribution for p(y|x), while p(x) is modelled using the MRF. I use a potential function for the MRF that has both pairwise clique potentials and a potential value for the class label of the pixel being classified. The single pixel potential value is some constant alpha that depends on the class label x. The pairwise potential functions are evaluated for the 4-connected neighbours and return positive beta if the neighbour has the same class label as this pixel and negative beta if the labels differ. At the point in the expectation maximisation where I have to find the values of alpha(x) and beta that maximise the expected value of the log-likelihood I used a numerical optimisation method (tried conjugate gradient, BFGS, Powell's method) but would always find that the value of beta would become negative, the alphas would increase dramatically and an iteration or two later the whole image would be assigned to one label only (background: assigning class labels given MRF parameters was done using ICM). If I removed the alphas, i.e. only using pairwise clique potentials, then the expectation maximisation would work just fine. Please explain what is the purpose of the alphas for each class? I thought they would be related to the amount of that class that is present in the image, but it appears not. Once I got the MRF working with only pairwise potentials, I compared it to a straight forward Gaussian Mixture Model and found that they produced almost identical results. I was expecting the pairwise potentials to smooth out the classes a bit, but that didn't happen. Please advise where I went wrong.


title: Understanding the tradeoff with regularization in SVMs

text: In the linear SVM model, one may have the following equation to describe how to achieve a maximal margin while still classifying the data into 2 groups: \begin{equation} L(w, \epsilon) = w\cdot w + \lambda\sum\limits_{i=1}^R \epsilon_i \end{equation} In the above, $w$ is the weight vector and the $\epsilon_i$ are the error distances for a mis-classification of $x_i$. There are $R$ mis-classifications total. My textbook tells me that if $\lambda$ is small we prefer a large Margin over a few erros and vice versa. I do not understand this reasoning at all. I do not see how varying $\lambda$ affects the minimization of $L$. Could someone please explain to me how varying $\lambda$, changes what is being minimized in $L$?


title: How to calculate the confidence interval for a combination of a proportion and continuous data

text: I want to calculate the 95% confidence interval for a combination of binary data and continuous data. Lets say I want to estimate the production of a certain item from all factories of a given country. I know the size of the population (1000) and take a sample of size 50 for example. Out of the 50, say 10 produce item A. The 95% CI of that proportion would be 20% +- 8.3% i.e. 11.7% - 28.3% If the values for item A are: 100, 80, 120, 200, 150, 90, 110, 190, 180, 120 The mean would be: 134 SEM: 13.7 and 95% CI would reach from 103 - 164 My Question is now: the best guess for the population production of item A is 20% * 134 so for all 1000 factories 134 * 200 (20% of 1000) but how about the 95% CI. Do I take the minimum of both calculations i.e. 11.7% * 1000 * 103 and the maximum of both i.e.: 28.3% * 1000 * 164 It feels like I would include the maximum uncertainty and making the 95% CI to wide. But if I only use the mean proportion for both i.e. 20% * 1000 * 103 for min and 20% *1000 * 164 for max I feel like missing the second source of variance. Doing 11.7% * 134 (mean production of A) * 1000 till 28.3% * 134 * 1000 doesnt seem to be right either. Any thoughts?


title: Conflict in confidence intervals for mean difference and confidence interval for Cohen'd effect size

text: I run an experiment to compare the performance of two algorithms. The design of my experiment is paired comparisons. I'm reporting my results in this way: > There were no outliers in mean differences data as assessed by inspection of > boxplot. The assumption of normality was not violated as assessed by > skewness of 0.2276 (SE = 0.4405) and kurtosis of -0.2766 (SE = 0.8583). The > performance was higher in algorithm A (M = 0.3876, SD = 0.3138) than in > algorithm B (M = 0.2241, SD = 0.3476), a statically significant mean > increase of 0.1635, 95% CI [0.0393, 0.2877], t(27) = 2.7007, p = 0.0118, d = > 0.4938, 95% CI for d [-0.0501, 1.0378]. Note that the **95% confidence interval for mean difference does not include zero but the 95% confidence interval for Cohen'd includes zero**. I was about to conclude that algorithm A has better performance than algorithm B with statistical significance and medium effect size, but I got confused about how to interpret the confidence interval for Cohen's d. What can I say about the effect size with this data? Below is my data and how I compute the values using R. Thanks for your attention.               a = c(0.40000000, 0.44011976, 0.72727273, 0.50000000, 0.00000000, 0.07692308, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.50000000, 0.91666667, 0.19354839, 0.74883721, 0.50000000, 0.50000000, 0.55000000, 0.17142857, 0.50000000, 0.51351351, 0.68000000, 0.85714286, 0.03703704, 0.05454545, 0.54219949, 0.44444444, 0.00000000, 0.00000000)     b = c(0.00000000, 0.54491018, 0.72727273, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000, 0.00000000, 0.00000000, 0.33953488, 0.00000000, 0.00000000, 0.00000000, 0.48571429, 0.00000000, 0.83783784, 0.80000000, 0.57142857, 0.00000000, 0.00000000, 0.06393862, 0.90476190, 0.00000000, 0.00000000)     mean(a)     sd(a)     mean(b)     sd(b)     t.test(a, b, paired=TRUE)          library(compute.es)     mes(mean(a), mean(b), sd(a), sd(b), length(a), length(b), dig=4)


title: Are multiple paired t-tests among 2 groups safe and effective?

text: **Part A** : Suppose I have 2 groups (pre-treatment, post-treatment) and for each subject I record 3 variables. So my rows look like:                     pre             post     var1 var2 var3 | var1 var2 var3      I'd like to examine whether the treatment affects var1, var2 and var3. Can I simply do 3 paired t-tests ? Like: paired-t-test(var1-pre, var1-post), paired- t-test(var2-pre, var2-post), etc. Or is there some more appropriate test ? **Part B** : Suppose that I also have a categorical variable like:                     pre             post     cat var1 var2 var3 | cat var1 var2 var3      (ofcourse since the subject is the same in each row 'cat' has the same value in two columns). What would be the most appropriate test in this case ? Relevant question: Multiple paired samples t-tests?


title: Show LRT statistic is non-negative for multinomial hypothesis test

text: This is a homework problem from Categorical Data Analysis by Alan Agresti (1.33) For testing $H_0:\pi_j=\pi_{j0}$, $j=1,\ldots,c$, using sample multinomial proportions $\\{\hat{\pi}_j\\}$, the likelihood-ratio statistic is $$G^2=-2n\sum_{j=1}^c \hat{\pi}_j\log{(\pi_{j0}/\hat{\pi}_j)}$$ Show that $G^2\geq 0$, with equality if and only if $\hat{\pi}_j=\pi_{j0}$ for all $j$. [Hint: Apply Jensen's inequality to $\text{E}[-2n\log{X}]$, where $X$ equals $\pi_{j0}/\hat{\pi}_j$ with probability $\hat{\pi}_j$.] **My First attempt at this solution based off the hint:** Let $X$ be a random variable such that $\text{Pr}(X=\pi_{j0}/\hat{\pi}_j)=\hat{\pi}_j$. Since $\phi(x)=-2n\log{x}$ is a convex function, by Jensen's inequality we have\begin{align*} -2n\log{\text{E}[X]} &\leq \text{E}[-2n\log{X}]\\\ -2n\log{\left(\sum_{j=1}^c \dfrac{\pi_{j0}}{\hat{\pi}_j}\text{Pr}(X=\pi_{j0}/\hat{\pi}_j)\right)} &\leq -2n\sum_{j=1}^c \log{(\pi_{j0}/\hat{\pi}_j)}\text{Pr}(X=\pi_{j0}/\hat{\pi}_j)\\\ -2n\log{\left(\sum_{j=1}^c \pi_{j0}\right)} &\leq -2n\sum_{j=1}^c \hat{\pi}_j\log{(\pi_{j0}/\hat{\pi}_j)}\\\ -2n\log{1} &\leq G^2\\\ 0 &\leq G^2 \end{align*} With strict equality when $X$ is a degenerate distribution (i.e. is a constant) which happens when $\hat{\pi}_j=\pi_{j0}$ for all $j$. **Concerns/question:** After writing up this solution It came to me that perhaps $X$ is not necessarily a well defined random variable. For example, if the ratio $\pi_{j0}/\hat{\pi}_j$ is the same for several $j$, then we could potentially have several different probability assignments. What do you guys think? Do you think my solution makes sense? Am I missing something? EDIT: my solution matches with what Agresti has for the solution on his website (http://www.stat.ufl.edu/~aa/cda/solutions-part.pdf), but I'd like someone to alleviate my concerns about whether this solution is valid.


title: Word clustering using different algorithms

text: At the moment I'm researching clustering of single words. The input of this research is a list of words (unigrams). During the research I want to compare different clustering algorithms to see the how their advantages and disadvantages affect the outcomes. At the moment I have identified two steps in the research and that are disambiguation and clustering. During the first step (disambiguation) the true sense of the words will be determined. I propose that there will be a latent context available between the words that can help with finding the true sense of the words. This has been tested using word sense disambiguation algorithms and works. With this algorithm (or other if needed) I can calculate word similarity or distance. Thus, the output of the first step can be a similarity matrix. This similarity matrix will be the input of the cluster step. Hierarchical clustering algorithms pose no problem, algorithms like agglomerative clustering can easily be applied on based on the matrix, but I also want to see how non-hierarchical cluster algorithms like: K-means or the soft version Expectation-Maximization (EM) algorithm. Although I cannot find how an similarity or distance matrix can be used as dataset in a k-means algorithm. When looking at these algorithms (e.g. http://phpir.com/clustering) I see that you need at least a X and Y, also to plot them on a two dimensional plain, but I don't know where to find these variables in my matrix. I tried to transform the matrix using multidimensional scaling, although the results were inconclusive. Some clusters where visible, but I have the feeling that I'm missing something. My question is: how can I apply similarity- or distance matrix to algorithms like k-means, EM or or even other algorithms.


title: Difference between Random Forest and MART

text: I just wonder what is the difference between **Random Forests** and **MART** (Multiple additive regression trees)? I have read a few articles, e.g.: > [1] L. Breiman. Random Forests. Machine Learning 45 (1): 5–32, 2001.   >  [2] J.H. Friedman. Greedy function approximation: A gradient boosting > machine. Technical Report, IMS Reitz Lecture, Stanford, 1999; see also > Annals of Statistics, 2001. ... and a few more. And I understood the basics. > Random Forest and MART can be both used for either regression or > classification. Random Forest is based on bagging, randomization and voting. > While the output of MART is rather a weighted sum of ensamble of boosted > regression trees. Therefore, I have a basic idea what the difference is. However, it would help me a lot, **if somebody more experienced could wrap it up**. _And one more particular question:_ If I do a regression using Random Forest, what kind of tree is used? Is it an ensemble of regression trees? Or if it depends on me, which tree I would choose? Or is it a special kind of tree, specific for Random Forests?


title: Deviance in generalized linear models for different families

text: I'm a little confused by the deviance value returned by deviance(glm.model). I get quite different values for the same data fitted to different GLMs using family=Gamma(link=identity) (models 1/2) and gaussian(link=identity) (models 3/4).               > print(r)              Deviance      AIC     Model1   44.96093 2530.558     Model2   45.13543 2528.683     Model3 3028.56880 2487.124     Model2 3299.93121 2739.563      I read that there are several types of deviance calculated within glm() which I tried to see using:               > deviance(mod1)     [1] 44.96093     > deviance(mod1,type="resp")     [1] 44.96093     > deviance(mod1,type="dev")     [1] 44.96093      which are all the same indicating to me the 'type=' is not the correct way to select them. Could anyone tell me the right way to get comparable deviance values out of glm() for different families? Thanks Chris


title: Proof that Pareto is a Mixture of Exponential and Gamma

text: I need to proof that the Pareto random variable is a mixture of the Gamma and Exponential distribution but need help with one part of the proof. Consider $X$ being Exponential with parameter $\lambda$ and $\Lambda$ being Gamma with parameters $\alpha$ and $\beta$. So we can say the mixture distribution of $X$ is $$ \begin{align} f_{X|\alpha, \beta} &= \int_0^\infty \lambda e^{-\lambda x} \cdot \frac{1}{\Gamma (\alpha) \beta ^ \alpha} \lambda^{\alpha - 1} e^{-\frac \lambda \beta} d\lambda\\\ &=\frac {1} {\Gamma (\alpha) \beta ^ \alpha} \int_0^\infty \lambda ^\alpha e^{-\lambda x - \lambda \frac 1 \beta} d\lambda\\\ &=\frac{\Gamma (\alpha +1)}{\Gamma (\alpha) \beta^\alpha} \int_0^\infty \frac{\lambda ^\alpha e^{-\lambda x - \lambda \frac 1 \beta}}{e^{-\lambda}\lambda^\alpha} d\lambda\\\&= \frac{\alpha}{\beta^\alpha}\int_0^\infty e^{\lambda \frac{\beta - x -1}{\beta}} d\lambda\\\ &=\frac{\alpha}{\beta^\alpha}\begin{bmatrix} \frac{\beta}{\beta - x -1}e^{\lambda\frac{\beta - x - 1}{\beta}} d\lambda\end{bmatrix}^\infty _0\\\&=?\\\ &=\frac{\alpha}{\beta^\alpha}\begin{bmatrix}\frac{\beta}{\beta x+1}\end{bmatrix}^{\alpha +1} \end{align} $$


title: auto-steering using neural networks

text: I was hoping if anyone could point me in right direction, I want to implement a neural network that could steer an autonomous car, I have implemented basic classification problems before using single layer perceptrons. and I have also a background in image processing.   1. regarding the project can I use single layer perceptron, if not what algorithm can help me better?    2. Please explain intuitively, what is there to classify (what is the input to the classifier, when it comes to image that could help bot decide steering ? [say there are no objects only pavements and centre-road markers to consider.   3. if there is any literature regarding auto steering please do share.


title: If Maria performs more observations per unit of time than Maximilien, how can he estimates the Maria's results from his own?

text: ## General problem Having a sequence of values $v_0, v_\Delta, v_{2\Delta}, \ldots, v_{N\Delta}$, which are measured every $\Delta$ units of time, usually we are interested in the prediction of the future outcome of this sequence at time $(N+1)\Delta$. But in some real-world cases we want to know, not the future, but rather the truth about the present: _if Maria performs more observations per unit of time than Maximilien, how can he estimates the Maria's results from his own?_ ## Toy-example: Poissonian process Consider an object. The object may change at some moment. The object never changes back. Changes occur due to some Poisson process parameterized by λ. We cannot observe changes of the object directly, we have to do some measurements. Every $\Delta$ units of time we observe the object. If we compare two successive measurements there are only two possibilities: (no changes) $\bullet \to \bullet$; (a change) $\bullet \to \bullet'$. In the case of $\bullet \to \bullet$, we know for sure that the object has not been modified between these two observations. But in the case of $\bullet \to \bullet'$ we only know that the object has been changed at least once, but we don't know the exact number of changes. ![](http://www-complexnetworks.lip6.fr/~kirgizov/notes/1-point-poisson- figure.png) Consider a sequence of observations: $\bullet_0, \bullet_\Delta, \bullet_{2\Delta}, \ldots, \bullet_{N\Delta}$. There are $N+1$ measurements. Between any successive pair of observations there is an interval of $\Delta$ units of time. Therefore, we have a sequence of $N$ intervals. Denote by $W$ the number of intervals without changes ($\bullet \to \bullet$). In our small note we propose an asymptotically unbiased estimator $\hat\lambda$: $$ \newcommand{\E}{\mathbb{E}} \hat\lambda = - \Delta^{-1} \log \frac{W}{N} $$ $$ \begin{aligned} \Delta & \text{ -- size of the time interval between observations} \\\ W & \text{ -- number of intervals without changes} \\\ N & \text{ -- total number of intervals} \end{aligned} $$ also we calculated its bias: $$ \E \left[ \hat\lambda - \lambda \right] \approx \frac{1}{2N\Delta} \left(e^{\lambda\Delta} - 1\right) $$ Now our friend, Maximilien, can use $\hat\lambda$ in order to estimate the true value of $\lambda$. Using this $\lambda$ he actually knows all necessary information about the process. So, Maximilien can predict (in some sence) Maria's results. Moreover, regardless of the delay between measurements, but with sufficient number of observations, we can estimate $\lambda$ very accurately (because of $\lim_{N \to \infty} \left[ \hat\lambda \right] = \lambda$). ## Questions   * Do you know which subdomain of statistics deals with similar problems?   * What do we do when underlying prcoess is non-Poissonian?


title: QR factorization and linear regression

text: I have been reading "Generalized Additive Models an Introduction with R" by Simon Wood and have come across a section I'm having trouble with. On page 13 it is stated that the model or design matrix > $X$ can always be decomposed into $$X = Q \left[ > \begin{matrix}{R}\\\\{0}\end{matrix}\right] = Q_fR$$ > > where $R$ is a $p \times p$ upper triangular matrix and $Q$ is an $ n \times > n $ orthogonal matrix, the first p columns of which form $Q_f$. Recall that > orthogonal matrices rotate vectors, but do not change their length. > Orthogonality also means that $QQ^t = Q^tQ=I_n$. Applying $Q^t$ to > $y-X\beta$ implies that > $$\left\|y-X\beta\|^2=\right\|Q^ty-Q^tX\beta\|^2=\left\|Q^ty-\left[\begin{matrix}R\\\0\end{matrix}\right]\beta\right\|^2 > $$ > > Writing $Q^ty=\left[\begin{matrix}f\\\r\end{matrix}\right]$, where $f$ is > vector of dimension $p$, and hence $r$ is a vector of dimension $n − p$, > yields > > $$ \|y-X\beta\|^2=\left\|\left[\begin{matrix}f\\\r\end{matrix}\right] - > \left[\begin{matrix}R\\\0\end{matrix}\right]\beta > \right\|^2=\|f-R\beta\|^2+\|r\|^2 $$ > > The length of $r$ does not depend on $\beta$ , while $\|f − R\beta\|^2$ can > be reduced to zero by choosing $\beta$ so that $R$ equals $f$ . Hence $$ > \hat\beta=R^{-1}f $$ is the least squares estimator of $\beta$. Notice that > $\|r\|^2=\|y-X\hat\beta\|^2$, the residual sum of squares of the model fit. I find this a bit confusing for a number of reasons. First, it seems to me that $Q$ will only be $n\times n$ if X is $n\times n$. However in general $Q$ will be $n\times p$ where the number of parameters $p$ is less than the number of subjects/cases $n$. It is also stated that the first $p$ columns of $Q$ make up $Q_f$. In order to multiply $Q$ and $R$, $Q$ must contain $p$ columns. Second, if we had a simple design matrix of dimension 2x4 (i.e. we have an intercept and slope paramater and 4 subjects) y would be 4x1 and $Q^t$ 2x4. $Q^ty$ will be 2x1. Thus it is not clear to me how $Q^ty =\left[\begin{matrix}f\\\r\end{matrix}\right]$ where $f$ is vector of dimension $p$, and hence $r$ is a vector of dimension $n − p$. In this case $f = 2$ and $n-p=2$ which would give $Q^ty$ a dimension of 4x1. I have checked the books errata and there are no mistakes listed for this section, so I must not be thinking about this quite right. If anyone could shed some light on what piece of this puzzle I'm missing, I'd greatly appreciate some enlightenment.


title: Binary classification on imbalance data with missing values and almost equal dimensions?

text: I have a data set with following **qualities** :   * Attributes: 500   * Instances: 1500   * Class Ratio: 1:15   * Missing Values: Yes ~5% How should I perform `classification` on this data which can give me best results. The following things I have tried so far so far:   1. Replaced missing values with mean (For time being, will be performing NA imputation using sophisticated methods)   2. `PCA` (dimensionality reduction, reduced no. of columns from 500 to 30 with almost 90% data quality)   3. `Oversampling` the minorities   4. `Ada boost` algorithm for classification Am I correct in my analysis? Although I am getting **100% accuracy** overall and for both classes also. **Is there anyway I can`oversample` this data without doing `PCA`?** I have tried `SMOTE` but that is throwing some colnames() error which I have no clue of. **Is there any package in R for any algorithm like`ada boost`, `random forest`, `ipred` which can take care of cost-sensitivity for such analysis? Or may be a `tuning parameters` for these algorithms by which I can set the cost parameter for minority class?**


title: Hypothesis testing of latent factors - CFA or something else?

text: **Summary:** I constructed a conceptual model with few factors, affecting on 1 DV. All factors and effects are based on theory; however, most of them are vaguely defined and they seem hard to be quantified. In my try to quantify the factors, I am using observed variables or measures. Some factors are measured by only 1 observed var., whereas others are measured by 2 or 3 variables. **Questions** When trying to prove/reject my hypotheses how should I proceed with the factors? E.g. my hypotheses are: * Factor 1 (measured by var.1 and var.2) positively affects DV. * Factor 2 (measured by var.3) positively affects DV. Should I go for Confirmatory Factor Analysis (SEM) and then use the factor scores to test my hypotheses? If yes, is it correct to include in the SEM factors constructed by only 1 observed variables? Thanks for the support


title: Paired-sample t test for usability scores

text: I have conducted an empirical study as a part of my master thesis. Each test subject filled out a questionnaire where some questions were regarding usability of different features. What is the best/correct way to calculate if a feature got a **_significant_** better usability score than another feature? (I used a 10-point Likert scale). I had 14 test subjects, and I want to compare the usability score for a feature called "bubble" and a feature called "hull". This is the scores:               Bubble: 6,  3,  7,      9, 7, 8, 7,  9, 8, 8, 6, 10, 3, 7     Hull:   9, 10, 10, 8.6154, 9, 6, 8, 10, 8, 8, 6, 10, 9, 9      _(One subject did not rate "hull", I have therefore inserted the mean (8.6154) in its place.)_ I have run a paired-sample t test which indicates that there is a significant difference in the usability scores for bubble (M=7.0, SD=2.04) and hull (M=8.62, SD=1.33); t(13)=2.43, p=0.0305. Is this the correct way to do it? My statistical knowledge is quite limited, but I think I should say if there was a _significant_ difference or not.


title: Random effects assumption and testing level 2 predictor variables

text: I was wondering if someone has any advice about the analysis I’m carrying out, or just give a recommended reference? I’m using a random effects modelling approach to account for clustering of patients within clinicians’ practices. I want to test whether the (binary) treatment outcome is associated with a level 2 (clinician-level) variable, i.e., a binary variable representing whether or not the clinician is experienced, with the resulting coefficient estimate indicating the population average difference in the treatment outcome variable between experienced vs. inexperienced clinicians when everything else in the model is held constant. The level 2 variable is associated with unobserved clinician factors, as is the case with common level 2 variables tested in random effects models, e.g., testing private vs. public schools when modelling the exams grades of pupils clustered in schools. I have read that the random effects approach relies on the assumption that the unobserved clinician-level variables are uncorrelated with the level 1 covariates, i.e., the random effects estimator is biased if any of the unobserved characteristics of the clinicians are correlated both with a level 1 covariate and the treatment outcome variable, because this induces correlation between the covariate and the treatment outcome variable. My question concerns how this problem effects level 2 variables. Most of the sources on the random effects assumption that I’ve come across refer to the fact that there must not be correlation with level 1 variables included in the model, but they don’t explain further about the implications for level 2 variables being tested in the model. I’ve come across these two quotes, but I was wondering if someone could expand on them?:   * “A level 2 explanatory variable correlated with level 2 random effects – add second equation with variable as response to avoid bias (+ identifiability constraints)”   * “The random vs. fixed distinction for variables and effects is important in multilevel regression. In multilevel regression models, both level-1 and level-2 predictors are assumed to be fixed. However, level-1 intercepts and slopes are typically assumed to vary randomly across groups.” Thanks!


title: Machine Learning to Predict Class Probabilities

text: I am looking for classifiers that output probabilties that examples belong to one of two classes. I know of logistic regression and naive Bayes, but can you tell me of others that work in a similar way? That is, classifiers that predict not the classes to which examples belong, but the probability that examples fit to a particular class? Bonus points for any thoughts you can share on the advantages and disadvantages of these different classifiers (including logistic regression and naive Bayes). For example, are some better for multi-class classification? Thank you!


title: How can one convert multiple DVs from long to wide?

text: Consider the following example long data.frame with two dependent measures, "Score" and "NewVariable", 1 between subjects variable "Prep" (3 levels), 2 within subjects variables "Day" (3 levels) and "Experiment" (2 levels), and a subject identifier "SID".               example <- structure(list(SID = structure(c(1L, 8L, 12L, 13L, 5L, 6L, 1L,      8L, 12L, 13L, 5L, 6L, 1L, 8L, 12L, 13L, 5L, 6L, 1L, 8L, 12L,      13L, 5L, 6L, 1L, 8L, 12L, 13L, 5L, 6L, 1L, 8L, 12L, 13L, 5L,      6L), .Label = c("S1", "S10", "S11", "S12", "S13", "S14", "S15",      "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9"), class = "factor"),          Prep = structure(c(2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L,          3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 2L,          2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L), .Label = c("Group work only",          "Lecture only", "No instruction"), class = "factor"), Day = structure(c(1L,          1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,          3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L,          3L, 3L, 3L, 3L, 3L), .Label = c("Day1", "Day2", "Day3"), class = "factor"),          Score = c(14, 14, 16, 18, 11, 12, 13, 15, 17, 15, 12, 11,          18, 17, 18, 17, 10, 12, 15, 15, 17, 19, 12, 13, 14, 16, 18,          16, 13, 12, 19, 18, 19, 18, 11, 13), NewVariable = c(-0.887056411864653,          -0.480360621343027, -0.490415963314823, 1.3654758915317,          -1.90913204292831, 0.0300532242614742, -1.84822348735206,          -0.3813757992351, -1.70572162896999, 1.00321046322335, -0.758813794873949,          -0.966033445643038, 0.11876111343571, -1.2312333727132, -0.836123526615442,          0.137868615951057, -1.05143917652043, -0.556162009526374,          0.112943588135347, 0.519639378656973, 0.509584036685177,          2.3654758915317, -0.909132042928306, 1.03005322426147, -0.848223487352064,          0.6186242007649, -0.705721628969986, 2.00321046322335, 0.241186205126052,          0.0339665543569619, 1.11876111343571, -0.231233372713198,          0.163876473384558, 1.13786861595106, -0.0514391765204339,          0.443837990473626), Experiment = c(1, 1, 1, 1, 1, 1, 1, 1,          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,          2, 2, 2, 2, 2, 2, 2, 2, 2)), .Names = c("SID", "Prep", "Day",      "Score", "NewVariable", "Experiment"), row.names = c("1", "2",      "6", "7", "13", "14", "16", "17", "21", "22", "28", "29", "31",      "32", "36", "37", "43", "44", "11", "23", "61", "71", "131",      "141", "161", "171", "211", "221", "281", "291", "311", "321",      "361", "371", "431", "441"), class = "data.frame")      My best attempt so far would involve creating the two corresponding wide datasets, e.g.               library(reshape2)     dcast(example,SID+Prep~Day+Experiment,value.var="Score")     dcast(example,SID+Prep~Day+Experiment,value.var="NewVariable")      ... then manually gluing them back together. Is there a better way?


title: Scaling/Normalisation or Standardization

text: I'm working on SVM and ANN classification tools. In order to improve the classification accuracy, I want to know the best or the recommended data- preprocessing, is it scaling/normalisation or standardization ? and why ? Many thanks


title: Intraclass correlation to validate a measure

text: I would like to understand whether the output of an algorithm I developed agrees with the ratings of external observers. The number of observers is 10 and each of them rated on 2 3-point-Likert items 10 stimuli. The algorithm provides one measure (a real number) for each stimulus. I am thinking about the use of intraclass correlation among observers and then to use this possible agreement (but I don't know if this sounds good and how to do this) to validate the output of the algorithm. Thanks


title: R: Defining a new continuous distribution function to use with Kolmogorov-Smirnov Test

text: I wish to run the Kolmogorov-Smirnov test on my data to determine how well it conforms to a specific continuous distribution function I have in mind. If my understanding is correct, the Kolmogorov-Smirnov test can be performed in R with something along the lines of:               ks.test(my_data, "pweibull", shape=2, scale=1)      where in the case above `pweibull` is indicating a test against the _Weibull_ function. The issue I am having is that the function I want to test against (the BiPareto function defined here for what it's worth) is not a standard part of R best I can tell. Is it possible to define one's own function with parameters in R that can be passed to `ks.test()`? ie. so I can have a call like:               ks.test(my_data, "bipareto", alpha=a_param, beta=b_param, c=c_param, k=k_param)


title: Using LASSO to rank features

text: Suppose I want to answer the following question: Rank the explanatory variables from most important to least important. I want to do the following: Pick lambda from largest to small, and report in a table : when lambda is 100, 80, 60, 40, 20, 10, the explanatory variables are such and such and their loading are such and such. 1.Is it legal to do this? Each lambda represent a Laplace prior, reporting multiple lambda values would be equivalent to have multiple priors. It seems weird we can have multiple priors. 2.If it is not legal, what method should I use to answer my original question.


title: Which non-parametric test for difference between three groups, of which two are paired?

text: I want to compare if three groups are different by a non-parametric test. Now..., the problem is that two groups are paired (prior to treatment and one- year after treatment), and the (reference) group is other healthy individuals. Is there a simple solution? I am using R,


title: Can empirical orthogonal function (EOF) analysis be used as a predictive model?

text: I am wondering if I can use the results of EOF to "predict" the EOF loadings of a second set of data? I have a data set $X$ that consists of sampled light spectra (each column represents a sample and each row is the amplitude of the spectra at a given wavelength). So I first construct a covariance matrix of my data: $C=XX^T$ This is then decomposed using an eigendecomposition: $C=E \lambda E^T$ with $E$ being the EOF loadings. I then project my data $X$ onto the loadings to get the principal components $Z$ : $Z = XE$ Now my question is if I can use the results of the EOF to project a second set of spectra data $Y$ into EOF loadings $E_{\mathrm{pred}}$ using the principal components $Z$ : $ZY^{-1} = E_{\mathrm{pred}}$ I haven't come across EOF being used in such a way (e.g. as a predictive model) but, given the amount of literature around EOF and PCA, I'm probably just missing something. Is there a problem with this approach? Does anyone know of examples using such EOF as a predictive model? Many thanks for any advice.


title: Median Averaging and Error Analysis

text: I am not sure if i should be asking this on the Mathematics of Physics Stack Exchange site. I am taking data for a physics research application and have a question about using median vs. mean averaging. I am concerned with the amount it takes to acquire each data point (location). For each point I currently take multiple measurements, average them, and use the standard deviation of the mean as my error calculation for the point. My problem is that it is only reasonable for me to take a relatively small number of measurement(under 10) for each data point. It seems to me; if any of the measurement vary drastically from the actual value, the median of the measurements will converge to a value close to the actual value in fewer data points than the mean. What do you think of this assumption? If I use the median of my measurements what would be the correct way to calculate the error for each data point?


title: Computing precision of continuous classifier with ordinal outcome

text: I have data where a rater predicts quality on a continuous scale, and quality is then validated in buckets. For example, one rater might give scores of 0.2, 0.4, 0.7, and 0.6 for cases with outcomes of LOW, MEDIUM, HIGH, LOW (respectively). This would be an overall decent rater aside from the fact that the fourth case received a higher rating than the second, yet experienced a worse outcome. I'm looking for a way to compute a precision metric for this type of rating problem. Hand & Till's multi-outcome AUC (R package) measure is helpful context, but insufficient since the raters aren't predicting each individual outcome. Rather, I'm looking to leverage the ordinal nature of the data.


title: Logistic Regression with 1 independant categorical variable, goodness of fit

text: Suppose I fit a logistic regression model on some data with 1 categorical variable, suppose this variable can take 5 different values eg: Coutries, or cities Eg: NY, NJ, London, ..... Now, if I fit this model and convert to odds ratio and finnally to probabilities. Will the probabilities be exactly the same as if I just did proportions or will the difference between the logistic probabilities and actual proportions depend on the goodness of fit ? Many thanks,


title: Data analysis for stratified data

text: I have a sample of 300 companies classified according to different economic activities. For example, computer parts comprises 20% of the sample, medical services 5%, wood processing 15%, etc. In terms of data analysis, what is the best approach to follow? If conducting a general analysis, results will be biased due the differences in proportions among sectors. From this standpoint, which is the recommendation to analyze and present results?


title: How to calculate the probability of absence for a certain category of artefacts from a sample, given prior knowledge about its abundance?

text: In archaeology, artefacts are commonly classified in categories according to certain criteria (those may include manufacturing technique, decoration, function, chronology, etc). I am trying to estimate the probability that objects of a certain category $C$ are **absent** from the site, even after group $G$ has been collected and no category $C$ objects were found. Group $G$ is a “sample” from the group of all objects that were used at the site in a certain moment (actual examples are a waste dump from a house, and domestic tools "frozen" under a collapsed building). Archaeologists will recover all artefacts from the soil but it is an established fact that not all objects will survive, and that we will never recover all of them. The data I have include some prior knowledge about the (relative and absolute) abundance of $C$: the amount and proportion of $C$ in all groups **where it is present** , and the total amount of artefacts in all groups, including ones where $C$ is not found (I have other related data but these variables seem the relevant ones). Groups aren't generally very large, around 300-500. Objects of category $C$ account for 1-3% of artefacts when they are found. I expect that the larger the “sample” the higher the probability of absence will be. A professor of statistics suggested me that methods used by ecologists might be appropriate, e.g. the work of McBride & Johnstone found here: http://www.nzes.org.nz/nzje/contents.php?volume_issue=j35_2 (Credible Interval Value). Is is feasible to apply these methods to the above problem using R, and how? If not, what other methods are appropriate?


title: Cost function for validating Poisson regression models

text: For count data that I have collected, I use Poisson regression to build models. I do this using the `glm` function in R, where I use `family = "poisson"`. To evaluate possible models (I have several predictors) I use the AIC. So far so good. Now I want to perform cross-validation. I already succeeded in doing this using the `cv.glm` function from the `boot` package. From the documentation of `cv.glm` I see that e.g. for binomial data you need to use a specific cost function to get a meaningful prediction error. However, I have no idea yet what cost function is appropriate for `family = poisson`, and an extensive Google search did not yield any specific results. My question is anybody has some light to shed on which cost function is appropriate for `cv.glm` in case of poisson glm's.


title: Present numerical results with variable system parameters

text: I've been working on a problem where I ended up with a potential field expressed in terms of the following integral: $$I(u)=\int_1^u \frac{e^{-x} (2 x-1)}{\sqrt{x~(A~e^{-x}+1)-B \sqrt{x}}}dx \,,$$ where $A,B\in \mathbb{R}$ are constants such that the integrand has no singularities. The integral, unfortunately, can't be solved analytically so I was hopping to show the behavior of the potential in a plot. My problem are the constants $A$ and $B$ which represent some parameters of the system that I'm studying. How can I plot the function $I$? Should I make different plots for various values of the constants? Should I present one plot for various values of $A$ and fixed $B$ and vice-versa? I would be really grateful if someone could share their experience on how to best present such information.


title: Constraints on ML for mixture of Gaussians

text: I have some data sampled from a mixture of two Gaussians where one of them is known, and the density function is as follows: $f(x, \mu, \sigma) = \frac{1}{2}\frac{1}{\sqrt{2\pi}\sigma} exp(-\frac{(x-\mu)^2}{2\sigma^2})+ \frac{1}{2}\frac{1}{\sqrt{2\pi}} exp(-\frac{x^2}{2})$ Now, for cetain values of $\mu$, the likelihood for given data $x_1, ... , x_n$ grows to $\infty$ as $\sigma$ tends to zero. Can we have a lower bound on $\sigma$, $\sigma > \sigma(n)$, s.t. we'll have a good approximation for $\mu, \sigma$? Thanks


title: Select the covariance structure or the model first for GEE?

text: I have data with repeated measurements and a binary outcome (Yes/No). I want to use GEE to model the data with a logit link function. Let $\pi=P(yes)$, then I want to consider the following five models: \begin{align} \text{Model 1:}\quad log(\frac{\pi}{1-\pi}) & = \beta_0 + \beta_1 X_1 + \beta_2X_2 + \beta_3X_3\\\ \text{Model 2:}\quad log(\frac{\pi}{1-\pi}) & = \beta_0 + \beta_1 X_1 + \beta_2X_2 + \beta_3 \log(X_3) \\\ \text{Model 4:}\quad log(\frac{\pi}{1-\pi}) & = \beta_0 + \beta_1 X_1 + \beta_2X_2 + \beta_3 X_4 + \beta_5X_5 \\\ \text{Model 5:}\quad log(\frac{\pi}{1-\pi}) & = \beta_0 + \beta_1 X_1 + \beta_2X_2 + \beta_3 X_4 + \beta_4log(X_5) \\\ \text{Model 5:}\quad log(\frac{\pi}{1-\pi}) & = \beta_0 + \beta_1 X_1 + \beta_2X_2 + \beta_3 X_4 + \beta_4X_5 + \beta_5 X_4X_5 \\\ \end{align} I want to choose among these models using the QIC. Additionally, I want to use the QIC to choose the working correlation structure (independence, unstructured, exchangeable, AR(1)). My question:What should I do first, find the "best" model or find the "best" correlation structure. Does it matter? If I try to find the best among the 5 models first, what correlation structure should I pick? On the other hand, if I try to pick the correlation structure first, based on what model should I do this? My first intuition was to pick the correlation structure first using the following model with all covariates: $$log(\frac{\pi}{1-\pi}) = \beta_0 + \beta_1 X_1 + \beta_2X_2 + \beta_3 X_3 + \beta_4 X_4 + \beta_5X_5$$ However, I wasn't sure if I should include the log- transformations and the interactions, too. On the other hand, since GEEs yield consistent estimates even if the covariance structure is misspecified and the variance sandwich estimator is also robust against misspecification, it seems to be easier to choose among the five models first, e.g. using an exchangeable structure. I tried to look for some literature regarding this, but I couldn't find any, so if anybody knows of any, I would appreciate it, too.


title: Wavelet analysis of EEG

text: I want to do a time-frequency analysis of an EEG signal. I found the GSL wavelet function for computing wavelet coefficients. How can I extract actual frequency bands (e.g. 8 - 12 Hz) from that coefficients? The GSL manual says: > For the forward transform, the elements of the original array are replaced > by the discrete wavelet transform $f_i \rightarrow w_{j,k}$ in a packed > triangular storage layout, where $j$ is the index of the level $j = 0, > \ldots, J-1$ and $K$ is the index of the coefficient within each level, $k = > 0 ... (2^j)-1$. The total number of levels is $J = \log_2(n)$. The output > data has the following form, $(s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0}, > \ldots, d_{j,k}, \ldots, d_{J-1,2^{J-1}-1})$ If I understand that right an output array data[] contains at position 1 (e.g. data[1]) the amplitude of the frequency band 2^0 = 1 Hz, and               data[2] = 2^1 Hz       data[3] = 2^1 Hz       data[4] = 2^2 Hz        until               data[7] = 2^2 Hz       data[8] = 2^3 Hz      and so on ... That means I have only the amplitudes for the frequencies 1 Hz, 2 Hz, 4 Hz, 8 Hz, 16 Hz,...   How can I get for example the amplitude of a frequency component oscillating at 5.3 Hz?   How can I get the amplitude of a whole frequency range, e.g. the amplitude of 8 - 13 Hz?


title: What exactly does it mean to 'pool data'?

text: I thought that 'pooling data' simply meant combining data that was previously split into categories...essentially, ignoring the categories and making the data set one giant 'pool' of data. I guess this is a question more about terminology than application of statistics. For example: I want to compare 2 sites, and within each site I have two year- types (good and poor). If I want to compare the 2 sites 'overall' (that is, ignoring the year-types), is it correct to say that I'm pooling the data within each site? Further to that, since several years of data comprise the good and poor year-types, is it also correct to say that I am pooling the data among years to achieve the 'good year' and 'poor year' data set within each site? Thanks for your help! Mog


title: Computing Gaussian mixture model probabilities

text: I was looking over the solution to this question on SO and it got me thinking about computing probabilities for a Gaussian mixture model. Let's assume you've fit some Gaussian mixture model so that it results in a mixture of three normals: \begin{equation} X_1 \sim \mathcal{N}(\mu_1,\sigma_{1}^{2}), \quad X_2 \sim \mathcal{N}(\mu_2,\sigma_{2}^{2}), \quad X_3 \sim \mathcal{N}(\mu_3,\sigma_3^2) \end{equation} with respective weights $\lambda_1, \lambda_2$, and $\lambda_3$. From here on out, take $\mathbf{X}=[X_1,X_2,X_3]$ and $\mathbf{\lambda}=[\lambda_1,\lambda_2,\lambda_3]$. Typically, to find the probability that this model is less than some value $x$, we find \begin{equation} \mathbf{P}(\mathbf{\lambda}\mathbf{X}^T\leq x) = \sum_{i=1}^{3} \lambda_i \mathbf{P}(X_i\leq x) \end{equation} It's possible that I've made a coding error, but it seems that the probability obtained from the above formula is different from the probability obtained if we compute the probability in a different way: \begin{equation} \mathbf{P}(\mathbf{\lambda}\mathbf{X}^T\leq x)=\mathbf{P}(Y\leq x) \end{equation} where $Y\sim\mathcal{N}(\mathbf{\lambda}\mathbf{\mu}^T,\sum_{i=1}^{3} \lambda_{i}^{2}\sigma_{i}^{2})$ with $\mathbf{\mu}=[\mu_1,\mu_2,\mu_3]$. If it's a coding error, please just leave a comment and I will delete this question.


title: Calculate F-statistic / p-value for subset of co-efficients in R

text: I'm wondering if there's an easy way of calculating an F-statistic / p-value for a subset of model coefficients. Particularly in R? I'm not sure what test would be needed to calculate this. For example,               summary(lm(a~w+x+y+z))      will give me the F-statistic and p-value for the whole model, but is it possible to extract the F-statistic and p-value for just the terms y and z? A faster way of doign something like this:               matrixOfResponses <- cbind(c(1,2,3,4,5), c(4,3,2,4,5), c(5,3,2,23,4), c(1,2,4,3,1,))     pValsOut <- numeric()     for(i in 1:ncol(matrixOfResponses))     {     pValsOut[i] <- anova(lm(matrixOfResponses[,i]~mMat), lm(matrixOfResponses[,i]~mMatReduced))$'Pr(>F)'[2]     }      So basically something like the above becomes very slow when "matrixOfRespones" contains a huge number of variables, for example a huge number of gene expression levels.


title: Calculation of confidence interval of a population parameter

text: Consider that $P$ is the _water pressure_ coming out from a valve A. Let $P_{dif}$ be defined as the difference between the maximum and the minimum pressure of valve A: $$P_{dif}:= P_{max} - P_{min}$$ Now, what I want to do is to estimate $P_{dif}$. In order to do that, I take a number of water pressure samples from valve A. Let $S$ be a set of 3 measured samples: $$S = \\{X_1 = 5, X_2 = 7, X_3 = 1\\}$$ That is, $S$ contains 3 random samples drawn from the population. Therefore, I then say that $\hat{P}_{dif} = 7 – 1 = 6$. **Fisrt question:** Considering a Gaussian distribution of the population parameter, how can I find the pdf, cdf of the parameter? **Second question:** Would a $95\%$ confidence interval (CI) of my estimation be defined as following?: $$(pdf, mean) - 2\cdot (pdf, std) \leq (pdf, mean) \leq (pdf, mean) + 2\cdot (pdf, std) $$ where (pdf, mean) is the mean of the pdf and (pdf, std) its standard deviation. If yes, how can I derive the (pdf, mean) and the (pdf, std)? **Edit:** As it seems quite difficult for me to find the CI of the range estimation, I am wondering whether I could use an alternative approach. _whuber_ suggested that successive estimate should not depend strongly on the preceding ones. Say that I split my data samples into equal parts (e.g. 3) and I find one estimate for each of these sub-sets. Next, I find the 95% CI of the $\hat{P}_{dif}$ by that: $$mean - 2\cdot std \leq mean \leq mean + 2\cdot std$$ where _mean_ and _std_ are the mean and std of the 3 different estimates, respectively. Would that be correct?


title: On connection weights in an Artificial Neural Network

text: I play a lot with PyBrain \-- Artificial Neural Network implementation in Python. I have noticed that in all the models that I receive the weights of the connections are roughly normally distributed around zero with a pretty low standard deviation (~0.3), which means that they are effectively limited within the [-1, 1] range. What does this mean? Is it a requirement of ANN? An outcome of backpropagaion learning? A sign for network health? Or just a random observation? ![enter image description here](http://i.stack.imgur.com/5q8YS.png)


title: From standard HMM to Bayesian HMM

text: I'm trying to understand what the difference between a standard HMM and a Bayesian HMM is. Wikipedia just briefly mentions how the model looks like but I need a more detailed tutorial. Does someone know of a paper or an implementation which I can look at? Also I have issues with the terminology used. What does it practically mean if you "place/put a Dirichlet prior on a distribution"?


title: How to calculate the difference of two standard deviations?

text: ## Background I have two estimates of variance and their associated standard errors calculated from sample sizes of $n=500$ and $n=10,000$ the results are $\hat{\sigma^2} (sd_{\hat{\sigma^2}})$: $$\hat{\sigma^2}_{n=500}=69 (6.4)$$ $$\hat{\sigma^2}_{n=10,000}=72 (1.5)$$ ## Question If I say that variance increased by 3, what is the standard deviation around this estimate? ## Notes   1. $SD$ of var calculated using $SD_{\hat{\sigma^2}}=\sqrt{s^4(2/(n-1) + k/n)}$   2. I suspect that the fact that I am estimating the sd of a variance is not relevant to the calculation, but may help in the interpretation of what I am doing.


title: Data normalization and standardization in neural networks

text: I am trying to predict the outcome of a complex system using neural networks (ANN's). The outcome (dependent) values range between 0 and 10,000. The different input variables have different ranges. All the variables have roughly normal distributions. I consider different options to scale the data before training. One option is to scale the input (independent) and output (dependent) variables to [0, 1] by computing cumulative distribution function using the mean and standard deviation values of each variable, independently. The problem with this method is that if I use the sigmoid activation function at the output, I will very likely miss extreme data, especially those not seen in the training set Another option is to use a z-score. In that case I don't have the extreme data problem; however, I'm limited to a linear activation function at the output. What are other accepted normalization techniques that are in use with ANN's? I tried to look for reviews on this topic, but failed to find anything useful.


title: How to equalize lower accuracy on larger sample sets?

text: I have a collection of image samples that I want to OCR. I also have the ground truth of the text on these images. Now I OCR them and eg. 10 of 100 samples are wrong. So I got an accuracy of 90%. Naturally the accuracy decreases, when I increase the size of the sample set. So if I increase it to 100 samples I got an accuracy of only 70%. Another set of images has maybe only 2 images samples. Even with very bad and specific OCR settings the accuracy could be 100% (pretty much by accident), while increasing the sample size of this this could bring the accuracy down to 2% (if these 2 samples are the only correct ones). I would like to make these values comparable, so that a higher accuracy on a small sample set, results in an equally low value as a low accuracy on a larger sample set, because it is a more "useless accuracy".   1. **So a measured accuracy has a higher value when it is measured on a larger sample set?**   2. **Do I have a more accurate accuracy, when the sample set is larger?**   3. **Can I equalize this accuracy of accuracy, to a common value?** I tried to find the answser here: http://en.wikipedia.org/wiki/Accuracy_and_precision but, I can't apply this because I have no binary classification (it is more unary). I don't have false positives, to calculate a precision. I hope someone can point me in the right direction.


title: L1-norm cost function for Neural Network. (Regression)

text: I am trying to build a regression model using a neural network. The final cost measure is the mean absolute error (MAE) on the output (one output unit, 200 input units). Right now all my hidden units have rectifier activation. The output unit is just a linear unit with pass-through activation. It seems the network cannot learn efficiently, the error (even on can't find a value that makes the error go down monotonically). I suspect the cost function (L1-norm) might be the culprit. Right now, when taking the gradient, I either pass 1 or -1 depending on predicted value vs actual output value. Is this the right way? (Since L1 is not smooth at 0, would this be the reason why the learning is not smooth/effective?) What is the right way to handle a L1-norm cost function? Thanks, any help is appreciated!


title: Tests of within subjects contrasts in R

text: I have run a split splot model in SPSS (via repeated measures function) and I would like to reproduce my results using R. To do so I used ezANOVA function from ez package to obtain sphericity tests and correction and type III SS. I have read that aov function does not give type III SS. So far the results are identical with those from SPSS but I cannot tell so when I run tests of within subjects contrasts. I use polynomial contrasts and my factor has six levels. What I manage to get is t-values instead of f-values or f-values but not the same as in SPSS. Here is my data and code               id    subj  treat  m1    m2    m3    m4    m5   m6     1       1     1   455   460   510   504   436   466     2       2     1   467   565   610   596   542   587     3       3     1   445   530   580   597   582   619     4       4     1   485   542   594   583   611   612     5       5     1   480   500   550   528   562   576     6       6     2   514   560   565   524   552   597     7       7     2   440   480   536   484   567   569     8       8     2   495   570   569   585   576   677     9       9     2   520   590   610   637   671   702     10     10     2   503   555   591   605   649   675     11     11     3   496   560   622   622   632   670     12     12     3   498   540   589   557   568   609     13     13     3   478   510   568   555   576   605     14     14     3   545   565   580   601   633   649     15     15     3   472   498   540   524   532   583      m1-m6 represent measurements in 6 different time points. My code for the contrasts.               long.df<- melt(data, id=c('subj','treat'))     long.df<- long.df[order(long.df$subj)]     names(long.df)<- c('subj','treat','time','meas')     mod<- lm(meas~time + time^2 + time^3 + time^4 + time^5, long.df)     summary(mod) # this is how I obtain t-values for the polynomial contrasts #     mod1<- aov(meas~time + time^2 + time^3 + time^4 + time^5, long.df)     summary(mod1, split=list(time=list("Linear"=1, "Quad"=2,'q'=3,'f'=4,'fif'=5)))      # this is how i obtain f-values #      I think that the F-values correspond to type 1 SS (as I use aov function). How am I wrong? Is there a way to conduct the Tests of within subjects contrasts and obtain type III SS so I will have identical results to SPSS? Here are the SPSS results. Many thanks in advance!               Source     time    Type III SS/   df/    Mean Square/        F/   Sig.       time     Linear     123662.881/    1/     123662.881/   83.591/   .000             Quadratic       5928.007/    1/       5928.007/   18.968/   .001                 Cubic      10462.676/    1/      10462.676/   28.075/   .000               Order 4        798.193/    1/        798.193/    4.010/   .068               Order 5       1702.743/    1/       1702.743/    4.878/   .047


title: Modeling Difference-in-Difference with Multiple Treatments

text: I ran a field experiment across 100 cities. The cities belong to two types: Type_A(50 cities) and Type_B(the other 50 cities). From January 1 to February 30(The Pre-Treatment Period), I only collected data on my variable of interest y. From March 1 to April 30 (The Post-Treatment Period), I simultaneously applied two treatments on two different sets of cities: I randomly chose 10 Type_A cities and 10 Type_B cities for Treatment_1, and another 10 Type_A cities and another 10 Type_B cities for Treatment_2. I plan to use DiD (Difference-in-Difference) for modeling. How does the specification look like? I am confused because there are several things going on together: different types of cities, multiple treatments, before and after periods, etc. In addition, if I want to implement the specification in Stata, how to do that? More importantly, how to compute the marginal effect of the treatments in Stata? Thanks a bunch in advance.


title: Panel Cointegration, Moderating Effects and Multicollinearity

text: I am running a panel fixed effects regression on 21 countries and 16 years. Its a secondary data taken from OECD website mostly. My model looks like this: $$ \log{(GDP/Labor)}_{i} = \beta_0 +\beta_1 \log{(RnD)}_{i,t-1}+\beta_2 \log{(Human Capital)}_{i,t-1}+\beta_3 \log{(x_j*(RnD_j)}_{t-1}+\epsilon_t $$ Where,   Rndi = Research and Development expenditure of country i   Rndj = Research and Development expenditure of country j   HC = Index of Human Capital 1). I have following problems with this model: 1) all variables have strong trend and they are expected to be non-stationary however, Im-Pesaran-Shin and Levin-Lin-Chu tests are giving contradictory results about non-stationarity so I cannot conclude about non-stationarity. In your opinion, which test is better for panel regressions? Also, is there any test I can do for cointegration in stata for panel regressions in addition to Westlund test? if I dont find evidence for cointeration, do I first difference the data? I lose ALOT of information and everything becomes insignificant if I do that. 2) One of the primary goals of my paper is to test for moderation effects of human capital on knowledge spillovers and productivity nexus. Problem is that there is a serious multicollinearity problem. VIF goes as high as 1000! How can I test for moderating effects other than using interactions? Moreover, in some cases I also find strong multicollinearity between my key variables and it really makes everything significant at < 0.001 with R2 less than 50. 3) How important is it to include year dummies in regressions? I lose all significance of my variables and things become very bad in terms of signs as well for my core variables when I use year dummies. I will be extremely grateful if you could please advice. Thanks alot Ali


title: Wald test for gamma distributed data?

text: One of the assumptions of Wald test is that the data is normally distributed. If my data follows Gamma distribution. And I have estimated the two parameters of the gamma distributions. Is there a test like Wald test I can use?


title: Regression of a dataset using median

text: I have a dataset containing ( _x_ , _y_ ) information from the years 1985–2013. The _x_ data is the pick number and the _y_ data is performance. Each year has 60 ( _x_ , _y_ ) data points. When looking at my data for individual years, I personally think the median of performance of each pick number is a better representation of performance than the mean. I want to run a regression to look at the relationship between pick number and performance. Is it 'wrong' statistically to first find the median of each pick's performance data and then run a regression in order to find a relationship between pick and performance? Any help would be greatly appreciated.


title: Should a non-significant adjustment variable be kept in a regression model?

text: I'm working with a structural equation model to study influenza infection risk. As age is a known risk factor to explain infection, I therefore adjusted my infection outcome on the subjects age class. After regressing my outcome on several latent variables, the age class is not significant anymore. My question is : should a non-significant adjustment variable be kept in a (structural) regression model or removed ?


title: Removing insignificant variables?

text: Suppose you fit a linear regression model on some data with 10 variables. The _F_ -statistic shows that 3 of them are significant ( _p_ < 0.05) , 2 are within trend (0.05 < _p_ < 0.10), and the other 5 are _p_ > 0.10. Is it OK to throw out the other 5 ( _p_ > 0.10) and retrain another model? Perhaps this time the 2 within trend will become significant.


title: General Linear Model with insignificant predictors

text: I have run a GLM- I have 1 DV and 9 predictors (4 of which being demographic variables such as age, gender etc). Results show that only 5 predictors are significant predictors with an R square of 34.1%. When I re-run the model after removing the insignificant predictors, I get an R square of 31.5% which is less than before. I understand that this means that I have lost some of the explanatory power of the model. Is it the right decision to remove the insignificant variables? Is a loss of (34.1-31.5%) 2.6% in R square low enough to accept? The purpose of the GLM is to explain an observed behaviour, which is the dependent variable. Something apart from this, I have been trying to assess the homogeneity of variables of this model. It was not so clear to me from a scatterplot of residual vs predicted values. However, SPSS is able to generate a Levene test with the GLM. These are the results: F= .684 df1= 7 df2= 141 Sig= .686 From my understanding, this means that my predictors are not heteroscodastic and that the assumption of homoscadasticity is met. Can someone please confirm this?


title: Is it necessary to do $k$-fold cross validation for decision trees in random forests?

text: Consider the following data set **train** :                   z   a   b   c         0   1   40  185         0   1   25  128         0   0   32  100         0   0   29  100         1   1   30  107         0   0   30  133         1   1   38  132         1   1   37  127         1   0   30  184         1   0   40  199         1   1   26  185         0   1   21  185         0   0   21  134         0   0   20  137         1   1   22  135         0   0   23  189         1   0   32  109         1   0   31  152         1   0   38  130         1   1   37  191         0   1   39  168         1   0   28  183         0   1   26  171         1   1   23  164         0   1   32  111         0   0   34  131         1   0   30  121         1   0   27  195         1   1   29  117         1   0   26  187         1   0   34  183         0   0   28  189         0   1   34  150         0   1   34  176         0   1   24  140         1   0   37  181         0   1   36  109         1   0   39  198         0   0   32  164      where **z** is a binary variable with predictors **a,b,c**. Suppose that there is some other **test** set with the same variables as the **train** data set and we want to predict **z**. For decision trees, is it better to use the full **train** data set to construct the tree? What would the purpose of $4$-fold cross validation be for example? In a random forest, is $k$-fold cross-validation necessary? I thought you could use OOB error?


title: Intuition related to standard deviation as threshold

text: I have a set of input output training data, few samples are               Input  output     [1 0 0 0 0]  [1 0 1 0 0]     [1 1 0 0 1]  [1 1 0 0 0]     [1 0 1 1 0]  [1 1 0 1 0]      and so on. I need to apply standard deviation on the entire output as a threshold. So, I calculate the mean standard deviation for the output. The application is that the model when presented this data should be able to learn and predict the output. There is a condition in my objective function design which is the distance = sum of the sqrt of the euclidean distance between model output and the desired target, corresponding to an input should be less than a threshold. My question is how should I justify the use of threshold? Is it justified ? I read this article article which says that it is common to take standard deviation as the threshold. For my case, what does it mean taking the standard deviation of the output of the training data?


title: How to compare a measure reflected by a regression residual in a pre-post design?

text: Maybe a simple question, but my head is spinning: I want to see if X increased after a certain treatment in a pre-post study. Unfortunately, I have only K which is a noisy measure of X. To "clean" K from potential influences unrelated to X, it is common in the respective literature to estimate a (first-stage) regression and then use the residual as a better measure K* of X, e.g., K = a + b*Influence + e; e = K*. Now my question:   1. In the literature the regression K = a + b*Influence + e is usually estimated by year... so how can I compare the resulting measure K* (= e, which should be by construction mean 0) through time then? Is it feasable to just suppress the intercept in the first-stage regression? How can I interpret the difference between the two time-periods with a suppressed intercept?   2. If i want to take the absolute value of e, |K*|, for each year (another measure of an alternative construct, say X2) and compare between the two years, how does the inclusion/exclusion of the intercept in the first-stage regression affect my comparison? Thanks!


title: e1071 svm queries regarding plot and tune

text: I am new to R and I am learning the e1071 packages' svm function. Following are the few questions I have.   1. How does the plot function work? I cannot understand the plotting case with more than 2 factors in the class variable. Like in the line below, why `Petal.Width` and `Petal.Length` are chosen as the two dimensions for plotting and how other dimensions affect the result. Moreover I am also not clear with the slice parameter, what is it and why the values 3 and 4 are set. `plot(model,iris, Petal.Width ~ Petal.Length, slice= list(Sepal.Width=3,Sepal.Length=4))` 2 What is tune function and how we choose the gamma and cost variables as shown in this example in cran documentation.               obj <- tune.svm(Species~., data = iris, gamma = 2^(-1:1), cost = 2^(2:4))


title: Two class unpaired test using LIMMA

text: I want identify the set of genes that are significantly-differentially expressed in two different classes of samples. For Example, Normal and Disease. Both classes have samples from different stages, some stages are common and some not. So, I was thinking of performing a 'two class unpaired test' similar to the SAM (Significance analysis of Microarrays). Is there any similar test in Limma package (Bioconductor)?


title: How would I correlate many variables to a few coefficients?

text: I'm trying to predict the "strength" vs "tempertature" and "time" curves of some chemical compounds as a function of the concentrations of their component substances. I have 20 different substances. From these, I created 550 chemical compounds by combining up to 7 of the substances. Among the different compounds, the concentrations of the various substances varies.. For each compound, I experimentally measured strength as a function of temperature and time. I have selected a particular family of functions to fit the experimental data. Each member of this family is specified by specifying the values of 4 constants. I would like to predict the value of the 4 constants as a function of compound's composition - i.e. the substances that are in it and their concentrations. If I had one substance at varying concentration, I know how I'd go about finding it's correlation to the coefficients. I'd find the trend of: $coefficient_i(substance,concentration)$ But I have no idea how I'd go beyond just one substance, and I've got a lot more than two! I'm hoping someone can point me in the right direction.


title: What does the tests of "scale" check?

text: I have found this test: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/mood.test.html unfortunately I didn't find any information on wikipedia, but what "scale" means? What do those tests check? Thanks


title: How to test which combination of product attributes generates better sales performance?

text: I am looking for help in identifying what is the best statistical approach to address the following market research question: (1) I have a dataset from an online website that helps users rent their vacation homes. For each vacation home I have information on:   * Their attributes (predictor variables): location (city and distance from city center), # of rooms, guest rating (10 stars),price, etc.   * Their sales performance (dependent variables): conversion rate, # of nights sold, sales volume in $ (2) The website rents vacation homes mainly for US and Canada customers **Question:** I would like to be able to identify which combinations of predictor variables are associated with better performance, in order to be more targeted in gathering new homes for the site. Because most of the website users (renters) originate from the US and Canada (and I have a sense that they have different preferences) I would be especially interested in understanding which combinations of attributes work best for each country? I imagine that for Canadian customers you could have 3 combinations of attributes that would work best. For instance, one combination could be: price $50-$75/day, 2 rooms, far from the city center, etc... Any help would be really appreciated.


title: Finding top n across multiple sets

text: Suppose we have a set of key value pairs, such that the keys can be ordered based on the value . We can easily extract and report the top $n$ keys. Now suppose the we havet S such sets. From each set we are supposed to extract a set of $p$ top keys and merge them to get the top $n$ keys from all the sets. How large should $p$ be? What is the relationship between the $S$, $n$ and $p$ and the error? I wanted to know if there is any literature on this subject. Clarification: Its possible to have the same keys in multiple sets. The merge and find top n in the ideal case would construct a combined map of key-value pairs (by adding all values corresponding to the same key) and finding the top n values. In my data the distribution of values in each set and the combined set has a long tail.


title: Results of Wilcoxon test are backwards when compared to R, don't understand why

text: I have my own implementation of a Wilcoxon signed rank test. In my code, I of course calculate W- (sum of ranks of negative numbers) and W+ (sum of ranks of positive numbers); the function provides a parameter for the user to select the alternate hypothesis, just as I have seen done in R:   * lesser - the true location is less than zero   * two-sided - the true location is not zero   * greater - the true location is greater than zero I also calculate the distributions of ranks as follows: if we have n ranked (non-zero) values, let a[i] (i from 0 to n(n-1)/2) be the outcomes for n-1. Then the outcomes for n are b[i+n] = a[i] + a[i+n] (with a[j] = 0 for j out of bounds) Now for lesser I calculate the pvalue as sum of b[i=0..W-], for greater b[i=0..W+], and two sided as 2 b[i=0..min(W-,W+)]. But it turns out that my results are opposite of R -- i.e. my 'greater' p-value is R's 'lesser' p-value, and vice versa. Is there some aspect I am not understanding with how to calculate the p-value?


title: How do I weigh several variables at two levels to create an overall composite variable score?

text: ### Context I am trying to find the **correlation between two latent variables**. Let's call them A and B.   * A has two dimensions. Let's call them Dimension 1 and 2.    * B also has two dimensions. Let's call them Dimension 3 and 4.    * The four dimensions (1, 2, 3 and 4) have several indicators each. Let's call them (a) and (b) for Dimension 1, (c), (d) and (e) for Dimension 2, (f) and (g) for Dimension 3 and (h), (i) and (j) for Dimension 4.   * I have measured the ten indicators (a, b, c, d, e, f, g, h, i and j) on a likert scale or else using cardinal numbers (e.g. one, two, three, four, etc.).   * I have then weighted the indicators using the likert scale codes that I created (e.g. 0 = not important, 1 = low importance, 2 = medium importance, 3 = high importance and 4 = absolute importance)   * I have then added the scores for all indicators to arrive at the dimension scores. This is my first weighting.   * I have then said that each dimension contributes 50% to the composite score of A or B and have standardised the dimension scores to 50% each. This is my second weighting.    * Finally I added the dimension scores to arrive at the composite score.  **To give an example:** If I have a score of 2 for (a) and 3 for (b) then my total for Dimension 1 is 5. Similarly, if I have a score of 3 for (c), 4 for (d) and 2 for (e), then my total for Dimension 2 is 9. The score of 5 and 9 are the dimension scores, based on likert scale weights. As I said that each dimension contributes 50%, I have divided the obtained score by the total score for each dimension and multiplied it by 50%. So I have got the following: (2+3)/8 (remember this is based on two likert scale responses so the maximum score is 4 plus 4) x 50% = 0.31 (3+4+2)/12 (remember this is based on three likert scale responses so the maximum score is 4 plus 4 plus 4) x 50% = 0.38 I then added 0.31 and 0.38 and got 0.69. This is my composite score for Composite Variable A (which has 2 dimensions and 5 indicators; 2 for the first dimension and 3 for the second dimension). To make it comprehensible, I have multipled 0.69 by 100 to get 69 out of 100 scores. I want to run corrleation on the composite scores of A and B to see the relationship between them. As you would see, I have **not done any complex statistical modeling** using factor analysis etc but arrived at my final score using a system of subjective weighting based on relative importance (which is aligned to the likert scale coding in the first instance) and equal weighting in the second instance. I have used this approach to keep everything simple. I know it is subjecive weighting but I can justify it to some extent with theory and the rest with commonsense! ### Questions   * Am I doing this right? Is this approach robust enough? **Remember, I don't want to complicate things!** I am not statistically minded!   * Can I do correlation analysis on the composite scores of A and B? They are rank scores (based on weights) so is the Spearman Rank Correlation the way to go for the analysis?


title: How do I create a predictor for a time series once I've confirmed Granger-causality?

text: I have a set of time series data that I've found granger-causality (i.e. regressed Y vs. X, X-1, Y-1), and am wondering how I can create a predictor from this linear model? Is it simply the coefficients I get back from my model? (I'm trying to code this programmatically, and am still a stats rookie, so if this is obvious, forgive me)


title: Inference to the population when the survey response rate is only 30%

text: I have conducted a survey in which the questionnaires were sent out to 450 individuals, but only 30% of them answered the questionnaires.   * Is it still valid to interpret the usual inference analysis (i.e., the inference analysis developed under the assumption of random sampling)?    * Is it correct to do logistic regression analysis with these data?    * If not, how to proceed?


title: Role of kalman filter prior, the 'right' prior?

text: We can solve the Riccati equation and get the steady state priors for a Kalman filter. So why are people still allowed to have any prior they want? They can solve Riccati equation _ex-ante_ and use the steady state prior as the initial value, and reach steady state at time 0 rather than time infinity. And this is the only 'right' prior any rational individual would use. I understand this is a philosophical question, any comment is welcome.


title: Different results by using chi square test and logistic regression

text: I am working on a problem to identify the risk factors to infection after operation. So there are risk factors such as age, pre existing condition, cause of infection etc. Since the dependent variable and most of variables are categorical data, I used logistic regression first, in which the factor A is not significant. Just want to double confirm, I used chi square test just between the dependent variable and A, this time the P value is 0.03 which means they are correlated. Can anyone give me a hint? Does this imply the correlation between the factors?


title: Clustering with the same number of objects per group

text: How can I cluster while forcing the final clusters to have more or less the same number of objects? I have just tried `kmeans` in R, and it does not take this into account.


title: Descriptive measures for 2x2 crosstabs with zeros

text: For 2x2 contingency tables I already know the RR (Relative Risk) and the OR (Odds Ratio). However none of these two is good dealing with empty cells in the table.                  |  A  | -A  |          B | n11 | n12 | n1.     -B | n21 | n22 | n2.        | n.1 | n.2 |  N       OR does not like n12 or n21 to be zero. RR won't handle n11+n21=0 and n12=0 if I am correct. I am looking for some indicator, that does not take infinite values (e.g., when dividing by zero) like the very simple percentage. My idea is to take the relative difference from the expected value for n11. Relative Difference: RD = ((n.1 * n1.) / N - n11) / (n.1) (If n11+n21 = 0 => RD=0) So the question is, if this is already used somewhere so I can quote something or if I have to "proof" that this is a useful indicator. Or can you come up with some other ideas? I will use this indicator for some visual representation, so it is quite useful if it has an upper and lower bound.


title: Is WEKA's logistic regression tool slow? If yes, are there faster alternatives?

text: I have to run logistic regression several thousand times, within a Java/Scala program. A general spec for one of the problems is: ~300 continuous attributes ~<10000 positives, and randomly pick double the number of negatives as positives The attributes are the same for all problems, and the positives/negatives are always drawn from the same pool, of about 100,000. I'm finding that naively plugging this into WEKA takes a few minutes per problem, which is unsuitable given the number of tasks which need to be performed. Is there a faster library I can use, or is there a way to parallelize these problems, or something else?


title: Can multinomial regression do what I want?

text: I have a list of observations of females and their responses, categorized into 5 behaviors, to a potential threat. I'm wondering if, for each response type, whether the presence of an infant makes it more or less likely a female would perform that response. (E.g., we might hypothesize that females with infants are more likely to hide and less likely to go looking for food.) The following was suggested:   1. Bin females by social group (i.e., each bin would only contain females who associated with each other).   2. Calculate the proportion of each response type observed in each group (all these, for each group, would add to 1).   3. Arcsine-transform the proportion data.   4. Perform a t-test, for each behavior, to see whether the proportions of females responding that way differs depending on the presence of an infant. This sounded sketchy to me, and it seems that others agree. The consensus seems to be that, in such cases, it's better to use multinomial regression (in this case, BEHAVIOR ~ INFANT_STATUS) to determine whether infant presence has an effect. However, I was wondering whether I could use those results determine whether (and in what direction) the presence of an infant affects the probability of each response behavior. Also, would that be possible if I were to include additional independent variables? Any advice you can give, on analytical design or actual coding, is much appreciated. I've been doing this in R, but I'm more comfortable in MATLAB if that works as well. Thanks in advance.


title: Kmeans algorithm cyclical solution

text: I am currently implementing a Kmeans clustering algorithm in R. I am not using any packages and I wrote it from scratch. I am using only one set of initial guesses, and my action upon finding an empty cluster is to select a new data point randomly and use that as the new mean for the empty cluster. I have gathered from reading online that the solution does not always converge, and it is highly sensitive to the initial means, so when I see that behavior I am not surprised. But I am finding that sometimes my solution is actually cycling between two or more different solutions. So I have two questions associated with this observation: 1) Within a solution cycle, one solution is always better than the others as measured by the total sum of squared distances of all points to their nearest clusters. So this implies that not only does the algorithm not necessarily find the global optimum, but also it sometimes does not even improve the total sum of squared distances from one iteration to the next? I thought the solution was at least always improving... 2) What is the best way to get around this problem? Do I have to program it to recognize cycles and then select the iteration in the cycle with the lowest total distance? Or is there an easier way? Any help would be greatly appreciated.   Thanks.