title: SVM prediction sensitivity when compared to neural networks and logistic regression

text: Basically I want to classify a rather rare status (about 2% of the 2000) with some predictors. I have used logistic regression, neural network, and Support Vector Machines to do it. All the predictors in the logistic regression are statistically significant. And to avoid overfitting, I have self-implemented a 10-fold CV for all of the methods. For each iteration, I used the training dataset to fit the model and find out the fitted values. Than I use the `ROCR` package in R to find out the decision criterion to achieve 70% sensitivity for the training dataset. Then I use the model and the criterion to predict the status of the testing dataset and compute the testing sensitivity and positive predictive value. After 10 iterations I got 10 testing sensitivities and PPVs. My finding is: the logistic regression did the best. The testing sensitivities were roughly 70% and PPVs were around 16%. But very surprisingly, the performance of SVM is very poor: mean testing sensitivity = 43%, PPV = 11%. I am not very familiar with the theory behind SVM so I tried both `kernlab` and `e1071` in R. I have also experimented with `C-svc`, `nu-svc`, `C-bsvc`, as well as tuning the svm using `tune.svm` in `e1071` but the performance was similar. So my question is: was I doing something wrong, or was I missing something when fitting a svm?


title: Poisson regression with (auto-correlated) time series

text: I have a time series dataset which shows, for each day, the number of complaints received by an organization about a particular problem. I also have a number of other time series for the same period (mostly environmental variables like weather, chemistry etc.) which may help to explain the pattern of complaints. My response variable is therefore discrete (number of complaints), while my possible explanatory variables are all continuous. Based on the advice of a colleague and a bit of Googling, it seems that some kind of Poisson regression model should be appropriate here. However, I'm running into a few difficulties:   1. The variance of my count data (0.76) is greater than the mean (0.33). I think this indicates that my response variable is "over-dispersed" and that Poisson regression is therefore inappropriate?   2. My count data is overwhelmingly dominated by zeros (i.e. most days have no complaints). It seems that this may be a problem for Poisson regression?   3. Virtually all of my time series (both the response and explanatory variables) are auto-correlated.   4. Some of my explanatory variables are correlated with each other. For example, stream flow is closely related to rainfall etc.  If I've understood what I've been reading correctly, all of this suggests that I'll need to do something a bit more involved than a "standard" Poisson regression. **My question:** What is the **_simplest_** way to meaningfully analyze this kind of data, please? What are my options and what techniques should I be researching? I've been reading a little bit about Box-Jenkins models for count data (?), but I'm getting out of my depth. Is there anything more straightforward that's still rigorous? The ultimate aim is to come up with some suggestions as to which factors _might_ be causing the problem. If I could end up saying something like, " _it appears that explanatory variables x, y and z are significantly correlated with the number of complaints_ " or " _no combination of the explanatory variables can adequately explain the pattern of complaints_ ", that would be useful. Thanks very much!


title: Using days of the week in regression

text: I am trying to recreate what my predecessor has done, with the only knowledge being he used "regression." I'm predicting my dependent variable based on the day of the week, but also on two other quantitative variables. I initially intended to use indicator variables for the week day, but then remembered that I will not be able to get a coefficient for each of the seven days, though my predecessor was able to somehow get a model with coefficients for all of Sunday-Saturday. Is there another way I can tackle this problem, or is this method the best way? Is there a method I can use to get separate coefficients for each day of the week?


title: Comparing two dependent variables with different regression models

text: My two separate response (dependent) variables are counts of two mesopredator grouper species on coral reefs. I have tried various models (LM, GLM Poisson and Gaussian and zinf and zinb) both with untransformed response variables and with square root transformed response variables. The best fit model for one species (due to the distribution of the data) is a Gaussian GLM with square root transformed counts and for the other species is a zero-inflated Poisson model with untransformed counts.   * Is it appropriate to use different models if I want to compare the effects of the explanatory variables on each of the species?   * If not, which model can I use to fit both species count data?


title: Logistic Regression with dependent observations

text: I have a dataset that contains 100 different patients over 5 year’s period. Every patient is examined each month with regard to particular illness and marked as healthy or ill (0 or 1). Every person appears 60 times in my sample (5 * 12 = 60). Every month patient provides A = Average blood pressure in that month, B = Average daily exercise hours and C = Average number of Cigarettes smoked in that month. The layout of the dataset is as follows:               ID (Unique Patient Identifier)     Month (1 to 60)     A (Average blood pressure in that month)     B (Average daily exercise hours)     C (Average number of Cigarettes smoked in that month)     Ill (Yes, No)      I was thinking of using Logistic Regression which uses information from last three months and gives a probability for patient to be flagged as Ill in next 2 months. My problem is that logistic regression assumes that observations are independent whereas in my case they are obviously not. What should I do? Should I use something like GEE or GLMM or something else?


Which of the following is a data mining task? Select one: A. Search Google for “Indian National Flag”. B. Giving a 100% correct prognosis of the future success or failure for a new business. C. Grouping search results by theme of the results. To the extend I have understood about DM I have given option B as the answer. But I am not sure about my choice so i am writing here to get the actual right answer.


title: Additive VS dominant model through logistic regression

text: I am working on a regression modelling approach to predict the health status of samples (disease vs no disease). I already have some evidence of the number of genes involved (numeric variable) is significant in predicting the status, and of one particular gene (gene G) being particularly involved in determining disease. The question that I am trying to answer is: since the number of genes is important, is the mechanism of disease determination additive (higher the number of genes altered, higher the likelihood of disease, with gene G being just "stronger" than the others), or dominant (any gene, especially gene G, if altered, could cause disease on its own and other genes hit beside it only act as modifiers)? I was thinking of a logistic regression approach to predict `Status`, checking for interaction between the variable _n.of.genes_ and, for example, gene G being hit or not (dummy variable, T/F). My feeling is that, if the interaction term `n.of.genes*geneG.hit` achieves significance, I could state that the model is dominant...but I am really puzzled and not sure about this being an appropriate strategy. Any suggestion on the way to go? Thank you!


title: Multilevel models: Meaning of random intercept

text: I am quite confused at the moment regarding multilevel models as implemented in SAS PROC MIXED. I am new to mixed models. Let's say I have 2 treatments, a new treatment and a control. And let's say that I have a sample of n subjects, each is contributing 1, 2 or 3 samples, i.e., the treatment (new or control) is being applied in 1 to 3 places in the subject's body. I can't understand the difference between:               Proc mixed data=...     class Treatment SubjectID;     model Y = Treatment;     Random SubjectID;     run;      and               Proc mixed data=...     class Treatment SubjectID;     model Y = Treatment;     Random intercept / subject = SubjectID;     run;      I can't figure out this intercept / slope issue and its meaning.


title: Factor Scores Level 2 in a Multilevel Model

text: My question is twofold and concerns weighting in multilevel models and using factor scores as predictors on level 2. Here is my problem: I want to estimate (using the meologit command in Stata 13) a multilevel model with a ordered categorical dependent variable. My data are persons (N=34500) nested in zip-codes (N=6348) obtained by using a complex survey design. The range of number of observation within level 2-units is [1-45] meaning I have many zip-codes with just 1 participant and several with 20+ person and some up to 45. The average is about 5 persons per zip-code. I combined the data from the interviews with offical statistics which I want to use as level-2 predictors. For instance, I calculated a factor ARMUT (poverty) consisting of 3 items which are theoretical meaningful. To do so I calculated a weight (1/"obs per zip-code") to take into account the distribution of persons over zip-codes. Hence, using the predict command, my factor has a mean=0 and a SD=1. However, when I use unweighted data (which I would prefer and following the _model based_ -approach) to estimate my two-level model, the mean and SD of my level-2 factors are obviously slightly off. OK, now my questions:   1. Do you think it is a good idea to calculate a weight to take into account the size of level-2 units in my factor analysis? If not my problem is solved, but I need a good argument for that approach.   2. How do you think about using weights in multivariate regressions (multilevel modeling in my case)? Any advice would be highly appreciated.


title: Sequential or dynamic multidimensional scaling

text: I am looking for a way to conduct sequential (dynamic) MDS based on several variables coded for a sequence of stimuli. Could anyone direct me to software for this purpose (especially R packages / code)?


title: What kind of sample is this?

text: I took a sample from the university, but I stipulated that the sample students have to be at a specific mental age: 19 at the first year. The number of all students at the 1 year is known, but the size of the students at the age of 19 at the 1 year is unknown, I took the sample through searching for the students who I wanted by non random, every students I found at 1 year I was asking them about their age, and who were 19y I took them. my question is: what we call this sample?


title: Changing variable values and examine the outcome difference between the altered and original data

text: I recently read an approach which is used to find the effect of changing an independent variable. They are doing a classification problem, so each data row (or record) is associated with an outcome of YES or NO. They take one data row (i.e., test row), and then build a model using the other n-1 rows (more specifically, they use k nearest neighbour to build the model). They use the model to compute the probability of the outcome of that one particular test row. They get a value of, for example, 51% chance that the outcome for this data row is YES. They then modify the value of one of the independent variable in that test row by X, and the use the model to examine the outcome again. They get a value of 55%, for example. After the steps, they say that changing the independent variable by X can increase the probability of this test row being YES by 4% (55-51). Since I am new to the field of statistical analysis and data mining...I am not sure whether this analysis approach is sound. I tried to google for other references for this approach, but I could't find any... Could someone please help me explain whether this technique is valid or point me to some references? Thank you very much in advance!


I hope this one is self-explanatory, but let me know if something is unclear: Is there a multivariate version of the Weibull distribution?


title: nls curve fitting of nested/shared parameters

text: I'm trying to fit raw data to curves, which works well on an individual basis. However, I'd like to "share" parameters (sometimes referred as nested parameters) across more than one data series. Is there a way to do this in R?


title: AP Statistics exam

text: I have a good handle on probability theory, at least I understand the book I **1st course in probablity** by Ross. However, I am having a hard time understanding Statistics, I got a prep book to teach me but it seems to be just logical deductions are given as the reason for why $\frac{\sigma}{\sqrt{n}}$ for the st.dev of a sample for a random sample. I am asking for any good books online that would give some more mathematical reason why statistics work. My level: I read and mostly understood the book listed above. I am studying for the AP Statistics exam.


title: Why is the geometric median called the $L_1$ estimator?

text: My question is simply, why is the geometric median called the $L_1$ estimator? This always reminds of $L_p$ spaces but the distance being minimized in the geometric median's definition isn't $L_1$ but rather the $L_2$ (Euclidean) norm. What does the $L_1$ refer to? Is it just a misnomer/an accident/a freak of nature/something historical?


title: Increase in rates over time

text: I am looking for an appropriate method to analyse these data/situation: We sell electronic devices. Customers are encouraged to report defective devices. Around every months we calculate the cumulative rate of {defective_devices}/{all_sold_devices}. I would like to monitor for sharp increases in the rates of defectives; for example if I calculate the cumulative rates each months from 1 to 5 I will have 0.3, 0.55, 0.4, 0.78, 0.85... it appears that the slope of the curve sharply increases at 0.78 onwards. What statistical methods can I use to determine that the slope significantly increases by a certain predetermined value? I'm also worried about the arbitrary selection of monitoring every month, why not every week or even every day? How does one determine what is the proper monitoring period?


title: Marginal effects with logistic generalized additive model in R

text: I am currently working with a logistic semi-parametric model in R using the mgcv package. The output from the model gives the standard log-odds coefficients; however, reviewers have requested marginal effects (like the ones in Stata using the margins command). I would like to do average marginal effects (though, marginal effects at the means--modes for categorical covariates--would at least be a start). I was wondering if anyone has an implementation for this in R for models that use the gam() function. I have a rather large data set (1.2 million observations, a large number of discrete and continuous covariates, and fixed effects for 2,000 individuals). Are these estimates even tractable, given the non-parametric treatment of a couple of the continuous covariates? Any information would be helpful. Here is what I am working with, though I get errors when attempting to do the bootstrapped SEs for the effects (comes from this helpful site probitlogit- marginal-effects-in-r-2/). I am not sure how this treats the non- parametrically estimated smooths (but maybe these don't matter, since it is using the "predict" function?):               mfxboot <- function(modform,dist,data,boot=1000,digits=3){ #dist is the distribution choice of logit or probit           require(mgcv)          x <- gam(modform, family=binomial(link=dist),method="GCV.Cp",data)           # get marginal effects          pdf <- ifelse(dist=="probit",                     mean(dnorm(predict(x, type = "link")))                  mean(dlogis(predict(x, type = "link")))      marginal.effects <- pdf*coef(x)               bootvals <- matrix(rep(NA,boot*length(coef(x))), nrow=boot)       set.seed(1111)       for(i in 1:boot){         samp1 <- data[sample(1:dim(data)[1],replace=T,dim(data)[1]),]         x1 <- gam(modform, family=binomial(link=dist),method="GCV.Cp",samp1)         pdf1 <- ifelse(dist=="probit",                        mean(dnorm(predict(x1, type = "link"))),                        mean(dlogis(predict(x1, type = "link"))))            bootvals[i,] <- pdf1*coef(x1)          }          res <- cbind(marginal.effects,apply(bootvals,2,sd),marginal.effects/apply(bootvals,2,sd))          if(names(x$coefficients[1])=="(Intercept)"){             res1 <- res[2:nrow(res),]         res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep=""),res1)),nrow=dim(res1)[1])              rownames(res2) <- rownames(res1)             } else {         res2 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep="")),nrow=dim(res)[1]))            rownames(res2) <- rownames(res)         }          colnames(res2) <- c("marginal.effect","standard.error","z.ratio")            return(res2)     }


title: What do I do if my predicted values are out of the dependent variable range?

text: I have built a beta regression model with log link for predicting adherence. My dependent variable's range is 0 to 1.When I used a test set to calculate the predicted values with the parameter estimates from this model, the values were out of range (>1 and max was 2.3). Is it wrong if I scale them to 0 to 1 range? It does nt really alter the distribution graph much. What should I do if the predicted values are out of range in general? **Edit** : The Code I am using proc glmmix (sas).               proc glimmix data=_exp0_.Beta_training;      class gndr_cd medispan_drug_class_grp_desc ecom ;      model days_adhered=norm_age gndr_cd drug_type norm_days_until_refill total_fills/ dist=beta solution link=log s;      random residual;      nloptions technique=quanew update=bfgs;      run;


title: RandomForests in Matlab and outliers detection

text: I am solving some regression problem with RandomForests in Matlab, using it's default TreeBagger class for this task. While I managed to get reasonable result already, there are few questions which I can't find answers by simple google search. All questions below are for regression task.   1. `predict` method in TreeBagger class returns predicted value but also it returns standard deviations of separate trees values. It's obvious that the higher this deviation is, the less reliable is result. But I am not sure if it is correct to make any decisions based on these values. So the question is, how standard deviations over the ensemble of trees can be used in practice? Can we do something like prediction outlier detection based on them?   2. I've seen different papers which mention that RF can be used for outlier detection. Is it possible to do that with Matlab's TreeBagger and how? If I am solving regression problem using RF, can I add to this procedure additional step of outlier detection using RF?   3. Let's asume that there are positive answers to my previous question, and I managed to build RF and set up an outlier detection procedure. Now, is it correct to do the outlier test for some new data before trying to predict it's value using built RF model?   4. For my problem I get same result both for forests with 10 and 100 trees. Is there any good reason to choose model with higher number of trees? I always assumed that simpler model should always be chosen to avoid possible overfitting, but in light of my first question, the higher the number of trees the more precise and informative standard deviation of error over ensemble of trees from my first question is. **UPDATE:** Okay regarding my first question, I've tried and removed 25% cases from my validation set with highest stndard deviation of error over ensemble of trees. For remaining 75% of validation data MSE improved by 20%. This means that I am just skipping 25% percent of my data and don't make any prediction at all, but for other data this gives me 20% improvement in prediction quality, which is acceptable for my task. But I still belive that something more clever can be done with those standard deviations.


title: What does a chi-square value of zero means in CFA?

text: I am new to SEM and I am trying to run CFA for a 3-item Job satisfaction scale in Amos. I am getting the following as results.               Chi-square = .000     Degrees of freedom = 0     Probability level cannot be computed      I am not able to understand what does this mean. Please help me understand how do I interpret this. Here is the model: ![enter image description here](http://i.stack.imgur.com/NYtHn.jpg)


title: Terminology for how "grouped" an ordered data set is

text: Apologies in advance for poor terminology / description below; I'm trying my best but I do not know the correct wording. Let's say a "sample" is some fixed number $n$ of boolean values (I'll use 1=true, 0=false):               1 0 0 0 1 0 0 1 1 1 1 1 0      Now I have a list of samples all of size $n$, and the order of items in the list is significant. So, for example, list A:               Sample 1:  0 1 1     Sample 2:  1 0 0     Sample 3:  0 0 1     Sample 4:  0 1 0     Sample 5:  1 0 1     Sample 6:  0 1 0     Sample 7:  1 0 0      Consider also this list, which is the same data as above but in a different order, call it list B:               Sample 1:  1 0 0     Sample 2:  1 0 0     Sample 3:  1 0 1     Sample 4:  0 0 1     Sample 5:  0 1 1     Sample 6:  0 1 0     Sample 7:  0 1 0      Both list A and list B have the same sets of samples, but the order of the samples in list B minimizes "vertical gaps" of 1's ("true" items). My question is: Is there a word that expresses some quantification of how "grouped" individual elements across samples are for a given order of the list? Sorry if that was vague, but what I mean is, is there some word where I can say "List B has a higher/lower ______ than list A"? Another way to state my question is, if I was given list A, and wanted to reorder it to get list B, I'm looking for a word that fits "Reordering list A such that the ______ is minimized/maximized gives list B." From my understanding, the concept is not "entropy", as entropy is independent of order, and so the entropy of list A and list B is the same (if a given "sample" is a value of the variable that entropy is being calculated of). Unless there is some other type of entropy or entropy-like concept that is order-dependent (so XXXXYYYY has less of it than XYXYXYXY). I'm looking for an actual statistics word / concept that I can search for (because I need to analyze this quality of various lists, somehow), not just descriptive English words like "List B is _nicer looking_ than list A". I'm trying to research algorithms to quantify this but I am frustrated because I do not know what to Google for.


title: Bootstrap quadratic term

text: I have bootstrapped a linear model in PROC GLM. Now I want to construct CI's around the parameter estimates. For binary variables, it's straightforward: drop the reference category where the parameter estimates are zero. However, for each replicate, the interaction of the binary term and a quadratic term produces two parameter estimates:                         b2*a     Est            Std err     b2*a       0      -0.0312886      0.02925361          b2*a       1      -0.0339503      0.03308580      It's not clear to me how to add up the replicates to obtain the mean in this case. Normally I use PROC UNIVARIATE and capture the percentile points, but only when there is one value per replicate. I'm not sure what to do. Can someone enlighten me?


title: Distinguishing two datasets

text: I have two datasets from some Web store (like Amazon). Datasets have one and the same structure. Each record in these datasets has the following attributes:               ID - user ID     ProductCnt - count of products user bought     DepartmentCnt - count of departments user was shopping       PayTotal - sum of payments made by user        PayCnt - count of payments     MaxPay - maximum payment      The first dataset is a collection of records related to different users randomly selected. The second dataset is a collection of records related only to users who also clicked on a particular advertisement located on the same page as the product they where shopping for. Problem: Find dependencies that distinguish users in second dataset from users in the first dataset. To solve this I would calculate statistical parameters such as mean, expected value and standard deviation for all parameters in both datasets and compare them.   Any other ideas how to find characteristic features distinguishing these datasets? I am new to this kind of problems, so please bear with me, and also let me know if my question makes no sense at all! Thanks!


title: Classification tasks where machines achieve significantly higher accuracy than humans

text: Spam filtering is an example of classification task where machines achieve approximately similar accuracy as humans. On the other hand, Asirra (recognizing cats and dogs in images) is an example where humans significantly outperform machines. Are there any examples of classification tasks where machines achieve significantly higher accuracy than humans?


title: Interpretation of error bars with covariance

text: I am working on decision making under uncertainty, in which I have four technologies which I should compare them in terms of some decision criteria. These decision criteria are evaluated as probability distribution functions from which I extract two values including modal value and confidence interval values for each technology (min/modal/max). At the end for each criterion I would have an error bar for each technology that I should interpret them. However I have faced a problem regarding interpretation of these error bars. Since the uncertain parameters are changing from one technology into another, I can not claim that all minimum values on these error bars would happen at the same time, and all maximum values would happen at the same time. There is possibility that minimum value of one technology happens when the maximum value of another is happening. As I understood it means that I have covariance. I do greatly appreciate if you can let me know if there is any technique for interpreting these error bars when we do have covariance. Thank you very much in advance for giving your time to help me. Best regards Shabnam


title: Interpreting Linear Model Coefficients with Lagged Variables

text: Let's say I have a data set which looks like below and I'm running a linear model to predict income on two predictors.               date    date-1    income   region    age     mar       apr      50        2        55     apr       may      10        3        40     may       jun      35        1        35       ....          Income = B0 + B1(Region) + B2(Age)      Let's say that I get the following coefficients from the model. How would I interpret such results.                            Est       Intercept    130     Region       0.05     Age          0.10      So with a one month lag, how would I go about interpreting this model. What if it were a two month lag?


title: Robust nonparametric estimation of hazard/survival functions based on low count data

text: We're trying to use a Gaussian process to model h(t) -- the hazard function -- for a very small initial population, and then fit that using the available data. While this gives us nice plots for credible sets for h(t) and so on, it unfortunately is also just pushing the inference problem from h(t) to the covariance function of our process. Perhaps predictably, we have several reasonable and equally defensible guesses for this that all produce different result. Has anyone run across any good approaches for addressing such a problem? Gaussian-process related or otherwise?


title: How Large a Difference Can Be Expected Between Standard GARCH and Asymmetric GARCH Volatility Forecasts?

text: I have been using various GARCH-based models to forecast volatility for various North American equities using historical daily data as inputs. Asymmetric GARCH models are often cited as a modification of the basic GARCH model to account for the 'leverage effect' i.e. volatility tends to increase more after a negative return than a similarly sized positive return. What kind of a difference would you expect to see between a standard GARCH and an asymmetric GARCH forecast for a broad-based equity index like the S&P 500 or the NASDAQ-100? There is nothing particularly special about these two indices, but I think it is helpful to give something concrete to focus the discussion, as I am sure the effect would be different depending on the equities used.


title: Number of cluster samples in GEE

text: In GEE, if you increase the number of clusters would this increase or decrease the standard errors of the regression coefficients? Would the estimated std error of the regression coefficients be smaller or larger than the true std error of the coefficients?


title: Bidding Collusion Question: Probability of Identical Bids

text: I was asked this question. My initial conclusion is that there isn't enough information to calculate the probability. I would appreciate it if anyone could provide provide their insight. Scenario:   Department of Defense is looking to purchase a good. Department of Defense has determined that a fair price is 1.0000 and told vendors of this value. Eight vendors submitted independent bids to provide the desired good. What are the odds/probability that three vendors would submit the exact same bid (see below)?               0.7700     0.7999     1.1120     0.7000     1.4262     0.6840     0.7000     0.7000


title: Interpretation of Type 3 analysis of effects in proc logistic

text: I am trying to best understand how to interpret the type 3 analysis of effects output from `proc logistic` in SAS. I understand it to be interpreted in the same way you would interpret a type 3 sum of squares in a linear regression. For example, I want to model gas price by season (winter,fall,summer,spring) and oil price. I would use the ss3 option in the model statement. I would then get the test for the overall effect of season on gas price given oil price was already in the model and the test for the effect of oil price on gas price given season was already in the model. Can I interpret the type 3 output from proc logistic in the same way? For example, if I wanted to model whether an arcade machine would eat my money or not based on if the machine was a .25, .50, or .75 cents machine and the humidity in the building where the machine was housed.


title: Converting similarity matrix to (euclidean) distance matrix

text: In Random forest algorithm, Breiman (author) constructs similarity matrix as follows:   1. Send all learning examples down each tree in the forest   2. If two examples land in the same leaf increment corresponding element in similarity matrix by 1   3. Normalize the matrix with number of trees He says: > The proximities between cases n and k form a matrix {prox(n,k)}. From their > definition, it is easy to show that this matrix is symmetric, positive > definite and bounded above by 1, with the diagonal elements equal to 1. It > follows that the values 1-prox(n,k) are squared distances in a Euclidean > space of dimension not greater than the number of cases. Source In his implementation, he uses _sqrt(1-prox)_ , where _prox_ is a similarity matrix, to convert it to distance matrix. I guess it has something to do with the "sqaured distances in a Euclidean space"-quoted above. Can somebody shine a little light on why it follows that 1-prox are squared distances in a Euclidean space and why he uses squared root to get distance matrix?


title: Almost sure convergence

text: The problem (not homework) I practiced is Consider the probability space $([0,1], B_{[0,1]}, P)$ where $B_{[0,1]}$ is the Borel set and $P$ is Lebesgue measure on $[0,1]$. For any integer $n>0$, there exist $m$ and $k$ such that $n = 2^m-2+k$ and $0\leq k \leq 2^{m+1.}$ Define $$X_n(\omega)=\left\{\begin{array}{clcr} 1, \mbox{} \frac{k}{2^m}\leq \omega\leq\frac{k+1}{2^m}\\0, \mbox{ otherwise} \end{array}\right.$$ for any $\omega\in[0,1]$. Is the statement, $X_n \to 0$ a.s., correct? I know the solution can be found from p.36, Exercise 47, Mathematical Statistics: Exercises and Solutions by Jun Shao. That is, $X_n(\omega)$ does not converge to $0$ a.s. The reason is following. For any fixed $\omega \in[0,1]$ and m, there exists $k$ with $1\leq k\leq 2^m$ such $\frac{k-1}{2^m}\leq \omega \leq \frac{k}{2^m}$. Let $n_m =2^m-2+k$. Then $X_{n_m}(\omega)=1$. Since m is arbitrary selected we can find an infinite sequence $\{n_m\}$ such that $X_{n_m}(\omega)=1$. This implies $X_n{\omega}$ does not converge to 0. Since $\omega$ is arbitrary, $X_n$ does not converge to 0 a.s. However, I got the different answer. My answer is : Since $P(|X_n| > \epsilon) \leq \frac{1}{2^m}$ for any $\epsilon >0$, $\sum_n P(|X_n| > \epsilon) < \infty$. By the first Borel-Cantelli lemma, $P(\{|X_n| > \epsilon\} \mbox{ i.o.}) = 0$. Since $\epsilon$ is arbitrary, $X_n \to 0$ a.s. If my answer is wrong, I'm curious why the 1st Borel Cantelli lemma does not work or where something wrong is. Thanks for help.


title: What does Determinant of Covariance Matrix give?

text: I am representing my 3d data in covariance matrix. I just want to know what the determinant of a covariance matrix gives. If the determinant is positive, zero, negative, high positive, high negative, what does it mean or represent? Thanks EDIT: Covariance is being used to represent variance for 3d coordiantes that I have. For example, covariance matrix A determinant is +100, and covariance matrix B determinant is +5. Which of these values show if the variance is more? Which value tells that datapoints are more dispersed? Which value shows that readings are further away from mean?


title: Multinomial logit regression with ordinal independent variables

text: I've got some ordinal variables b and a and a categorized variable c. I would like to fit a multinomial logit regression from the library car. I tried to ignore the ordinal scale. I have the following data:                a<-c( 3, 4,   4,   4,   3, 4,   3, 3, 4,   2, 2, 4,   3, 3, 3, 1,   3, 2, 2, 3, 3, 1,   3, 2, 2, 3, 2, 3, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 3, 3, 2, 2, 3, 1,   2, 2, 2, 2, 3, 2, 3, 4,   4,   3, 3, 2, 2, 3, 3, 3, 2, 1,   1,   1,   1,   1,   1,   2, 3, 4,   3, 3, 4,   3, 4,   3, 2, 3, 3, 3, 3, 3, 4,   3, 4,   3, 3, 2, 3, 3, 3, 3, 2, 2, 3, 3, 3, 2, 3, 2, 3, 1,   2, 2, 1,   1,   4,   3, 3, 2, 2, 2, 2, 2, 2, 3, 4,   4,   4,   3, 3, 3, 3, 3, 4,   4,   3, 3, 2, 3, 3, 3, 3, 4,   3, 4,   2, 2, 3, 3, 3, 2, 2, 3, 2, 4,   2, 2, 2, 2, 2, 1,   2, 2, 1,   1,   3, 4,   3, 3, 2, 3, 3, 3, 3, 2, 3, 3, 3, 2, 1)                c<-c(5 ,3 ,4 ,3 ,4 ,4 ,2 ,2 ,3 ,4 ,2 ,5, 3, 5, 4 ,3 ,2 ,4 ,4 ,4, 4 ,4, 4, 2 ,3, 4 ,2 ,3 ,3 ,3 ,4 ,3 ,3 ,2 ,2 ,3 ,3 ,3 ,3 ,4 ,2 ,4 ,3, 3, 3, 4, 4, 3, 3 ,2 ,3 ,3 ,3 ,3, 4 ,4, 4, 3, 2, 2 ,4 ,4 ,3 ,3 ,2 ,2 ,1 ,2 ,2 ,2 ,1 ,2 ,5 ,2 ,3 ,3 ,2, 4 ,3 ,1 ,2 ,3 ,2 ,3 ,3 ,3 ,3 ,3 ,3 ,2 ,2 ,2 ,2 ,3 ,2 ,4 ,3, 3 ,2 ,3, 2, 4, 3, 3, 3 ,3 ,4 ,2 ,2 ,4 ,3 ,3 ,3 ,3 ,3 ,2 ,3, 3 ,3, 3, 4 ,4 ,4 ,1 ,3 ,3 ,3 ,4 ,4 ,4 ,3 ,2 ,4 ,4 ,2 ,4 ,4 ,4 ,4 ,2 ,3 ,3, 2, 2 ,3 ,2 ,3 ,4 ,5 ,2, 3 ,3 ,2 ,3 ,2 ,2 ,3 ,2 ,2 ,4 ,4 ,3 ,3 ,2, 4 ,4 ,2 ,4 ,3 ,4, 4, 3 ,2 ,3)     b<-c(3 ,2, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 4, 2, 3, 1, 3, 1, 2, 4, 2, 1, 3, 2, 2, 2, 1, 3, 3, 3, 2, 2, 2, 2, 1, 3, 1, 3, 2, 3, 1 ,3 ,3 ,2, 2, 3, 1, 3, 2, 2, 2, 2, 2, 2, 3, 4, 3, 3, 2, 1, 4, 3 ,3 ,2 ,2, 1, 2, 2, 2, 2, 1, 2, 5, 3, 3, 4, 3, 4, 1, 2, 3, 3, 3, 3, 2, 3, 3, 3, 3, 2, 3, 2, 2, 3, 2 ,4 ,2 ,3 ,2 ,2 ,2 ,4 ,2 ,2 ,2 ,2 ,2 ,2, 1, 5 ,4 ,3 ,2, 2 ,2 ,2 ,2 ,4 ,2 ,2 ,4 ,3 ,3 ,1 ,2 ,2 ,2 ,2 ,2 ,2, 2, 2, 2, 2, 2 ,2 ,2 ,2 ,2 ,2 ,2, 2, 2, 2 ,2 ,2 ,2 ,2 ,5 ,4 ,3 ,2 ,1 ,1 ,1 ,4 ,3 ,2 ,2 ,3 ,3 ,3 ,2 ,2 ,2 ,2 ,2, 3, 2 ,2 ,2 ,2 ,2 ,1)      now I ignored the ordinal scale and treated them as factors to fit the multinomial logit regression               require(car)     a<-as.factor(a)         b<-as.factor(b)     c<-as.factor(c)     multinom(formula = a ~ b + c)          Call:     multinom(formula = a ~ b + c)          Coefficients:       (Intercept)        b2       b3       b4        b5         c2         c3         c4        c5     2   0.3410779  1.009797 41.80056 45.22081 -13.02923 -0.5229982  0.9216514  0.2170273 -18.03928     3  -1.4697131  2.698228 44.91938 47.04268 -16.24570 -0.7341395  0.7088424  1.2495310  20.70641     4 -46.0095393 33.603384 75.13911 79.00502  56.91264 -7.4198320 13.0220759 14.2526951  33.85774          Std. Errors:       (Intercept)        b2        b3        b4           b5           c2        c3        c4           c5     2   1.2654428 0.6530052 0.4659520 0.5495402          NaN 1.337075e+00 1.4180126 1.4993079 8.028986e-16     3   1.6649206 0.9361438 0.5123106 0.5879588 2.446562e-15 1.640462e+00 1.7003411 1.7558418 8.601766e-01     4   0.3399454 0.4767032 0.3699569 0.4144527 3.321501e-11 6.973173e-08 0.6549144 0.6953767 8.601766e-01          Residual Deviance: 328.1614     AIC: 382.1614        I think I found the mistake....the column b5 is empty for a1 and a2.               table(b,c,a)     , , a = 1             c     b    1  2  3  4  5       1  0  3  2  2  0       2  1  7  1  0  0       3  0  0  0  0  0       4  0  0  0  0  0       5  0  0  0  0  0          , , a = 2             c     b    1  2  3  4  5       1  1  5  2  2  0       2  1 12 21  4  0       3  0  1  6  1  0       4  0  2  1  1  0       5  0  0  0  0  0      But do you know how to solve this problem?


title: Handling missing data (holidays) in multiple time series (historical simulation VaR)

text: I have two time series of daily returns on two stock indices (S&P 500 and BOVESPA) that I would like to estimate the portfolio value at risk (VaR) for. Since these are indices from two different regions they do not naturally line up smoothly as their business days are unequal due to some country specific holidays. An example would be on the 4th of July, when the return of the S&P would be 0 and the BOVESPA a positive or negative number. Including this sample in the empirical returns distribution would include a 0 return which never happened. As I am interested in the portfolio VaR I would prefer returns on days where both markets were open so as to catch the implicit correlations in the returns of the indices. However, this would mean I will remove 250 paired observations (over a 20 year period) where either one of them have a zero value. Are there any better ways to handle the data?


title: How can I fit the parametres of t-distribution if I know the number of degrees of freedom in R?

text: I have some data. Does anybody know a function which can fit the parameters of a t-distribution of this data when I know the numbers of degrees of freedom? Thanks!


title: Can you use the chi-squared test when the expected values are not determined?

text: I have six ratios (females-males). These are the observed values: 1.78 - 1.17 - 0.53 - 1 - 0.85 - 0.56 I need to run a chi-squared test to show that these ratios are significantly different. I expect all these ratios to be equal. But as for the “equal” value, there is no expectation. It could be any value, providing that this value is the same for all the cases. How can I run a chi-squared test in this condition? What do I plug into the “expected value” of the formula? Thanks These are the data:               Group    1    2    3    4    5    6        Total          Men :    9   17   13   12   11   19         81     Women:  16    9   11   14   11    5         66          Total:  25   26   24   26   22   23        147


title: Which W to report for the Wilcoxon test?

text: Imagine that I have the following data:               GroupA <- c(2,4,3,1,2,3,3,2,3,1)     GroupB <- c(3,5,4,2,4,3,5,5,3,2)      Then when I call in R               wilcox.test(GroupA, GroupB)      and then               wilcox.test(GroupB, GroupA)      I get different Ws (but the same p-value). So my questions are: Why is this happening? (am I doing something wrong here?) What is the correct interpretation for this? Which W should I choose when reporting the results? Thanks


title: Correlation of categories

text: I am comparing the Facebook posts of a TV news organization's Facebook Page to the contents of the TV news broadcast. (This is called "agenda setting" in communications research.) I have used a grounded theory approach to code hundreds of Facebook posts and segments of the TV show. Here is the problem: There are about 1500 Facebook posts with 500 code categories and 51 TV segments with 300 code categories. When I test the correlation, I have been testing the top 20 Facebook post categories against the number of times the same categories appear on the TV show. Some people have suggested that I need to test all of the Facebook categories against all of the TV show categories. Is this true? They are not equal in number. As a compromise I have converted the number of times a code appears in Facebook or on TV to percentage - percentage out of 1500 for Facebook and percentage out of 51 for the TV show. Does this seem reasonable? Is there a better way to do this without comparing the 500 to the 300 codes. I am using Pearson's correlation coefficient.


title: Why to use fixed regressor in Breusch-Pagan?

text: I need to use Breusch-Pagan to check heteroscedasticity of my timeseries. I found the explanation of the function, here: http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/lmtest/html/bptest.html Looking at the example, I saw that the regressor is a simple serie with 1 -1 repeated values. Could someone explain me why ? Do I have to use fixed values on the regressor? However if I have a vector(R vector) that I would test for heteroscedasticity I could do:               x <- rep(c(-1,1), 375)     mySerie <- rnorm(750)     mySerie <- 1 + x + mySerie     bptest(y ~ mySerie)      Right?   Thank you!


title: How do I compare bootstrapped regression slopes?

text: Let us assume I have two data sets with _n_ observations of data pairs of independent variable _x_ and dependent variable _y_ each. Let us further assume I want to generate a distribution of regression slopes for each data set by bootstrapping the observations (with replacement) _N_ times and calculating the regression _y_ = _a_ + _bx_ each time. How do I compare the two distributions in order to say the slopes are significantly different? A U-test for testing the difference between the medians of the distributions would be heavily dependent on N, that is, the more often I repeat the bootstrapping the more significant will be the difference. How do I have to calculate the overlap between the distributions to determine a significant difference?


title: Using statistical significance test to validate cluster analysis results

text: I am surveying the use of statistical significance testing (SST) to validate the results of cluster analysis. I have found several papers around this topic, such as "Statistical Signiﬁcance of Clustering for High-Dimension, Low–Sample Size Data" Yufeng LIU, David Neil HAYES, Andrew NOBEL, and J. S. MARRON "On some significance tests in cluster analysis" H. Bock but I am interested in **finding some literature arguing that SST is NOT appropriate to validate results of cluster analysis**. The only source I have found claiming this is a web page of a software vendor: http://www.statsoft.com/textbook/cluster-analysis/ To clarify: I am interested in testing whether a significant cluster structure has been found as a result of cluster analysis, so, I'd like to know of papers supporting or refuting the concern "about the possibility of post-hoc testing of the results of exploratory data analysis used to find clusters". I've just found a paper from 2003 "Clustering and classification methods" by Milligan and Hirtle saying, for example, that using ANOVA would be an invalid analysis since data have not have random assignments to the groups.


title: Inferring a Gaussian from noisy data

text: Assume a noise comes from a specific point on a line, noise which I can detect but not completely accurately. My uncertainty is assumed to be Gaussian. I want to gather evidence about the real location of the noise source from the noise I can perceive. So, I represent my accumulating knowledge with a Gaussian variable with an initial value of $\mathcal{N}(0, 1000)$ to represent the mean (starting at $0$, as it must be a positive number) and variance ($1000$, to represent lack of initial knowledge of the actual noise source). The real location of the noise source will be on a point within the noises I can detect. Assuming the Gaussian variable represents my growing knowledge of the noise position on the line, what is the best way to infer the actual noise position from the perceived data?


title: Multiple Regression and number of parameters to include for a learning algorithm

text: I am quite new to Machine Learning and come from a computing background.   1. I have a quite big set of features (~50) with about 4k observations. Is it correct thinking to include all of them in a multiple linear regression model? Wouldn't that possibly introduce overfitting? If that's the case, is it safe to manually identify which of the features would matter the most for the model and for future predictions, and include only those?   2. Given prior experience with Java, would you recommend any other tools than Java for Machine Learning computations? I am aware of some Java libraries such as Weka, Mahout and other languages such as Python, MATLAB, and R (which I hear is quite slow).


title: How to calculate marginal effects for categorical covariates using mlogit in R

text: I am trying to use the `mlogit` package in R and have been following the vignette trying to figure out how to get the marginal effects for my data. The example provided uses continuous variables, but I am wondering how to do this with categorical explanatory variables. I have a value of risk which is continuous as a covariate, but I also have age, class, and gender as covariates. I want to see the marginal effects of "females" only or of "Young - females" in regard to risk. How would I do this? The help documents say:               z <- with(Fish, data.frame(price = tapply(price, index(m)$alt, mean),                                    catch = tapply(catch, index(m)$alt, mean),                                income = mean(income)))     # compute the marginal effects (the second one is an elasticity     effects(m, covariate = "income", data = z)     effects(m, covariate = "price", type = "rr", data = z)     effects(m, covariate = "catch", type = "ar", data = z)      I'm not sure how to manipulate the `z` data frame to get the mean risk for females or young females to then be able to calculate the marginal effects. Would I do them all separately? Do I somehow divide the data frame by age class (say I have just 2 age classes: young and old) so that I have 1 data frame for the young, and a separate new data frame for old, then calculate mean risk? What I am hoping to get from my own data is to be able to interpret the magnitude of the likelihood in producing my categories of offspring. As an example, what I want to say is that if there is a 1 unit increase in risk, it is 10% more likely for older females to produce 2 offspring. As there is a 1 unit increase in risk, younger females are 15% more likely to produce 2 offspring. I am not sure how to calculate the marginal effects by hand, and therefore am confused as to how to get a package to do it for me.


title: lme4 or other open source R package code equivalent to asreml-R

text: I want to fit mixed model using lme4, nlme, baysian regression package or any available. _**Mixed model in Asreml- R coding conventions_** before going into specifics, we might want to have details on asreml-R conventions, for those who are unfamiliar with ASREML codes.               y = Xτ + Zu + e ........................(1) ;       the usual mixed model with, y denotes the n × 1 vector of observations,where τ is the p×1 vector of ﬁxed eﬀects, X is an n×p design matrix of full column rank which associates observations with the appropriate combination of ﬁxed eﬀects, u is the q × 1 vector of random eﬀects, Z is the n × q design matrix which associates observations with the appropriate combination of random eﬀects, and e is the n × 1 vector of residual errors.The model (1) is called a linear mixed model or linear mixed eﬀects model. It is assumed ![enter image description here](http://i.stack.imgur.com/gxdur.jpg) where the matrices G and R are functions of parameters γ and φ, respectively. The parameter θ is a variance parameter which we will refer to as the scale parameter. In mixed eﬀects models with more than one residual variance, arising for example in the analysis of data with more than one section or variate, the parameter θ is ﬁxed to one. In mixed eﬀects models with a single residual variance then θ is equal to the residual variance (σ2). In this case R must be correlation matrix. Further details on the models are provided in the Asreml manual (link). Variance structures for the errors: R structure and Variance structures for the random eﬀects: G structures can be specified. ![enter image description here](http://i.stack.imgur.com/or4Gj.jpg)![enter image description here](http://i.stack.imgur.com/oXTgc.jpg) variance modelling in asreml() it is important to understand the formation of variance structures via direct products. The usual least squares assumption (and the default in asreml()) is that these are independently and identically distributed (IID). However, if the data was from a field experiment laid out in a rectangular array of r rows by c columns, say, we could arrange the residuals e as a matrix and potentially consider that they were autocorrelated within rows and columns.Writing the residuals as a vector in field order, that is, by sorting the residuals rows within columns (plots within blocks) the variance of the residuals might then be ![enter image description here](http://i.stack.imgur.com/SPE5b.jpg) ![enter image description here](http://i.stack.imgur.com/IcikW.jpg) are correlation matrices for the row model (order r, autocorrelation parameter ½r) and column model (order c, autocorrelation parameter ½c) respectively. More specifically, a two-dimensional separable autoregressive spatial structure (AR1 x  AR1) is sometimes assumed for the common errors in a field trial analysis. _**The example data:_** nin89 is from asreml-R library, where different varities were grown in replications / blocks in rectangular field. To control additional variability in row or column direction each plot is referenced as Row and Column variables (row column design). Thus this row column design with blocking. Yield is measured variable. **Example models** I need something equivalent to the asreml-R codes: The simple model syntax will look like the follows:                rcb.asr <- asreml(yield ∼ Variety, random = ∼ Replicate, data = nin89)        .....model 0      The linear model is specified in the fixed (required), random (optional) and rcov (error component) arguments as formula objects.The default is a simple error term and does not need to be formally specified for error term as in the model 0. here the variety is fixed effect and random is replicates (blocks). Beside random and fixed terms we can specify error term. Which is default in this model 0. The residual or error component of the model is specified in a formula object through the rcov argument, see the following models 1:4. The following model1 is more complex in which both G (random) and R (error) structure are specified. **Model 1:**               data(nin89)                # Model 1: RCB analysis with G and R structure          rcb.asr <- asreml(yield ~ Variety, random = ~ idv(Replicate),            rcov = ~ idv(units), data = nin89)      This model is equivalent to above model 0, and introduces the use of G and R variance model. Here the option random and rcov specifies random and rcov formulae to explicitly specify the G and R structures. where idv() is the special model function in asreml() that identifies the variance model. The expression idv(units) explicitly sets the variance matrix for e to a scaled identity. _**# Model 2: two-dimensional spatial model with correlation in one direction_**                 sp.asr <- asreml(yield ~ Variety, rcov = ~ Column:ar1(Row), data = nin89)      experimental units of nin89 are indexed by Column and Row. So we expect random variation in two direction - row and column direction in this case. where ar1() is a special function specifying a first order autoregressive variance model for Row. This call specifies a two-dimensional spatial structure for error but with spatial correlation in the row direction only.The variance model for Column is identity (id()) but does not need to be formally specified as this is the default. _**# model 3: two-dimensional spatial model, error structure in both direction_**                sp.asr <- asreml(yield ~ Variety, rcov = ~ ar1(Column):ar1(Row),        data = nin89)     sp.asr <- asreml(yield ~ Variety, random = ~ units,       rcov = ~ ar1(Column):ar1(Row), data = nin89)      similar to above model 2, however the correlation is two direction - autoregressive one. I am not sure how much of these models are possible with open source R packages. Even if solution of any one of these models will be of great help. **_Even if the bouty of +50 can stimulate to develop such package will be of great help !_** _**See MAYSaseen has provided output from each model and data (as answer) for comparision._** _**Edits: The following is suggestion I received in mixed model discussion forum:_** " You might look at the regress and spatialCovariance packages of David Clifford. The former allows fitting of (Gaussian) mixed models where you can specify the structure of the covariance matrix very flexibly (for example, I have used it for pedigree data). The spatialCovariance package uses regress to provide more elaborate models than AR1xAR1, but may be applicable. You may have to correspond with the author about applying it to your exact problem."


title: Why am I getting this result in modeling a Pareto Type II distribution in Excel?

text: In Excel for a project I'm trying to model the density, distribution and survival function ($1-F(X)$) and I can't get the density to sum to one and I can't get the distribution to go to one. For the parameters I'm using $\alpha=3$ and $\theta=10$. I'm getting extremely strange results such that halfway through the distribution it levels of at $.36366$ until the end. The density function is summing to $1.15...$. After checking the formula that I input, I know it's right, but is there something I should know about modeling this in excel that maybe the program is picky about? I'm not sure if this is enough information, if not comment on what you need to know and I'll get back to you.


title: k-means clustering on percentages

text: Can we do _k_ -means clustering on percentage data (like 56%, 44%, 22%, 13%, etc.)?   There is a data set, and data in various parts are measured in percentages.


title: Trusting the results of a binary vector comparison

text: I have two pretty long binary vectors ($\sim 3$ mi positions), and I found the distance between then using Jaccard's Index $J=\frac{a}{a + b + c}$. After getting the result, I asked myself how much can I trust this value. I thought about two alternatives, and would like to hear your advice to be sure that they are appropriate.   1. Randomize both vectors multiple times and check with z-score or t-test, where my original result is falling within this random distribution.   2. Resample with replacement (Bootstrap) the vectors and verify if my original result is always within the 95% Confidence Interval. The issue with the first alternative is that even after 10.000 randomizations, many times the original result is very far from the distribution mean (z-score $= \sim 20$). =========================================== Here are more details about the problem at hand. The binary vectors are indicators of change in my samples - one indicates a change, zero otherwise. One vector represents one characteristic, and the other vector another characteristic. They are in principle varying independently and subject to distinct external influence. However, my objective is to find out if these two characteristics are in fact changing together - i.e., when the characteristic represented by one vector changes, the other characteristic changes as well. And importantly, how often could this " _vector similarity_ " occur by chance. For instance, please consider: Vector A - `0001111` Vector B - `0001011` and Vector C - `01010` Vector D - `10100` In this case, vectors A and B seem to be pretty related, but vectors C and D not so much. Now I'm trying to figure out rigorous ways to assess how this vector similarity could arrive by chance.


title: Case-control study with data collected in batches

text: I have (matched) case-control data. The data is collected in batches in such a way that the batch determines some quality of the data (there is a 'batch' effect). Also, the cases and controls were not randomly distributed among these batches. So batch 1 contains only cases, batch 2 contains only controls, batch 3 is mostly cases (a couple controls) and batch 4 is mostly controls (a couple cases). I'm trying to do logistic regression to determine if some covariates are associated with disease. I need to control for the batch effect. Will the fact that the cases/controls are not randomly distributed among the batches have a negative impact on my model if I include the batch as a covariate? EDIT: From my reading it appears that can be somewhat of a common problem called Quasi or Complete separation.


title: Logistic Regression - Predicting an event with a couple time-related issues

text: I'm using logistic regression to predict the occurrence of tree carcasses falling after a mortality event. Data include a variety of topographic and tree characteristic variables and also time since death for many individuals over a number of sites. Data were taken repeatedly over several years. The way things were presented to me, most variables are treated as fixed effects, with 'site' as the only random effect. However, I am wondering if I need to treat the time-related variables in a special way. Specifically,   1. Whether I should include time as a random effect, or in some sort of repeated-type setup.   2. Also, data exist for trees after they fell. I am unsure if I need to come up with a special to deal with this. Obviously a tree stayed on the ground once it got there. I'm not sure if it matters, but I will be running this in glimmix in sas.


title: Confusion related to the bagging technique

text: I am having a bit of confusion. I was reading this paper where it explained that bagging technique greatly reduces variance and only slightly increases bias. I didn't get it how come it reduces variance. I know what variance and bias is. Bias is the inability of the model to learn the data. Variance is something similar to overfitting. I just don't get it how bagging reduces variance.


title: Specification and interpretation of interaction terms using glm()

text: I am fitting a logistic model to data using the glm function in R. I have attempted to specify interaction terms in two ways:               fit1 <- glm(y ~ x*z, family = "binomial", data = myData)      fit2 <- glm(y ~ x/z, family = "binomial", data = myData)       I have 3 questions: 1) What is the difference between specifying my interaction terms as x*z compared to x/d? When I call summary(fit1) the report includes results for the intercept, x, z, and x:z while summary(fit2) only includes results for intercept, x, and x:z. I did look at Section 11.1 in "An Introduction to R" but couldn't understand the meaning. 2) How do I interpret the fit equation mathematically? Specifically, how do I interpret the interaction terms formulaically? Moving to math instead of R, do I interpret the equation as: logit(y) = (intercept) + (coeff_x)*x + (coeff_z)*x + (coeff_xz)*x*z ? This interpration may differ in the two specifications fit1 and fit2. What is the interpretation of each? 3) Assuming the above interpretation is correct, how to I fit the model of x*(1/z) in R? Do I need to just create another column with these values?


title: Interpret regression coefficients after WLS

text: I have a couple question that have arisen studying Weighted Least Squares from Woolridge, all centering around corrections for heteroskedisticity in cross sectional data.   1. Woolridge first shows an example of weighted least squares assuming the variance function is known. He divides the equation $$sav_{i}=\beta_{0}+\beta_{1}inc_{i}+u_{i}$$ by $\sqrt{inc_{i}}$ assuming that $\text{var}(u|x) = \sigma^2inc_{i}$ to arrive at the equation $$\frac{sav_{i}}{\sqrt{inc_{i}}}=\beta_{0}(\frac{1}{\sqrt{inc_{i}}})+\beta_{1}\sqrt{inc_{i}}+u_{i}$$ **Question:** Do you now have to interpret $\beta_{1}$ in terms of the square root of income on the ratio of savings to the square root of income?   2. When using FGLS (not knowing the variance function) to correct for heteroskadasticity, is the interpretation of the betas the same as in the ordinary regression (just now they are not unbiased) - i.e. with $sav_{i}=\beta_{0}+\beta_{1}inc_{1}+\beta_{2}X_{2}+u_{i}$ that for each unit increase in income, savings increases by $\beta_{1}$, holding $X_{2}$ (whatever it is, but was used in the weighting - see FGLS algorithm below) constant?   3. Are most econometricians these days using OLS for the unbiased point estimates and making the correction for heteroeskadacity by using robust standard errors?  **NOTES:** In both #1 and #2, you are using weighted least squares and getting there by dividing both sides of the equation by $\sqrt{weight_{i}}$. **In case #1** , you assume the variance is known and is proportionate to income (i.e. inc). **In case #2** , you are estimating the variance by the following process: ![enter image description here](http://i.stack.imgur.com/NFNa6.jpg)


title: Matched binomial pairs

text: I am trying to analyse the following table: 5 5 5 5 4 5 5 5 5 3 4 3 3 4 4 3 5 5 Each number is a success score out of 5 attempts, so we can assume binomial (p, 5), for unknown p. Each column pertains to one person. The bottom row contains the "pre" test results and the top row is the "post" test results. At a quick glance, every item in the top row is at least as large as the corresponding item in the bottom row, so by a quick sign test, the hypothesis that there is a difference pre and post should be significant. If I do a paired t-test on these data (holding my nose, because each item is binomial (p,5)) I get a significant result. If I put the data into a 2 x 2 contingency table (pre post x right wrong) and ignore the matching, I get a significant result. But when I do a contingency test of the table above, testing for independence conditional on the margins, the result is non-significant (p = 0.97). I used fisher.test in R, the Fisher exact test. So what is the best way to test the hypothesis that there is a significant difference, pre and post, taking matching into account? And why does the analysis of independence not apply here?


title: Text feature vector extraction

text: I have a class assignment to implement a couple existing ways to extract feature vectors from a given set of texts, so they can be used to classify those texts using k-nearest neighbour algorithm. The texts are newspaper articles from Reuters, and they are to be sorted based on the country related to the story, and its subject (economy, politics, etc.). What existing text feature extraction algorithms would work well for this task?


title: Multiple comparisons for parameter significance in a sparse high dimensional regression model

text: I have a regularized linear regression model with a large number of parameters (~100,000) that has been optimized using a sparse fitting algorithm (coordinate descent with early stopping). The model parameters and prediction accuracy were cross-validated using repeated random subsampling. I want to assess the statistical significance of individual parameters using the t-ratio (mean coefficient divided by standard deviation of the resampled estimates). Is correction for multiple comparisons appropriate in this case and if so is it appropriate to select only the non-zero weights for significance testing? Correction for multiple comparisons of individual regression parameters does not seem to be a common practice (e.g., correcting for testing significance of multiple main effects in ANOVA).


title: Balancing a skewed dataset in R

text: I am trying to classify a dataset which has 20 features and a class (target) feature. The problem is that "class-1" is found in 25% of cases while the rest 75% is "class-2" examples. My question is : How can i under sample instances of class2 (or oversample class1 instances) in R? I already tried to look for some answers but i would not like to use weighting. Thank you


title: How can a standardized canonical coefficient be greater than 1.0 in a discriminant function analysis?

text: I have a data that I have used discriminant function analysis with. In the results, one variable has a standardized canonical function coefficient that is greater than 1.0. I didn't think these could be greater than 1 and I am not sure what to make of it.


title: Analysis of a complex sample with replicate weights: error when using survey package

text: I have a data base from a stratified multi-stage survey. Due to confidentiality reasons, they do not gave me the strata neither the psu, ssu… variables. Instead they gave me the:   * “Final weights” which sum is approximately equal to N. these weights include the deff, a non response weight and they are also calibrated for sex and strata. (I will call these weights “FinalW”)   * “Replicate weights” for Jacknife method (I will cal this file “repw”)   * “Coeficientes” for Jacknife method (I will cal this file “jkcoef”)   * “Normalized weights” with mean equal to 1 (I don’t know why it is equal to 1). According to the available information, these weights are the “Final weights” times the ratio between sample size and population size by strata and sex. (I will cal these weights “NormW”) To define the design I used the svrepdesign: design_name<-svrepdesign(data=name_data_base,repweights= repw, weights=~FinalW, type="JKn",rscales=jkcoef,scale=1,combined.weights=F) I tried to compute several t-test for independente samples using svyttest and svy.glm. Example: svyttest(NumericVariable~FactorVariable, design_name) summary(svyglm(NumericVariable~FactorVariable, design=design_name)) With both instructions I get the error:   “Error in crossprod(x, y): requires numeric/complex matrix/vector arguments” Actualy I get this error also with svychisq, svyglm,… Can anybody tell me what I’m doing wrong? Do these functions work with replicate weights? Or should I use the usual instructions (outsider survey package) with normalized weights? What are the consequences of such a procedure or it is equivalent to use the functions of survey package? I will be very grateful if somebody could help me. Thanks in advance. Please, apologize my English but I’m not native in English.


title: Simple Linear Regression in R

text: I have answered all of the following questions I need someone to verify me or if there is a better approach I would like to know about. > Q1) Look at your model summary to find the x variable whose model > coefficient is most significantly different from 0. (You don't have to write > R code to find this other variable–just read your model summary.)   >  Q2) Make a simple linear regression model for PBE vs. this x.   >  Q3) Make a scatterplot of PBE vs. this x.   >  Q4) Add the simple regression line to your scatterplot.   >  Q5) Include a reasonable title and axis labels. * * *               A0) beeflm = lm(PBE ~., data=beef)          A1) (coef(beeflm))       (Intercept)          YEAR           CBE      2693.01348650   -1.28693774   -1.84919910                PPO           CPO           PFO        -0.99901169   -1.73045916    1.27410503               DINC           CFO         RDINC        -2.49792219    1.04485422    1.32154103                RFP        -0.01729997           A2) PBEvsDINC = lm(PBE~DINC, data=beef)          A3,A5) plot(beef$PBE, beef$DINC, main="Beef PBE vs DINC", xlab="DINC", ylab="PBE")        (Am I right?)      Error in plot.new() : figure margins too large          A4) abline(PBEvsDINC)      And finally here's the question I don't know the answer: > Q: Are the coefficients (y-intercept and slope in the x direction) the same > for this second simple linear regression model as they are in the first > multiple regression model? Here's the data for beef data.frame:               YEAR    PBE     CBE     PPO     CPO     PFO     DINC    CFO     RDINC   RFP     1925    59.7    58.6    60.5    65.8    65.8    51.4    90.9    68.5     877     1926    59.7    59.4    63.3    63.3      68    52.6    92.1    69.6     899     1927      63    53.7    59.9    66.8    65.5    52.1    90.9    70.2     883     1928      71    48.1    56.3    69.9    64.8    52.7    90.9    71.9     884     1929      71      49      55    68.7    65.6    55.1    91.1    75.2     895     1930    74.2    48.2    59.6    66.1    62.4    48.8    90.7    68.3     874     1931    72.1    47.9      57    67.4    51.4    41.5      90      64     791     1932      79      46    49.5    69.7    42.8    31.4    87.8    53.9     733     1933    73.1    50.8    47.3    68.7    41.6    29.4      88    53.2     752     1934    70.2    55.2    56.6    62.2    46.4    33.2    89.1      58     811     1935    82.2    52.2    73.9    47.7    49.7      37    87.3    63.2     847     1936    68.4    57.3    64.4    54.4    50.1    41.8    90.5    70.5     845     1937      73    54.4    62.2      55    52.1    44.5    90.4    72.5     849     1938    70.2    53.6    59.9    57.4    48.4    40.8    90.6    67.8     803     1939    67.8    53.9      51    63.9    47.1    43.5    93.8    73.2     793     1940    63.4    54.2    41.5    72.4    47.8    46.5    95.5    77.6     798      ##  1941 56 60 43.9 67.4 52.2 56.3 97.5 89.5 830 Here is the t-value results, so is CBE the most significant variable according to this?:               summary(lm(formula = PBE ~ ., data = beef))          Call:     lm(formula = PBE ~ ., data = beef)     Residuals:     Min 1Q Median 3Q Max      -1.1986 -0.5658 -0.1001 0.7067 1.3251     Coefficients:     Estimate Std. Error t value     (Intercept) 2693.0135 1337.9523 2.013     YEAR -1.2869 0.6970 -1.846     CBE -1.8492 0.2834 -6.525     PPO -0.9990 0.4051 -2.466     CPO -1.7305 0.5685 -3.044     PFO 1.2741 1.9645 0.649     DINC -2.4979 2.6018 -0.960     CFO 1.0449 1.0678 0.979     RDINC 1.3215 1.7117 0.772     RFP -0.0173 0.1131 -0.153     Pr(>|t|)      (Intercept) 0.084020 .      YEAR 0.107337      CBE 0.000326 ***     PPO 0.043064 *      CPO 0.018743 *      PFO 0.537314      DINC 0.368995      CFO 0.360411      RDINC 0.465337      RFP 0.882717      ---     Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1     Residual standard error: 1.238 on 7 degrees of freedom     Multiple R-squared: 0.986, Adjusted R-squared: 0.968      F-statistic: 54.77 on 9 and 7 DF, p-value: 1.165e-05


title: Comparing power law fits with large uncertainties

text: I want to test how well my data fit the particular power law: $y=ax^b$ where $b$ , for physical reasons, should equal exactly $-0.5$. I would like to find the probability that the data do not obey this, or more correctly what the chance is of obtaining the data if it did equal $-0.5$. Here's some typical data with estimated $2\sigma$ uncertainties:                  x =(   5,  15,  13,  60, 125, 505)     delx =(   3,  14,   8,   2,   2,   3)        y =( 415, 316, 168,  59,  29,   5)     dely =(  33,   3,   6,   4,   3,   1)      I hate to fit a straight line on a loglog plot in this neighbourhood, but the errors are so poorly estimated it hardly matters. The logging also biases the calculated power, but that should be ok as it just means that the $p$-value will be a maximum. The problem is that most of the error is in the $x$-axis-- is it permissible to swap $x$ and $y$ so a weighted linear regression takes account of the bulk of the error? I.e. instead of $\log(y)=\log(a) + b \cdot log(x)$ use the model $\log(x)=\frac{1}{b} \log(y) -\frac{1}{b} \log(a)$ Does a $t$-test on the power and its resulting $p$-value actually mean anything when there are so many violated assumptions? Also, would leave-one- out-cross-validation be a defence against the small sample?


title: For the binomial distribution, why does no unbiased estimator exist for $1/p$?

text: Suppose that $X$ ~ $Binomial(n,p)$ for $0 < p < 1$ Why does no unbiased estimator exist for $1/p$? My approach: We try to find the structure of $E_p(U(x))$, where $U(x)$ is any estimator of $1/p$. Now, we will have: $\sum{U(x)\binom{n}{x}p^x(1-p)^{n-x}}<\sum{U(x)\binom{n}{x}}=M(n)<\infty$ so that the expectation is bounded above. So this is supposed to mean that if $p < 1/M(n)$, then the expectation cannot attain $p$ but I am not sure why the above argument even makes sense and what being bounded means for the expectation. Thanks!


title: Sufficient Statistic for Beta distribution

text: How can I show that the geometric mean $( \prod_{i=1}^{n} X_i )^{1/n}$ of a random sample of size $n$ from a distribution with pdf $f(x;\theta)=\theta x^{\theta-1},0<x<1,$, zero elsewhere, and $\theta>0$ is a sufficient statistic for $\theta$? By the Factorization Theorem we can see that: $$\prod_{i=1}^{n} f(x_i;\theta)=\theta^n \left[ \prod_{i=1}^n x_i \right]^{\theta-1} $$ and therefore $\prod x_i$ is a sufficient statistic for $\theta$. Now, I know that sufficient statistics are not unique and thus one does not exclude the other. This beta distribution does belong to the exponential class and as a result the above can be rewritten as: $$ \exp \{n\log\theta +(\theta-1) \sum \log x_i \}$$ but afterwards, in order to arrive at the geometric mean, there seems to be some subtle manipulation that I am missing. I would appreciate some help, thank you.


title: Symbolic matrix algebra software

text: just wondering is there a symbolic matrix algebra software like Mathematica that would solve complicated matrix expressions? Many thanks!


title: Frequency data differences test across two groups

text: I am looking into contrasting two groups with regards to frequency differences over a few time points. To make it clearer: I have for group 1, the frequency difference 5; 4; 1; 2; 3 and for group 2: 3; -3; -1; -2; 2 The negative differences are crucially important as it essentially shows the 'typical' effect is reversed by a certain magnitude. If I use a simple chi-squared test for this it will however treat my negative differences of let's say -3 just as my positive difference of 3, which is not really capturing the scale. So, I am a little unclear as to what statistical test I could perform here.


title: How adding covariance noise in Gaussian processes to prevent overfitting?

text: I am told, in Gaussian Processes, adding covariance function noise to others, say SEiso or Materns, cause a better result, since it prevents from over fitting. I appreciate if someone could put more light on it !


title: How to fix an overdispersion in a Poisson GLMM with lmer function in R?

text: I want to model counts as being dependent on two nominal variables, one continuous variable (all as fixed effects) with 3rd-order interactions and one grouping variable (as random effect). However, I have an overdispersion in outcomes (I used `lmer` function from `lme4` library). How should I manage this? I have found some solution for the problem (http://stats.stackexchange.com/a/9670/38080) but I am not able to incorporate that recommendation into my model. Here is my model:               m1<-lmer(dependent.var ~ cat.var1 * cat.var2 * contin.var + (1|group),              data = dat, family = "poisson")      Any suggestion? (I did it also like a marginal model with 'geeglm' function (library geepack), but **I would like to calculate R-squared of the model** , which is possible to obtain just from former GLMM (see Nakagawa & Schielzeth 2013; http://onlinelibrary.wiley.com/doi/10.1111/j.2041-210x.2012.00261.x/abstract).)


title: Segmented time-dependent covariate

text: I'm performing the cox regression analysis, where I have the model with 1 time-independent (treatment (yes/no)) and 1 time-dependent predictor (blood pressure level). Outcome is death (yes/no). Time until death is measured in days and the count starts after the date of treatment prescription (time=0 days) in patients. Time-independent predictor (treatment) have a value 1 (yes) or 0 (no) at time=0 days and remains the same until either death occur or time of observation is over (censored). Time-dependent predictor, in turn, may appear in the model at time 90 (on 90th day after treatment date). I compute the segmented time-dependent predictors in a way as it described in SPSS help (BP - blood pressure): COMPUTE T_COV = (T_ < BP1_time) * BP1 + (T_ >= BP1_time & T_ < BP2_time) * BP2 + (T_ >= BP2_time & T_ < BP3_time) * BP3 + (T_ >= BP3_time & T_ < BP4_time) * BP4.... **The question is:** how does SPSS deal with the time within the period from 0 to 90 days for the time-dependent predictor? Does it just ignore this time period and start to measure hazard at the day when the BP variable enters to the model?


title: What is concordance probability?

text: In one of my questions: How to measure test set error with logistic regression One of the answers mentioned "concordance probability" or c-index. Unfortunately I'm not familiar with the term and there seems to be no Wikipedia page and Google turns up nothing. It's my understanding that it relates to classification.


title: How to fit a multivariate lognormal distribution to a given dataset

text: Given a dataset X with N observations in 11 dimensions, where each variable is restricted to be >= 0, how is it possible to fit an 11-dimensional log-normal distribution to this dataset? I only found sources for fitting an univariate lognormal to data, but I didn't find anything for the multivariate case. I would be happy if this was possible in MATLAB, however python or R would also be fine.


title: Population size for Genetic algorithm - rule of thumb

text: I have a 4-dimensional function $F(a,b,c,d)$ which I need to optimize (find the minimum) Each of the 4 parameters of my function $(a,b,c,d)$ are made to vary in steps over a range, so each one can take only a certain finite number of values: $a:\{a_1,a_2, ...a_A\}\,;\, b:\{b_1,b_2, ...b_B\}\,;\, c:\{c_1,c_2, ...c_C\}\,;\, d:\{d_1,d_2, ...d_D\}$ The total number of possibles solutions is then the combination of all the values each parameter can take: $N=A*B*C*D$, and I need to find the optimal one among those $N$ solutions. Given $N$, is there a rule of thumb for the population size one should feed a genetic algorithm to ensure optimal coverage of the parameter space? Or is this decision too complicated for a simple rule of thumb to apply?


title: How do I interpret data from a 2x2 experiment with a floor effect that causes it to lack any variation in one cell?

text: I have performed a 2x2 experiment to measure the effects of two genetic factors on the occurrence and size of a rare birth defect. When an individual has the defect, we can measure its size (size is a continuous value greater than zero, with no upper bound.) Here is some size data similar to what we're seeing. When an individual doesn't have a defect it's shown as a 0 on the plot (So this is a case of severely zero-inflated data): ![http://i.stack.imgur.com/BJcjT.png](http://i.stack.imgur.com/74JYb.png) I would like to:   1. Determine whether there are differences in the incidence of defects between the four groups.   2. Determine whether there are differences in the size of the defects that do occur.  Ideally, I would like to demonstrate an interaction (or lack thereof) between the two factors. For incidence (1) I have binned these values into a very simple system such as the following:               +-----------+-----------+-----------+--------+------+--------+     | Condition |  Gene 1   |  Gene 2   | Normal | Mild | Severe |     +-----------+-----------+-----------+--------+------+--------+     | A         | Low Risk  | Low Risk  |      9 |    0 |      0 |     | B         | Low Risk  | High Risk |      5 |    5 |      1 |     | C         | High Risk | Low Risk  |      6 |   10 |      0 |     | D         | High Risk | High Risk |      1 |    5 |     16 |     +-----------+-----------+-----------+--------+------+--------+      and then I can perform Fisher's exact test between the relevant conditions: A vs B, A vs C, A vs D, B vs D, and C vs D. This works fine but it doesn't assess an interaction between the factors. For defect size (2) I initially wanted to do a two-way ANOVA (to assess interaction), but I don't think it's possible. Most individuals don't have a defect except in the most severe situation, thus I end up with zero-inflated data set. Not only that, in the Low-Low category there's no variation at all! If I drop all the zero values (which makes sense - the absence of a defect is not really the same thing as a defect of even negligible size) then I'm left with nothing at all in my low-low condition. I can make comparisons between the groups with rank-sum data, but only if I drop all the zeroes, otherwise the test is inappropriate. When I do this, I destroy 75% of my sample size (100% in the low-low case) Is there a more efficient way to look at this data?


title: Post processing random forests using regularised regression: what about bias?

text: I have been playing around with post processing the results of the random forest for regression machine learning algorithm in order to try and do better than the default mean of all trees prediction. Chapter 16 of Elements of statistical learning talks about how and why this can be done using regularised regression (in that case elastic net). On the other hand, I routinely re-calibrate the results of random forests using a logistic model, since I often see the kind of bias discussed in this question, and find that a simple linear logistic regression model tends to account for it. I have found that the calibration needs to be done on the ensemble prediction, attempting to calibrate individual trees doesn't work nearly as well. I am unsure then of how to proceed when you want to do both of these things. If you do the elastic net regression first, you are finding the optimal weighted mean/subset of the forest ensemble for the biased prediction. I see no reason to be sure that this will be the optimal combination for a re- calibrated prediction. I don't see how the elastic net analysis could be run simultaneously with re-calibration, as the calibration is non-linear. The non-linear problem of mutuality optimising the calibration parameters as well as weights and sub-set selection is doable is principle, but seem almost intractable computationally for large data sets. Can anyone suggest an efficient path for achieving both of these aims in an optimised way?


title: Inference on random graph, with an application to mobile sensors

text: I've attended a course on Machine Learning and another one in Network Analysis, and I wonder if this two topics already intersect, in particular I'm interested in the following model:   * we have a mobile network (that is, the nodes perform a random walk in some area, then sending a message within some assigned range, let's say with radius r, communicating with every other node that in that moment is in that range).   * we take the point of view of a node and would like to infer as much as we can about the whole network (under suitable assumption, for example knowing or not knowing the total number of nodes, etc.), possibly in order to perform some achieved probabilistic version of well known distributed algorithm. Are there results with respect to this scenario?


title: Information theoretic substitute for precision and recall?

text: Is there a substitute for measuring precision, recall of a classifier (binary or multi-class) to evaluate its performance using information-theoretic quantities like entropy, mutual information or normalized mutual information?


title: Johansen test loading matrix

text: I'm using the `URCA` package in `R` to test for cointegration by Johansen's method. **Can anyone tell me what the weights (loading matrix) tells me or what I would use it for?** I understand what the critical value, test statistic and eigenvectors are for but confused by the loading matrix. Here's an example I am working with:               library(urca)     da=read.table("m-bnd.txt")     aaa=da[,4]     baa=da[,5]          x=cbind(aaa,baa)     colnames(x) <- c("aaa","baa")     coeig=ca.jo(x,ecdet="const",type="eigen",K=7,spec="transitory")          summary(coeig)      The output is:               Test type: maximal eigenvalue statistic (lambda max),     without linear trend and constant in cointegration     Eigenvalues (lambda): [1] 3.963354e-02 5.378046e-03 1.853392e-19          Values of teststatistic and critical values of test:          test 10pct 5pct 1pct     r <= 1 | 3.25 7.52 9.24 12.97     r = 0 | 24.35 13.75 15.67 20.20          Eigenvectors, normalised to first column:     (These are the cointegration relations)     aaa.l1 baa.l1 constant     aaa.l1 1.00000000 1.0000000 1.0000000     baa.l1 -0.89285894 0.6561369 -0.1770799     constant 0.05492912 -13.9106270 4.6994175          Weights W:     (This is the loading matrix)          aaa.l1 baa.l1 constant     aaa.d -0.006761805 -0.003055965 1.480155e-19     baa.d 0.066353797 -0.002469648 -1.113343e-17


title: Estimating the covariance matrix in LDA

text: I was trying to derive the equations from page 109 in "elements of statistical learning" (image below) ![enter image description here](http://i.stack.imgur.com/qkaZi.png) To be honest, I am not sure how the covariance $\Sigma$ is estimated (the third bullet point in image). Can someone kindly show me and how $\hat{\Sigma}$ is derived (in particular I am not sure how $N - K$ appear in the denominator) ? Thanks UPDATE: I still cannot figure why $N-K$ (where $K$ is the number of parameters) shows up in the denominator, but I suspect it has to do with making the covariance matrix estimator unbiased. This is similar to $\frac{1}{N} \sum_i (x_i - \bar{x})^2$ is biased, while $s^2 = \frac{1}{N-1} \sum_i (x_i - \bar{x})^2$ is not. Please CORRECT ME IF I AM WRONG.


title: Correlation coefficient for non-dichotomous nominal variable and ordinal/numeric variable

text: I've already read all the pages in this site trying to find the answer to my problem but no one seems to be the right one form me... First I explain you the kind of data I'm working with... Let's say that I have an array vector with several names of city, one for each of 300 users. I also have another array vector with scores response to a survey of each user or a continuous value for each user. I would like to know if exist a correlation coefficient that compute the correlation between these two variables so, nominal and numeric/ordinal variables I've searched on the Internet and in some pages they suggest to use the contingency coefficient or Cramer's V or Lambda coefficient or Eta . For each of this measure the just say that they could be applied for such data in which we have a nominal variable and interval or numerical variable. The thing is that searching and searching, trying to understand every one of them, sometime is written or watching the examples that they are reasonable to use them if you have dichotomous nominal variable, except for Cramer's V, other time is not written any requirement for the type of data. A lot of other pages say that is right to apply regression instead, that is right, but I would just simply like to know if there is a coefficient like pearson/spearman for this kind of data. I also think that is no so properly to use Spearman Correlation coeff since the cities are not sortable. I have also built the function of Cramer'sV and Eta by myself (I'm working with Matlab) but for Eta they don't talk about any p-value to see if the coefficient is statistically significant... In the matlabWorks site there is also a nice toolbox that says to compute eta^2 but the kind of input it needs is not understandable. Is here someone that have done a test like mine? If you need more detail to understand the kind of data I'm using just ask me and I'll try to explain you better. Thanks a lot


title: Multivariate Gaussian v. PCA + Independent Gaussians for Anomaly Detection

text: In Andrew Ng's Machine Learning Coursera Class, he covers anomaly detection in multiple dimensions for both independent univariate Gaussians and multivariate Gaussians, the latter being more costly than the former. Would running independent anomaly detection after orthogonalizing the data produce the same results as a multivariate anomaly detection? Is PCA too costly for this to ever be worthwhile (assuming it works at all)?


title: Is a Chi-squared test sufficient to tell me if my process change was "good enough"?

text: **Experiment Design:** **Process** : Outputs the order in which to do the work (assumption: ordering of work items is most preferred/valuable by stakeholders) **Work items** : The process above produces a rank ordering of items. Individual teams select work items from this "prioritized" queue to do their work. I introduced a new process to perform the prioritization and we did notice some improvement in the selection of work items i.e., work items were more and more goal oriented and not as per "whims of individuals" or "star performers" only. **Question to answer** : Given 4 year experimental data I wish to know if my 'new' process was any better? **Current Approach** : For a set of work-items I compute how strongly did their selection/choice compare with the prioritization ordering of the new/old process. I compute a correlation metric (almost like an exam score in essence) which is a weighted average of the proportion of work-items completed. For example:   * Percentage of top 25% work items completed * $w_1$   * Percentage of 2nd 25% work items completed * $w_2$   * Percentage of 3rd 25% work items completed * $w_3$   * Percentage of bottom 25% work items completed * $w_4$ Weights are chosen as $w_1 = 8, w_2 = 4, w_3 = 2, w_1 = 1$ That is if highest priority tasks are selected we wish to give that a higher score. Weights are normalized on $[0,1]$. A simple coverage metric, like $10/22$ items done, would not do justice. Also 100% of the tasks can never be completed, since that makes the prioritization redundant, but it's good if that's the case, albeit rare. I create three categories based on the scores: _Strong, Medium and Weak_ if the scores are $> 0.5,[0.2, 0.5), < 0.2$, respectively. I performed a Chi-squared test as follows, based on the frequency counts of Strong/Medium/Weak correlations as above. The DV is the work-item-selection and the IV is the Prioritization Process. $H_0 = $ work-item-selection is not impacted by process change $H_1 = $ work-item-selection is impacted by process change ![enter image description here](http://i.stack.imgur.com/HtK4u.png) I conclude there's a high chance that the new process has impacted work item selection rejecting $H_0$ and accepting $H_1$. **Questions:**   1. Is this approach correct and are the inferences drawn from it valid?   2. Is $\chi^2$ the best that can be done given the above DV/IV? Are there any other tests? Would they add more value to the insight gained from $\chi^2$? **Note** : For simplicity we are considering only one IV since most of the other variables is what we controlled for keeping them almost constant (i.e., team size/composition etc.,) and we (people reviewing this experiment) have agreed that they don't impact the DV and can thus be removed from the model.


title: Multiple regression with small data sets

text: I have a dataset of project case studies for a new type of research method for Government agencies to support decision making activities. My task is to develop an estimation method based on past experience for future projects for estimation purposes. My dataset is limited to 50 cases. I have 30+ (potential) predictors recorded and one response variable (i.e. hours taken to complete the project). Not all predictors are significant, using step-wise selection techniques I'm expecting number of prediction variables is likely to be in the 5-10 variable range. Although I'm struggling to get a predictor set using the standard appraoches in tools like PASW (SPSS). I'm well aware of all the material talking about rules of thumb for sample sizes and predictor variable to case ratios. My dilemma is that it's taken close to 10 years to collect 50 cases as it is, so it's about as good as it will get. My question is what should I do to get the most out of this small sample set? That is any good references for dealing with small smaple sets? Changes in p-value significance? Changes to step-wise selection approaches? Use of transforms such as centre-ing or log? Any advice is appreciated.


title: What is the difference between using aov() and lme() in analyzing a longitudinal dataset?

text: Can anyone tell me the difference between using `aov()` and `lme()` for analyzing longitudinal data and how to interpret results from these two methods? Below, I analyze the same dataset using `aov()` and `lme()` and got 2 different results. With `aov()`, I got a significant result in the time by treatment interaction, but fitting a linear mixed model, time by treatment interaction is insignificant.               > UOP.kg.aov <- aov(UOP.kg~time*treat+Error(id), raw3.42)     > summary(UOP.kg.aov)          Error: id               Df  Sum Sq Mean Sq F value Pr(>F)     treat      1   0.142  0.1421  0.0377 0.8471     Residuals 39 147.129  3.7725                         Error: Within                 Df  Sum Sq Mean Sq  F value  Pr(>F)         time         1 194.087 194.087 534.3542 < 2e-16 ***     time:treat   1   2.077   2.077   5.7197 0.01792 *       Residuals  162  58.841   0.363                          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          > UOP.kg.lme <- lme(UOP.kg~time*treat, random=list(id=pdDiag(~time)),                          na.action=na.omit, raw3.42)     > summary(UOP.kg.lme)     Linear mixed-effects model fit by REML      Data: raw3.42             AIC      BIC    logLik       225.7806 248.9037 -105.8903          Random effects:      Formula: ~time | id      Structure: Diagonal             (Intercept)      time  Residual     StdDev:   0.6817425 0.5121545 0.1780466          Fixed effects: UOP.kg ~ time + treat + time:treat                       Value Std.Error  DF   t-value p-value     (Intercept)  0.5901420 0.1480515 162  3.986059  0.0001     time         0.8623864 0.1104533 162  7.807701  0.0000     treat       -0.2144487 0.2174843  39 -0.986042  0.3302     time:treat   0.1979578 0.1622534 162  1.220053  0.2242      Correlation:                 (Intr) time   treat      time       -0.023                   treat      -0.681  0.016            time:treat  0.016 -0.681 -0.023          Standardized Within-Group Residuals:              Min           Q1          Med           Q3          Max      -3.198315285 -0.384858426  0.002705899  0.404637305  2.049705655           Number of Observations: 205     Number of Groups: 41


title: Phi coefficient and odds ratios

text: I have something strange going on. My phi coefficient of two binary variables is .07 while my odds ratio of the same two binary variables is 1.80. How is this possible?


title: Determining statistical significance of difference between two binomial distributions

text: Basically, I have two distributions with the same mean value and same number of observations. One distribution is experimental data, the other being pre-generated. I need to be able to determine the probability these two distributions are of the same set, in order to determine if statistical significance is present, however so far I have been unable to find anything of the sort. As far as I can tell, the Kolmogorov-Smirnov test doesn't fit my needs, but I would be glad to be shown otherwise. One step further would be that I plan to implement this in Java, however I should be able to handle that end myself. Thankyou in advance!


title: Averaging correlation values

text: Let's say I test how variable `Y` depends on variable `X` under different experimental conditions and obtain the following graph: ![enter image description here](http://i.stack.imgur.com/hpihy.png) The dash lines in the graph above represent linear regression for each data series (experimental setup) and the numbers in the legend denote the Pearson correlation of each data series. I would like to calculate the "average correlation" (or "mean correlation") between `X` and `Y`. May I simply average the `r` values? What about the "average determination criterion", $R^2$? Should I calculate the average `r` and than take the square of that value or should I compute the average of individual $R^2$'s?